Mercurial > hg > octave-lyh
annotate scripts/polynomial/polyfit.m @ 14200:64d9f33313cc stable rc-3-6-0-1
3.6.0-rc1 release candidate
* configure.ac (AC_INIT): Version is now 3.6.0-rc1.
(OCTAVE_RELEASE_DATE): Now 2012-01-12.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 12 Jan 2012 14:31:50 -0500 |
parents | 72c96de7a403 |
children | f3d52523cde1 |
rev | line source |
---|---|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
14110
diff
changeset
|
1 ## Copyright (C) 1996-2012 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 |
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
14062
diff
changeset
|
24 ## @var{n} that minimizes the least-squares-error of the fit to the points |
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
14062
diff
changeset
|
25 ## @code{[@var{x}, @var{y}]}. |
3418 | 26 ## |
4491 | 27 ## The polynomial coefficients are returned in a row vector. |
28 ## | |
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
14062
diff
changeset
|
29 ## The optional output @var{s} is a structure containing the following fields: |
3426 | 30 ## |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
31 ## @table @samp |
4491 | 32 ## @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
|
33 ## Triangular factor R from the QR@tie{}decomposition. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
34 ## |
4491 | 35 ## @item X |
7501 | 36 ## 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
|
37 ## |
4491 | 38 ## @item df |
7501 | 39 ## The degrees of freedom. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
40 ## |
4491 | 41 ## @item normr |
7501 | 42 ## The norm of the residuals. |
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
43 ## |
4491 | 44 ## @item yf |
7501 | 45 ## The values of the polynomial for each value of @var{x}. |
4491 | 46 ## @end table |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
47 ## |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
48 ## 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
|
49 ## 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
|
50 ## |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
51 ## 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
|
52 ## 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
|
53 ## @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
|
54 ## 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
|
55 ## 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
|
56 ## stability of the fit. |
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
14062
diff
changeset
|
57 ## @seealso{polyval, polyaffine, roots, vander, zscore} |
3368 | 58 ## @end deftypefn |
2311 | 59 |
5428 | 60 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
2312 | 61 ## Created: 13 December 1994 |
62 ## Adapted-By: jwe | |
63 | |
4491 | 64 function [p, s, mu] = polyfit (x, y, n) |
2325 | 65 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
66 if (nargin < 3 || nargin > 4) |
6046 | 67 print_usage (); |
2261 | 68 endif |
2325 | 69 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
70 if (nargout > 2) |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
71 ## Normalized the x values. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
72 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
|
73 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
|
74 endif |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
75 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
76 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
|
77 error ("polyfit: X and Y must be vectors of the same size"); |
2261 | 78 endif |
2325 | 79 |
13279
984359717d71
Use common code idiom for checking whether a double value is an integer.
Rik <octave@nomad.inbox5.com>
parents:
11593
diff
changeset
|
80 if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == fix (n))) |
14062
5b49cafe0599
Use non-negative, non-positive with hyphens in error messages.
Rik <octave@nomad.inbox5.com>
parents:
13279
diff
changeset
|
81 error ("polyfit: N must be a non-negative integer"); |
2261 | 82 endif |
2325 | 83 |
14110
b8d9530e940e
undo non-documentation changes in changeset 614505385171
John W. Eaton <jwe@octave.org>
parents:
14104
diff
changeset
|
84 y_is_row_vector = (rows (y) == 1); |
b8d9530e940e
undo non-documentation changes in changeset 614505385171
John W. Eaton <jwe@octave.org>
parents:
14104
diff
changeset
|
85 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
86 ## 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
|
87 l = numel (x); |
9138 | 88 x = x(:); |
89 y = y(:); | |
2325 | 90 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
91 ## Construct the Vandermonde matrix. |
9138 | 92 v = vander (x, n+1); |
2261 | 93 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
94 ## Solve by QR decomposition. |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
95 [q, r, k] = qr (v, 0); |
9138 | 96 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
|
97 p(k) = p; |
3091 | 98 |
4491 | 99 if (nargout > 1) |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
100 yf = v*p; |
3105 | 101 |
14110
b8d9530e940e
undo non-documentation changes in changeset 614505385171
John W. Eaton <jwe@octave.org>
parents:
14104
diff
changeset
|
102 if (y_is_row_vector) |
4491 | 103 s.yf = yf.'; |
104 else | |
105 s.yf = yf; | |
3105 | 106 endif |
3091 | 107 |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
108 s.R = r; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
109 s.X = v; |
4491 | 110 s.df = l - n - 1; |
111 s.normr = norm (yf - y); | |
2261 | 112 endif |
113 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
114 ## Return a row vector. |
5395 | 115 p = p.'; |
116 | |
2261 | 117 endfunction |
7411 | 118 |
119 %!test | |
120 %! x = [-2, -1, 0, 1, 2]; | |
121 %! assert(all (all (abs (polyfit (x, x.^2+x+1, 2) - [1, 1, 1]) < sqrt (eps)))); | |
122 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
123 %!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
|
124 |
7411 | 125 %!test |
126 %! x = [-2, -1, 0, 1, 2]; | |
127 %! assert(all (all (abs (polyfit (x, x.^2+x+1, 3) - [0, 1, 1, 1]) < sqrt (eps)))); | |
128 | |
129 %!test | |
130 %! x = [-2, -1, 0, 1, 2]; | |
131 %! fail("polyfit (x, x.^2+x+1)"); | |
132 | |
133 %!test | |
134 %! x = [-2, -1, 0, 1, 2]; | |
135 %! fail("polyfit (x, x.^2+x+1, [])"); | |
136 | |
7500
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
137 ## 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
|
138 ## 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
|
139 ## 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
|
140 ## 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
|
141 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
142 %! 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
|
143 %! -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
|
144 %! 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
|
145 %! 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
|
146 %! 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
|
147 %! [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
|
148 %! [p2, s2, mu] = polyfit (x, y, 10); |
9290
40ba43a4745f
fix too optimistic polyfit test
Jaroslav Hajek <highegg@gmail.com>
parents:
9245
diff
changeset
|
149 %! 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
|
150 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
151 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
152 %! x = 1:4; |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
153 %! 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
|
154 %! y0 = polyval (p0, x); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
155 %! 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
|
156 %! 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
|
157 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
158 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
159 %! x = 1000 + (-5:5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
160 %! 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
|
161 %! pn = ones (1,5); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
162 %! y = polyval (pn, xn); |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
163 %! [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
|
164 %! [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
|
165 %! 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
|
166 %! 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
|
167 %! assert (mu, [mean(x), std(x)]) |
7547 | 168 %! 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
|
169 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
170 %!test |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
171 %! 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
|
172 %! 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
|
173 %! 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
|
174 %! 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
|
175 %! 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
|
176 |
2df882e69f13
use QR decomposition and normalization for polyfit; normalization for polyval
Ben Abbott <bpabbott@mac.com>
parents:
7411
diff
changeset
|
177 |