Mercurial > hg > octave-lyh
annotate scripts/polynomial/polyfit.m @ 13279:984359717d71
Use common code idiom for checking whether a double value is an integer.
* num2str.m, rotdim.m, get_first_help_sentence.m, ind2rgb.m,
commutation_matrix.m, figure.m, legend.m, polyfit.m, bartlett.m, blackman.m,
detrend.m, hamming.m, hanning.m, factorial.m, mode.m, skewness.m, statistics.m,
mcnemar_test.m: Use idiom 'x == fix (x)' to test for integerness.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Wed, 05 Oct 2011 13:10:10 -0700 |
parents | 1577c6f80926 |
children | 5b49cafe0599 |
rev | line source |
---|---|
11523 | 1 ## Copyright (C) 1996-2011 John W. Eaton |
2313 | 2 ## |
3 ## This file is part of Octave. | |
4 ## | |
5 ## Octave is free software; you can redistribute it and/or modify it | |
6 ## under the terms of the GNU General Public License as published by | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
2313 | 9 ## |
10 ## Octave is distributed in the hope that it will be useful, but | |
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 ## General Public License for more details. | |
14 ## | |
15 ## You should have received a copy of the GNU General Public License | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
2303 | 18 |
3368 | 19 ## -*- texinfo -*- |
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
9290
diff
changeset
|
20 ## @deftypefn {Function File} {@var{p} =} polyfit (@var{x}, @var{y}, @var{n}) |
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
9290
diff
changeset
|
21 ## @deftypefnx {Function File} {[@var{p}, @var{s}] =} polyfit (@var{x}, @var{y}, @var{n}) |
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
9290
diff
changeset
|
22 ## @deftypefnx {Function File} {[@var{p}, @var{s}, @var{mu}] =} polyfit (@var{x}, @var{y}, @var{n}) |
3368 | 23 ## 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
|
24 ## @var{n} that minimizes the least-squares-error of the fit. |
3418 | 25 ## |
4491 | 26 ## The polynomial coefficients are returned in a row vector. |
27 ## | |
7501 | 28 ## The second output is a structure containing the following fields: |
3426 | 29 ## |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
30 ## @table @samp |
4491 | 31 ## @item R |
11593
1577c6f80926
Use non-breaking spaces between certain adjectives and their nouns in docstrings.
Rik <octave@nomad.inbox5.com>
parents:
11587
diff
changeset
|
32 ## Triangular factor R from the QR@tie{}decomposition. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
33 ## |
4491 | 34 ## @item X |
7501 | 35 ## The Vandermonde matrix used to compute the polynomial coefficients. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
36 ## |
4491 | 37 ## @item df |
7501 | 38 ## The degrees of freedom. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
39 ## |
4491 | 40 ## @item normr |
7501 | 41 ## The norm of the residuals. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
42 ## |
4491 | 43 ## @item yf |
7501 | 44 ## The values of the polynomial for each value of @var{x}. |
4491 | 45 ## @end table |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
46 ## |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
47 ## The second output may be used by @code{polyval} to calculate the |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
48 ## 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
|
49 ## |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
50 ## When the third output, @var{mu}, is present the |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
51 ## 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
|
52 ## @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
|
53 ## 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
|
54 ## 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
|
55 ## stability of the fit. |
8286
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
7547
diff
changeset
|
56 ## @seealso{polyval, residue} |
3368 | 57 ## @end deftypefn |
2311 | 58 |
5428 | 59 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
2312 | 60 ## Created: 13 December 1994 |
61 ## Adapted-By: jwe | |
62 | |
4491 | 63 function [p, s, mu] = polyfit (x, y, n) |
2325 | 64 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
65 if (nargin < 3 || nargin > 4) |
6046 | 66 print_usage (); |
2261 | 67 endif |
2325 | 68 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
69 if (nargout > 2) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
70 ## Normalized the x values. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
71 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
|
72 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
|
73 endif |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
74 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
75 if (! size_equal (x, y)) |
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
10821
diff
changeset
|
76 error ("polyfit: X and Y must be vectors of the same size"); |
2261 | 77 endif |
2325 | 78 |
13279
984359717d71
Use common code idiom for checking whether a double value is an integer.
Rik <octave@nomad.inbox5.com>
parents:
11593
diff
changeset
|
79 if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == fix (n))) |
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
10821
diff
changeset
|
80 error ("polyfit: N must be a nonnegative integer"); |
2261 | 81 endif |
2325 | 82 |
3091 | 83 y_is_row_vector = (rows (y) == 1); |
84 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
85 ## 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
|
86 l = numel (x); |
9138 | 87 x = x(:); |
88 y = y(:); | |
2325 | 89 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
90 ## Construct the Vandermonde matrix. |
9138 | 91 v = vander (x, n+1); |
2261 | 92 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
93 ## Solve by QR decomposition. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
94 [q, r, k] = qr (v, 0); |
9138 | 95 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
|
96 p(k) = p; |
3091 | 97 |
4491 | 98 if (nargout > 1) |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
99 yf = v*p; |
3105 | 100 |
101 if (y_is_row_vector) | |
4491 | 102 s.yf = yf.'; |
103 else | |
104 s.yf = yf; | |
3105 | 105 endif |
3091 | 106 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
107 s.R = r; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
108 s.X = v; |
4491 | 109 s.df = l - n - 1; |
110 s.normr = norm (yf - y); | |
2261 | 111 endif |
112 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
113 ## Return a row vector. |
5395 | 114 p = p.'; |
115 | |
2261 | 116 endfunction |
7411 | 117 |
118 %!test | |
119 %! x = [-2, -1, 0, 1, 2]; | |
120 %! assert(all (all (abs (polyfit (x, x.^2+x+1, 2) - [1, 1, 1]) < sqrt (eps)))); | |
121 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
122 %!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
|
123 |
7411 | 124 %!test |
125 %! x = [-2, -1, 0, 1, 2]; | |
126 %! assert(all (all (abs (polyfit (x, x.^2+x+1, 3) - [0, 1, 1, 1]) < sqrt (eps)))); | |
127 | |
128 %!test | |
129 %! x = [-2, -1, 0, 1, 2]; | |
130 %! fail("polyfit (x, x.^2+x+1)"); | |
131 | |
132 %!test | |
133 %! x = [-2, -1, 0, 1, 2]; | |
134 %! fail("polyfit (x, x.^2+x+1, [])"); | |
135 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
136 ## 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
|
137 ## 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
|
138 ## 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
|
139 ## 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
|
140 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
141 %! 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
|
142 %! -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
|
143 %! 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
|
144 %! 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
|
145 %! 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
|
146 %! [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
|
147 %! [p2, s2, mu] = polyfit (x, y, 10); |
9290
40ba43a4745f
fix too optimistic polyfit test
Jaroslav Hajek <highegg@gmail.com>
parents:
9245
diff
changeset
|
148 %! 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
|
149 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
150 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
151 %! x = 1:4; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
152 %! 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
|
153 %! y0 = polyval (p0, x); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
154 %! 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
|
155 %! 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
|
156 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
157 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
158 %! x = 1000 + (-5:5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
159 %! 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
|
160 %! pn = ones (1,5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
161 %! y = polyval (pn, xn); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
162 %! [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
|
163 %! [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
|
164 %! 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
|
165 %! 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
|
166 %! assert (mu, [mean(x), std(x)]) |
7547 | 167 %! 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
|
168 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
169 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
170 %! 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
|
171 %! 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
|
172 %! 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
|
173 %! 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
|
174 %! 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
|
175 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
176 |