# HG changeset patch # User Ben Abbott # Date 1231758570 -3600 # Node ID 377d908f7e408b253f7d8eacd954f0dfdb5aa486 # Parent 97eab9de6981c2318404df977078b6eba5252ad2 use QR decomposition and normalization for polyfit; normalization for polyval diff --git a/scripts/polynomial/polyfit.m b/scripts/polynomial/polyfit.m --- a/scripts/polynomial/polyfit.m +++ b/scripts/polynomial/polyfit.m @@ -18,39 +18,38 @@ ## . ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{p}, @var{s}] =} polyfit (@var{x}, @var{y}, @var{n}) +## @deftypefn {Function File} {[@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 -## @iftex -## @tex -## $$ -## \sum_{i=1}^N (p(x_i) - y_i)^2 -## $$ -## @end tex -## @end iftex -## @ifinfo -## @code{sumsq (p(x(i)) - y(i))}, -## @end ifinfo -## to best fit the data in the least squares sense. +## @var{n} that minimizes the least-squares-error of the fit. ## ## The polynomial coefficients are returned in a row vector. ## -## If two output arguments are requested, the second is a structure +## The second output is a structured variable, @var{s}, ## containing the following fields: ## -## @table @code +## @table @samp ## @item R -## The Cholesky factor of the Vandermonde matrix used to compute the -## polynomial coefficients. +## Triangular factor R from the QR decomposition. ## @item X -## The Vandermonde matrix used to compute the polynomial coefficients. +## The Vandermonde matrix used to compute the polynomial coefficients. ## @item df -## The degrees of freedom. +## The degrees of freedom. ## @item normr -## The norm of the residuals. +## The norm of the residuals. ## @item yf -## The values of the polynomial for each value of @var{x}. +## 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. +## +## When the third output, @var{mu}, is present the +## coefficients, @var{p}, are associated with a polynomial in +## @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. +## @seealso{polyval, polyconf, residue} ## @end deftypefn ## Author: KH @@ -59,12 +58,17 @@ function [p, s, mu] = polyfit (x, y, n) - - if (nargin != 3) + if (nargin < 3 || nargin > 4) print_usage (); endif - if (! (isvector (x) && isvector (y) && size_equal (x, y))) + if (nargout > 2) + ## Normalized 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 @@ -74,17 +78,21 @@ y_is_row_vector = (rows (y) == 1); - l = length (x); + ## Reshape x & y into column vectors. + l = numel (x); x = reshape (x, l, 1); y = reshape (y, l, 1); - X = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0)); + ## Construct the Vandermonde matrix. + v = (x * ones (1, n+1)) .^ (ones (l, 1) * (n : -1 : 0)); - p = X \ y; + ## Solve by QR decomposition. + [q, r, k] = qr (v, 0); + p = r \ (y' * q)'; + p(k) = p; if (nargout > 1) - - yf = X*p; + yf = v*p; if (y_is_row_vector) s.yf = yf.'; @@ -92,15 +100,55 @@ s.yf = yf; endif - [s.R, dummy] = chol (X'*X); - s.X = X; + s.R = r; + s.X = v; s.df = l - n - 1; s.normr = norm (yf - y); - endif - ## Return value should be a row vector. - + ## Return a row vector. p = p.'; +## 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 +%! 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 (s1.normr, 0.11264, 0.1) +%! 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 +%! 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); +%! assert (p, pn, s.normr) +%! assert (s.yf, y, s.normr) +%! assert (mu, [mean(x), std(x)]) +%! assert (s.normr/s2.normr < 1e-9) + +%!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)) + + endfunction diff --git a/scripts/polynomial/polyval.m b/scripts/polynomial/polyval.m --- a/scripts/polynomial/polyval.m +++ b/scripts/polynomial/polyval.m @@ -18,15 +18,20 @@ ## . ## -*- texinfo -*- -## @deftypefn {Function File} {} polyval (@var{c}, @var{x}) -## Evaluate a polynomial. -## -## @code{polyval (@var{c}, @var{x})} will evaluate the polynomial at the -## specified value of @var{x}. -## -## If @var{x} is a vector or matrix, the polynomial is evaluated at each of +## @deftypefn {Function File} {@var{y}=} polyval (@var{p}, @var{x}) +## @deftypefnx {Function File} {@var{y}=} polyval (@var{p}, @var{x}, [], @var{mu}) +## Evaluate the polynomial at of the specified values for @var{x}. When @var{mu} +## is present evaluate the polynomial for (@var{x}-@var{mu}(1))/@var{mu}(2). +## If @var{x} is a vector or matrix, the polynomial is evaluated for each of ## the elements of @var{x}. -## @seealso{polyvalm, poly, roots, conv, deconv, residue, filter, +## @deftypefnx {Function File} {[@var{y}, @var{dy}] =} polyval (@var{p}, @var{x}, @var{S}) +## @deftypefnx {Function File} {[@var{y}, @var{dy}] =} polyval (@var{p}, @var{x}, @var{S}, @var{mu}) +## In addition to evaluating the polynomial, the second output +## represents the prediction interval, @var{y} +/- @var{dy}, which +## contains at least 50% of the future predictions. To calculate the +## prediction interval, the structured variable @var{s}, originating +## form `polyfit', must be present. +## @seealso{polyfit, polyvalm, poly, roots, conv, deconv, residue, filter, ## polyderiv, polyinteg} ## @end deftypefn @@ -34,30 +39,92 @@ ## Created: June 1994 ## Adapted-By: jwe -function y = polyval (c, x) +function [y, dy] = polyval (p, x, s, mu) - if (nargin != 2) + if (nargin < 2 || nargin > 4 || (nargout == 2 && nargin < 3)) print_usage (); endif - if (! (isvector (c) || isempty (c))) + if (nargin < 3) + s = []; + endif + + if (! (isvector (p) || isempty (p))) error ("polyval: first argument must be a vector"); endif + if (nargin < 4) + mu = [0, 1]; + endif + if (isempty (x)) y = []; return; endif - if (length (c) == 0) - y = c; + if (length (p) == 0) + y = p; return; endif - n = length (c); - y = c (1) * ones (rows (x), columns (x)); - for index = 2:n - y = c (index) + x .* y; - endfor + n = length (p) - 1; + k = numel (x); + x = (x - mu(1)) / mu(2); + A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0)); + y(:) = A * p(:); + y = reshape (y, size (x)); + + if (nargout == 2) + ## The line below is *not* the result of a conceptual grasp of statistics. + ## Instead, after reading the links below and comparing to the output of Matlab's polyval.m, + ## http://www.mathworks.com/access/helpdesk/help/toolbox/stats/index.html?/access/helpdesk/help/toolbox/stats/finv.html + ## http://www.mathworks.com/access/helpdesk/help/toolbox/curvefit/index.html?/access/helpdesk/help/toolbox/curvefit/bq_5ka6-1_1.html + ## Note: the F-Distribution is generally considered to be single-sided. + ## http://www.itl.nist.gov/div898/handbook/eda/section3/eda3673.htm + ## t = finv (1-alpha, s.df, s.df); + ## dy = t * sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df) + ## If my inference is correct, then t must equal 1 for polyval. + ## This is because finv (0.5, n, n) = 1.0 for any n. + dy = sqrt (1 + sumsq (A/s.R, 2)) * s.normr / sqrt (s.df); + dy = reshape (dy, size (x)); + endif endfunction + +%!test +%! fail("polyval([1,0;0,1],0:10)"); + +%!test +%! r = 0:10:50; +%! p = poly (r); +%! p = p / max(abs(p)); +%! x = linspace(0,50,11); +%! y = polyval(p,x) + 0.25*sin(100*x); +%! [pf, s] = polyfit (x, y, numel(r)); +%! [y1, delta] = polyval (pf, x, s); +%! expected = [0.37235, 0.35854, 0.32231, 0.32448, 0.31328, ... +%! 0.32036, 0.31328, 0.32448, 0.32231, 0.35854, 0.37235]; +%! assert (delta, expected, 0.00001) + +%!test +%! x = 10 + (-2:2); +%! y = [0, 0, 1, 0, 2]; +%! p = polyfit (x, y, numel (x) - 1); +%! [pn, s, mu] = polyfit (x, y, numel (x) - 1); +%! y1 = polyval (p, x); +%! yn = polyval (pn, x, [], mu); +%! assert (y1, y, sqrt(eps)) +%! assert (yn, y, sqrt(eps)) + +%!test +%! p = [0, 1, 0]; +%! x = 1:10; +%! assert (x, polyval(p,x), eps) +%! x = x(:); +%! assert (x, polyval(p,x), eps) +%! x = reshape(x, [2, 5]); +%! assert (x, polyval(p,x), eps) +%! x = reshape(x, [5, 2]); +%! assert (x, polyval(p,x), eps) +%! x = reshape(x, [1, 1, 5, 2]); +%! assert (x, polyval(p,x), eps)