Mercurial > octave
changeset 28502:39e6308e4251
polyfit.m: Overhaul function (bug #57964).
* polyfit.m: Rewrite documentation and add example code.
Reshape all inputs into column vectors for Matlab compatibility.
Validate number of points in x and y are the same.
Force logical input N to be a row vector (no test previously).
Validate that requested polynomial degree, N, can be achieved given
size of data; issue a warning if it cannot and calculate a solution
based on the lowest polynomial powers (x^0, x^1, ...).
Test "isargout (2)" to avoid unnecessarily calculating complex outputs.
Correctly size outputs s.R, s.C even when there is insufficient data.
Correctly restrict degrees of freedom output to range >= 0.
Add BIST tests for logical input N.
Add BIST test for orientation of output.
Add BIST test for insufficient data.
Add BIST test for bug #57964.
Add BIST tests for input validation.
author | Rik <rik@octave.org> |
---|---|
date | Thu, 25 Jun 2020 15:22:37 -0700 |
parents | e5ced6bd5ac0 |
children | f480103d8333 |
files | scripts/polynomial/polyfit.m |
diffstat | 1 files changed, 205 insertions(+), 71 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/polynomial/polyfit.m Wed Jun 24 07:42:37 2020 -0700 +++ b/scripts/polynomial/polyfit.m Thu Jun 25 15:22:37 2020 -0700 @@ -27,25 +27,33 @@ ## @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}]}. +## 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}(:)]}. ## -## 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. +## @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 a row vector. +## 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 R -## Triangular factor R from the QR@tie{}decomposition. +## +## @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 @@ -56,25 +64,37 @@ ## ## @item normr ## The norm of the residuals. -## -## @item yf -## The values of the polynomial for each value of @var{x}. ## @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 +## 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 (s.C)/s.df)*s.normr}. +## @code{sqrt (diag (@var{s.C})/@var{s.df}) * @var{s.normr}}. ## -## When the third output, @var{mu}, is present the coefficients, @var{p}, are -## associated with a polynomial in +## 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}). ## -## This linear transformation of @var{x} improves the numerical stability of -## the fit. +## 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 @@ -84,93 +104,118 @@ 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) - ## Normalized the x values. + ## Center and scale the x values. mu = [mean(x), std(x)]; x = (x - mu(1)) / mu(2); endif - if (! size_equal (x, y)) - error ("polyfit: X and Y must be vectors of the same size"); - 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; - ## n is the polynomial degree as given the polymask size; m is the - ## effective number of used coefficients. - n = length (polymask) - 1; m = sum (polymask) - 1; + 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 = logical (ones (1, n+1)); m = n; + polymask = true (1, n+1); + m = n; + pad_output = false; endif - y_is_row_vector = (rows (y) == 1); - - ## Reshape x & y into column vectors. - l = numel (x); - x = x(:); - y = y(:); + 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. - v = vander (x, n+1); + X = vander (x, n+1); + v = X(:, polymask); ## Solve by QR decomposition. - [q, r, k] = qr (v(:, polymask), 0); + [q, r, k] = qr (v, 0); p = r \ (q' * y); p(k) = p; - if (n != m) - q = p; p = zeros (n+1, 1); - p(polymask) = q; - endif - - if (nargout > 1) + if (isargout (2)) yf = v*p; - if (y_is_row_vector) s.yf = yf.'; else s.yf = yf; endif - s.X = v; + + s.X = X; - ## r.'*r is positive definite if X(:, polymask) is of full rank. - ## Invert it by cholinv to avoid taking the square root of squared - ## quantities. If cholinv fails, then X(:, polymask) is rank deficient - ## and not invertible. + ## 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+1, m+1); + C = NaN (m, m); end_try_catch - if (n != m) - ## fill matrices if required + if (pad_output) s.X(:, ! polymask) = 0; - s.R = zeros (n+1, n+1); s.R(polymask, polymask) = r; - s.C = zeros (n+1, n+1); s.C(polymask, polymask) = C; + 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 = l - m - 1; + + s.df = max (0, nx - m - 1); s.normr = norm (yf - y); endif - ## Return a row vector. - p = p.'; + 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, x.^2+x+1, 2), [1, 1, 1], sqrt (eps)) -%!assert (polyfit (x, x.^2+x+1, 3), [0, 1, 1, 1], sqrt (eps)) -%!fail ("polyfit (x, x.^2+x+1)") -%!fail ("polyfit (x, x.^2+x+1, [])") + +%!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 @@ -188,25 +233,29 @@ %! assert (s2.normr < s1.normr); %!test -%! x = 1:4; -%! p0 = [1i, 0, 2i, 4]; -%! y0 = polyval (p0, x); -%! p = polyfit (x, y0, numel (p0) - 1); -%! assert (p, p0, 1000*eps); - -%!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); -%! [p, s, mu] = polyfit (x, y, numel (pn) - 1); -%! [p2, s2] = polyfit (x, y, numel (pn) - 1); +%! 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]; @@ -214,4 +263,89 @@ %! expected = [0, 1, -14, 65, -112, 60] / 12; %! assert (p, expected, sqrt (eps)); -%!error <vectors of the same size> polyfit ([1, 2; 3, 4], [1, 2, 3, 4], 2) +## 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 polyfit () +%!error polyfit (1) +%!error polyfit (1,2) +%!error polyfit (1,2,3,4,5) +%!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");