Mercurial > hg > octave-nkf
annotate scripts/polynomial/polyfit.m @ 9290:40ba43a4745f
fix too optimistic polyfit test
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Wed, 03 Jun 2009 07:38:17 +0200 |
parents | 16f53d29049f |
children | a8ce6bdecce5 |
rev | line source |
---|---|
7017 | 1 ## Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, |
9245 | 2 ## 2007, 2008, 2009 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 -*- |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
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 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
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 ## | |
7501 | 27 ## The second output is a structure containing the following fields: |
3426 | 28 ## |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
29 ## @table @samp |
4491 | 30 ## @item R |
7501 | 31 ## Triangular factor R from the QR decomposition. |
4491 | 32 ## @item X |
7501 | 33 ## The Vandermonde matrix used to compute the polynomial coefficients. |
4491 | 34 ## @item df |
7501 | 35 ## The degrees of freedom. |
4491 | 36 ## @item normr |
7501 | 37 ## The norm of the residuals. |
4491 | 38 ## @item yf |
7501 | 39 ## The values of the polynomial for each value of @var{x}. |
4491 | 40 ## @end table |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
41 ## |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
42 ## The second output may be used by @code{polyval} to calculate the |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
43 ## statistical error limits of the predicted values. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
44 ## |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
45 ## When the third output, @var{mu}, is present the |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
46 ## coefficients, @var{p}, are associated with a polynomial in |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
47 ## @var{xhat} = (@var{x}-@var{mu}(1))/@var{mu}(2). |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
48 ## Where @var{mu}(1) = mean (@var{x}), and @var{mu}(2) = std (@var{x}). |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
49 ## This linear transformation of @var{x} improves the numerical |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
50 ## stability of the fit. |
8286
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
7547
diff
changeset
|
51 ## @seealso{polyval, 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 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
60 if (nargin < 3 || nargin > 4) |
6046 | 61 print_usage (); |
2261 | 62 endif |
2325 | 63 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
64 if (nargout > 2) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
65 ## Normalized the x values. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
66 mu = [mean(x), std(x)]; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
67 x = (x - mu(1)) / mu(2); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
68 endif |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
69 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
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 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
80 ## Reshape x & y into column vectors. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
81 l = numel (x); |
9138 | 82 x = x(:); |
83 y = y(:); | |
2325 | 84 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
85 ## Construct the Vandermonde matrix. |
9138 | 86 v = vander (x, n+1); |
2261 | 87 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
88 ## Solve by QR decomposition. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
89 [q, r, k] = qr (v, 0); |
9138 | 90 p = r \ (q' * y); |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
91 p(k) = p; |
3091 | 92 |
4491 | 93 if (nargout > 1) |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
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 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
102 s.R = r; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
103 s.X = v; |
4491 | 104 s.df = l - n - 1; |
105 s.normr = norm (yf - y); | |
2261 | 106 endif |
107 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
108 ## Return a row vector. |
5395 | 109 p = p.'; |
110 | |
2261 | 111 endfunction |
7411 | 112 |
113 %!test | |
114 %! x = [-2, -1, 0, 1, 2]; | |
115 %! assert(all (all (abs (polyfit (x, x.^2+x+1, 2) - [1, 1, 1]) < sqrt (eps)))); | |
116 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
117 %!error(polyfit ([1, 2; 3, 4], [1, 2, 3, 4], 2)) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
118 |
7411 | 119 %!test |
120 %! x = [-2, -1, 0, 1, 2]; | |
121 %! assert(all (all (abs (polyfit (x, x.^2+x+1, 3) - [0, 1, 1, 1]) < sqrt (eps)))); | |
122 | |
123 %!test | |
124 %! x = [-2, -1, 0, 1, 2]; | |
125 %! fail("polyfit (x, x.^2+x+1)"); | |
126 | |
127 %!test | |
128 %! x = [-2, -1, 0, 1, 2]; | |
129 %! fail("polyfit (x, x.^2+x+1, [])"); | |
130 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
131 ## Test difficult case where scaling is really needed. This example |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
132 ## demonstrates the rather poor result which occurs when the dependent |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
133 ## variable is not normalized properly. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
134 ## Also check the usage of 2nd & 3rd output arguments. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
135 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
136 %! x = [ -1196.4, -1195.2, -1194, -1192.8, -1191.6, -1190.4, -1189.2, -1188, \ |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
137 %! -1186.8, -1185.6, -1184.4, -1183.2, -1182]; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
138 %! y = [ 315571.7086, 315575.9618, 315579.4195, 315582.6206, 315585.4966, \ |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
139 %! 315588.3172, 315590.9326, 315593.5934, 315596.0455, 315598.4201, \ |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
140 %! 315600.7143, 315602.9508, 315605.1765 ]; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
141 %! [p1, s1] = polyfit (x, y, 10); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
142 %! [p2, s2, mu] = polyfit (x, y, 10); |
9290
40ba43a4745f
fix too optimistic polyfit test
Jaroslav Hajek <highegg@gmail.com>
parents:
9245
diff
changeset
|
143 %! assert (s2.normr < s1.normr) |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
144 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
145 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
146 %! x = 1:4; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
147 %! p0 = [1i, 0, 2i, 4]; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
148 %! y0 = polyval (p0, x); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
149 %! p = polyfit (x, y0, numel(p0)-1); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
150 %! assert (p, p0, 1000*eps) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
151 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
152 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
153 %! x = 1000 + (-5:5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
154 %! xn = (x - mean (x)) / std (x); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
155 %! pn = ones (1,5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
156 %! y = polyval (pn, xn); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
157 %! [p, s, mu] = polyfit (x, y, numel(pn)-1); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
158 %! [p2, s2] = polyfit (x, y, numel(pn)-1); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
159 %! assert (p, pn, s.normr) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
160 %! assert (s.yf, y, s.normr) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
161 %! assert (mu, [mean(x), std(x)]) |
7547 | 162 %! assert (s.normr/s2.normr < sqrt(eps)) |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
163 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
164 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
165 %! x = [1, 2, 3; 4, 5, 6]; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
166 %! y = [0, 0, 1; 1, 0, 0]; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
167 %! p = polyfit (x, y, 5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
168 %! expected = [0, 1, -14, 65, -112, 60]/12; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
169 %! assert (p, expected, sqrt(eps)) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
170 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
171 |