Mercurial > octave
view scripts/polynomial/polyfit.m @ 31221:f5755dbacd8d
maint: merge stable to default
author | Pantxo Diribarne <pantxo.diribarne@gmail.com> |
---|---|
date | Wed, 31 Aug 2022 22:04:02 +0200 |
parents | 796f54d4ddbf |
children | 597f3ee61a48 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 1996-2022 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{p} =} polyfit (@var{x}, @var{y}, @var{n}) ## @deftypefnx {} {[@var{p}, @var{s}] =} polyfit (@var{x}, @var{y}, @var{n}) ## @deftypefnx {} {[@var{p}, @var{s}, @var{mu}] =} polyfit (@var{x}, @var{y}, @var{n}) ## Return the coefficients of a polynomial @var{p}(@var{x}) of degree @var{n} ## that minimizes the least-squares-error of the fit to the points ## @code{[@var{x}(:), @var{y}(:)]}. ## ## @var{n} is typically an integer @geq{} 0 specifying the degree of the ## approximating polynomial. If @var{n} is a logical vector, it is used as a ## mask to selectively force the corresponding polynomial coefficients to be ## used or ignored. ## ## The polynomial coefficients are returned in the row vector @var{p}. The ## output @var{p} may be directly used with @code{polyval} to estimate values ## using the fitted polynomial. ## ## The optional output @var{s} is a structure containing the following fields: ## ## @table @samp ## ## @item yf ## The values of the polynomial for each value of @var{x}. ## ## @item X ## The @nospell{Vandermonde} matrix used to compute the polynomial ## coefficients. ## ## @item R ## Triangular factor R from the QR@tie{}decomposition. ## ## @item C ## The unscaled covariance matrix, formally equal to the inverse of ## @var{x'}*@var{x}, but computed in a way minimizing roundoff error ## propagation. ## ## @item df ## The degrees of freedom. ## ## @item normr ## The norm of the residuals. ## @end table ## ## The second output may be used by @code{polyval} to calculate the statistical ## error limits of the predicted values. In particular, the standard deviation ## of @var{p} coefficients is given by ## ## @code{sqrt (diag (@var{s.C})/@var{s.df}) * @var{s.normr}}. ## ## When the third output, @var{mu}, is present the original data is centered ## and scaled which can improve the numerical stability of the fit. The ## coefficients @var{p} are associated with a polynomial in ## ## @code{@var{xhat} = (@var{x} - @var{mu}(1)) / @var{mu}(2)} @* ## where @var{mu}(1) = mean (@var{x}), and @var{mu}(2) = std (@var{x}). ## ## Example 1 : logical @var{n} and integer @var{n} ## ## @example ## @group ## f = @@(x) x.^2 + 5; # data-generating function ## x = 0:5; ## y = f (x); ## ## Fit data to polynomial A*x^3 + B*x^1 ## p = polyfit (x, y, logical ([1, 0, 1, 0])) ## @result{} p = [ 0.0680, 0, 4.2444, 0 ] ## ## Fit data to polynomial using all terms up to x^3 ## p = polyfit (x, y, 3) ## @result{} p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ] ## @end group ## @end example ## ## @seealso{polyval, polyaffine, roots, vander, zscore} ## @end deftypefn function [p, s, mu] = polyfit (x, y, n) if (nargin < 3) print_usage (); endif y_is_row_vector = isrow (y); ## Reshape x & y into column vectors. x = x(:); y = y(:); nx = numel (x); ny = numel (y); if (nx != ny) error ("polyfit: X and Y must have the same number of points"); endif if (nargout > 2) ## Center and scale the x values. mu = [mean(x), std(x)]; x = (x - mu(1)) / mu(2); endif ## n is the polynomial degree (an input, or deduced from the polymask size) ## m is the effective number of coefficients. if (islogical (n)) polymask = n(:).'; # force to row vector n = numel (polymask) - 1; m = sum (polymask) - 1; pad_output = true; else if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == fix (n))) error ("polyfit: N must be a non-negative integer"); endif polymask = true (1, n+1); m = n; pad_output = false; endif if (m >= nx) warning ("polyfit: degree of polynomial N is >= number of data points; solution is not unique"); m = nx; pad_output = true; ## Keep the lowest m entries in polymask idx = find (polymask); idx((end-m+1):end) = []; polymask(idx) = false; endif ## Construct the Vandermonde matrix. X = vander (x, n+1); v = X(:, polymask); ## Solve by QR decomposition. [q, r, k] = qr (v, 0); p = r \ (q' * y); p(k) = p; if (isargout (2)) yf = v*p; if (y_is_row_vector) s.yf = yf.'; else s.yf = yf; endif s.X = X; ## r.'*r is positive definite if matrix v is of full rank. Invert it by ## cholinv to avoid taking the square root of squared quantities. ## If cholinv fails, then v is rank deficient and not invertible. try C = cholinv (r.'*r)(k, k); catch C = NaN (m, m); end_try_catch if (pad_output) s.X(:, ! polymask) = 0; s.R = zeros (rows (r), n+1); s.R(:, polymask) = r; s.C = zeros (rows (C), n+1); s.C(:, polymask) = C; else s.R = r; s.C = C; endif s.df = max (0, nx - m - 1); s.normr = norm (yf - y); endif if (pad_output) ## Zero pad output q = p; p = zeros (n+1, 1); p(polymask) = q; endif p = p.'; # Return a row vector. endfunction %!shared x %! x = [-2, -1, 0, 1, 2]; %!assert (polyfit (x, 3*x.^2 + 2*x + 1, 2), [3, 2, 1], 10*eps) %!assert (polyfit (x, 3*x.^2 + 2*x + 1, logical ([1 1 1])), [3, 2, 1], 10*eps) %!assert (polyfit (x, x.^2 + 2*x + 3, 3), [0, 1, 2, 3], 10*eps) %!assert (polyfit (x, x.^2 + 2*x + 3, logical ([0 1 1 1])), [0 1 2 3], 10*eps) ## Test logical input N %!test %! x = [0:5]; %! y = 3*x.^3 + 2*x.^2 + 4; %! [p, s] = polyfit (x, y, logical ([1, 0, 1, 1])); %! assert (p(2), 0); %! assert (all (p([1, 3, 4]))); %! assert (s.df, 3); ## Test difficult case where scaling is really needed. This example ## demonstrates the rather poor result which occurs when the dependent ## variable is not normalized properly. ## Also check the usage of 2nd & 3rd output arguments. %!test %! warning ("off", "Octave:nearly-singular-matrix", "local"); %! x = [ -1196.4, -1195.2, -1194, -1192.8, -1191.6, -1190.4, -1189.2, -1188, ... %! -1186.8, -1185.6, -1184.4, -1183.2, -1182]; %! y = [ 315571.7086, 315575.9618, 315579.4195, 315582.6206, 315585.4966, ... %! 315588.3172, 315590.9326, 315593.5934, 315596.0455, 315598.4201, ... %! 315600.7143, 315602.9508, 315605.1765 ]; %! [p1, s1] = polyfit (x, y, 10); %! [p2, s2, mu] = polyfit (x, y, 10); %! assert (s2.normr < s1.normr); %!test %! warning ("off", "Octave:nearly-singular-matrix", "local"); %! x = 1000 + (-5:5); %! xn = (x - mean (x)) / std (x); %! pn = ones (1,5); %! y = polyval (pn, xn); %! n = numel (pn) - 1; %! [p, s, mu] = polyfit (x, y, n); %! [p2, s2] = polyfit (x, y, n); %! assert (p, pn, s.normr); %! assert (s.yf, y, s.normr); %! assert (mu, [mean(x), std(x)]); %! assert (s.normr/s2.normr < sqrt (eps)); ## Complex polynomials %!test %! x = 1:4; %! p0 = [1i, 0, 2i, 4]; %! y = polyval (p0, x); %! n = numel (p0) - 1; %! p = polyfit (x, y, n); %! assert (p, p0, 1000*eps); ## Matrix input %!test %! x = [1, 2, 3; 4, 5, 6]; %! y = [0, 0, 1; 1, 0, 0]; %! p = polyfit (x, y, 5); %! expected = [0, 1, -14, 65, -112, 60] / 12; %! assert (p, expected, sqrt (eps)); ## Orientation of output %!test %! x = 0:5; %! y = x.^4 + 2*x + 5; %! [p, s] = polyfit (x, y, 3); %! assert (isrow (s.yf)); %! [p, s] = polyfit (x, y.', 3); %! assert (iscolumn (s.yf)); ## Insufficient data for fit %!test %! x = [1, 2]; %! y = [3, 4]; %! ## Disable warnings entirely because there is not a specific ID to disable. %! wstate = warning (); %! unwind_protect %! warning ("off", "all"); %! p0 = polyfit (x, y, 4); %! [p1, s, mu] = polyfit (x, y, 4); %! unwind_protect_cleanup %! warning (wstate); %! end_unwind_protect %! assert (p0, [0, 0, 0, 1, 2], 10*eps); %! assert (p1, [0, 0, 0, sqrt(2)/2, 3.5], 10*eps); %! assert (size (s.X), [2, 5]); %! assert (s.X(:,1:3), zeros (2,3)); %! assert (size (s.R), [2, 5]); %! assert (s.R(:,1:3), zeros (2,3)); %! assert (size (s.C), [2, 5]); %! assert (s.C(:,1:3), zeros (2,3)); %! assert (s.df, 0); %! assert (mu, [1.5, sqrt(2)/2]); %!test %! x = [1, 2, 3]; %! y = 2*x + 1; %! ## Disable warnings entirely because there is not a specific ID to disable. %! wstate = warning (); %! unwind_protect %! warning ("off", "all"); %! p0 = polyfit (x, y, logical ([1, 1, 1, 0 1])); %! [p1, s, mu] = polyfit (x, y, logical ([1, 1, 1, 0 1])); %! unwind_protect_cleanup %! warning (wstate); %! end_unwind_protect %! assert (p0, [0, -2/11, 12/11, 0, 23/11], 10*eps); %! assert (p1, [0, 2, 0, 0, 5], 10*eps); %! assert (size (s.X), [3, 5]); %! assert (s.X(:,[1,4]), zeros (3,2)); %! assert (size (s.R), [3, 5]); %! assert (s.R(:,[1,4]), zeros (3,2)); %! assert (size (s.C), [3, 5]); %! assert (s.C(:,[1,4]), zeros (3,2)); %! assert (s.df, 0); %! assert (mu, [2, 1]); %!test <*57964> %! ## Disable warnings entirely because there is not a specific ID to disable. %! wstate = warning (); %! unwind_protect %! warning ("off", "all"); %! [p, s] = polyfit ([1,2], [3,4], 2); %! unwind_protect_cleanup %! warning (wstate); %! end_unwind_protect %! assert (size (p), [1, 3]); %! assert (size (s.X), [2, 3]); %! assert (s.X(:,1), [0; 0]); %! assert (size (s.R), [2, 3]); %! assert (s.R(:,1), [0; 0]); %! assert (size (s.C), [2, 3]); %! assert (s.C(:,1), [0; 0]); ## Test input validation %!error <Invalid call> polyfit () %!error <Invalid call> polyfit (1) %!error <Invalid call> polyfit (1,2) %!error <X and Y must have the same number of points> polyfit ([1, 2], 1, 1) %!error <X and Y must have the same number of points> polyfit (1, [1, 2], 1) %!error <N must be a non-negative integer> polyfit (1, 2, [1,2]) %!error <N must be a non-negative integer> polyfit (1, 2, -1) %!error <N must be a non-negative integer> polyfit (1, 2, Inf) %!error <N must be a non-negative integer> polyfit (1, 2, 1.5) %!test <*57964> %! fail ("p = polyfit ([1,2], [3,4], 4)", "warning", "solution is not unique");