Mercurial > octave
annotate scripts/polynomial/polyfit.m @ 11917:deb777a926ee release-3-0-x
style fixes
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 12 Jan 2009 12:13:21 +0100 |
parents | 377d908f7e40 |
children | f2af2233ce7f |
rev | line source |
---|---|
7017 | 1 ## Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, |
2 ## 2007 John W. Eaton | |
2313 | 3 ## |
4 ## This file is part of Octave. | |
5 ## | |
6 ## Octave is free software; you can redistribute it and/or modify it | |
7 ## under the terms of the GNU General Public License as published by | |
7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
9 ## your option) any later version. | |
2313 | 10 ## |
11 ## Octave is distributed in the hope that it will be useful, but | |
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 ## General Public License for more details. | |
15 ## | |
16 ## You should have received a copy of the GNU General Public License | |
7016 | 17 ## along with Octave; see the file COPYING. If not, see |
18 ## <http://www.gnu.org/licenses/>. | |
2303 | 19 |
3368 | 20 ## -*- texinfo -*- |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
21 ## @deftypefn {Function File} {[@var{p}, @var{s}, @var{mu}] =} polyfit (@var{x}, @var{y}, @var{n}) |
3368 | 22 ## Return the coefficients of a polynomial @var{p}(@var{x}) of degree |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
23 ## @var{n} that minimizes the least-squares-error of the fit. |
3418 | 24 ## |
4491 | 25 ## The polynomial coefficients are returned in a row vector. |
26 ## | |
11917 | 27 ## The second output is a structure containing the following fields: |
3426 | 28 ## |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
29 ## @table @samp |
4491 | 30 ## @item R |
11917 | 31 ## Triangular factor R from the QR decomposition. |
4491 | 32 ## @item X |
11917 | 33 ## The Vandermonde matrix used to compute the polynomial coefficients. |
4491 | 34 ## @item df |
11917 | 35 ## The degrees of freedom. |
4491 | 36 ## @item normr |
11917 | 37 ## The norm of the residuals. |
4491 | 38 ## @item yf |
11917 | 39 ## The values of the polynomial for each value of @var{x}. |
4491 | 40 ## @end table |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
41 ## |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
42 ## The second output may be used by @code{polyval} to calculate the |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
43 ## statistical error limits of the predicted values. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
44 ## |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
45 ## When the third output, @var{mu}, is present the |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
46 ## coefficients, @var{p}, are associated with a polynomial in |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
47 ## @var{xhat} = (@var{x}-@var{mu}(1))/@var{mu}(2). |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
48 ## Where @var{mu}(1) = mean (@var{x}), and @var{mu}(2) = std (@var{x}). |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
49 ## This linear transformation of @var{x} improves the numerical |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
50 ## stability of the fit. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
51 ## @seealso{polyval, polyconf, residue} |
3368 | 52 ## @end deftypefn |
2311 | 53 |
5428 | 54 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
2312 | 55 ## Created: 13 December 1994 |
56 ## Adapted-By: jwe | |
57 | |
4491 | 58 function [p, s, mu] = polyfit (x, y, n) |
2325 | 59 |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
60 if (nargin < 3 || nargin > 4) |
6046 | 61 print_usage (); |
2261 | 62 endif |
2325 | 63 |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
64 if (nargout > 2) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
65 ## Normalized the x values. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
66 mu = [mean(x), std(x)]; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
67 x = (x - mu(1)) / mu(2); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
68 endif |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
69 |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
70 if (! size_equal (x, y)) |
2261 | 71 error ("polyfit: x and y must be vectors of the same size"); |
72 endif | |
2325 | 73 |
4030 | 74 if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == round (n))) |
2261 | 75 error ("polyfit: n must be a nonnegative integer"); |
76 endif | |
2325 | 77 |
3091 | 78 y_is_row_vector = (rows (y) == 1); |
79 | |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
80 ## Reshape x & y into column vectors. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
81 l = numel (x); |
2261 | 82 x = reshape (x, l, 1); |
83 y = reshape (y, l, 1); | |
2325 | 84 |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
85 ## Construct the Vandermonde matrix. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
86 v = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0)); |
2261 | 87 |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
88 ## Solve by QR decomposition. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
89 [q, r, k] = qr (v, 0); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
90 p = r \ (y' * q)'; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
91 p(k) = p; |
3091 | 92 |
4491 | 93 if (nargout > 1) |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
94 yf = v*p; |
3105 | 95 |
96 if (y_is_row_vector) | |
4491 | 97 s.yf = yf.'; |
98 else | |
99 s.yf = yf; | |
3105 | 100 endif |
3091 | 101 |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
102 s.R = r; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
103 s.X = v; |
4491 | 104 s.df = l - n - 1; |
105 s.normr = norm (yf - y); | |
2261 | 106 endif |
107 | |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
108 ## Return a row vector. |
5395 | 109 p = p.'; |
110 | |
11916
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
111 ## Test difficult case where scaling is really needed. This example |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
112 ## demonstrates the rather poor result which occurs when the dependent |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
113 ## variable is not normalized properly. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
114 ## Also check the usage of 2nd & 3rd output arguments. |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
115 %!test |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
116 %! x = [ -1196.4, -1195.2, -1194, -1192.8, -1191.6, -1190.4, -1189.2, -1188, \ |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
117 %! -1186.8, -1185.6, -1184.4, -1183.2, -1182]; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
118 %! y = [ 315571.7086, 315575.9618, 315579.4195, 315582.6206, 315585.4966, \ |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
119 %! 315588.3172, 315590.9326, 315593.5934, 315596.0455, 315598.4201, \ |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
120 %! 315600.7143, 315602.9508, 315605.1765 ]; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
121 %! [p1, s1] = polyfit (x, y, 10); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
122 %! [p2, s2, mu] = polyfit (x, y, 10); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
123 %! assert (s1.normr, 0.11264, 0.1) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
124 %! assert (s2.normr < s1.normr) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
125 |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
126 %!test |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
127 %! x = 1:4; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
128 %! p0 = [1i, 0, 2i, 4]; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
129 %! y0 = polyval (p0, x); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
130 %! p = polyfit (x, y0, numel(p0)-1); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
131 %! assert (p, p0, 1000*eps) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
132 |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
133 %!test |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
134 %! x = 1000 + (-5:5); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
135 %! xn = (x - mean (x)) / std (x); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
136 %! pn = ones (1,5); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
137 %! y = polyval (pn, xn); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
138 %! [p, s, mu] = polyfit (x, y, numel(pn)-1); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
139 %! [p2, s2] = polyfit (x, y, numel(pn)-1); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
140 %! assert (p, pn, s.normr) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
141 %! assert (s.yf, y, s.normr) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
142 %! assert (mu, [mean(x), std(x)]) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
143 %! assert (s.normr/s2.normr < 1e-9) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
144 |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
145 %!test |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
146 %! x = [1, 2, 3; 4, 5, 6]; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
147 %! y = [0, 0, 1; 1, 0, 0]; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
148 %! p = polyfit (x, y, 5); |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
149 %! expected = [0, 1, -14, 65, -112, 60]/12; |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
150 %! assert (p, expected, sqrt(eps)) |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
151 |
377d908f7e40
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7017
diff
changeset
|
152 |
2261 | 153 endfunction |