Mercurial > hg > octave-lojdl
changeset 14377:c87f927d047d
polyfit.m: add the ability to specify the polynomial template.
author | Pascal Dupuis <Pascal.Dupuis@uclouvain.be> |
---|---|
date | Sun, 05 Feb 2012 15:07:45 +0100 |
parents | 7dd6ac033e69 |
children | dbc99d17f0ad |
files | scripts/polynomial/polyfit.m |
diffstat | 1 files changed, 51 insertions(+), 9 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/polynomial/polyfit.m +++ b/scripts/polynomial/polyfit.m @@ -22,7 +22,9 @@ ## @deftypefnx {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 the least-squares-error of the fit to the points -## @code{[@var{x}, @var{y}]}. +## @code{[@var{x}, @var{y}]}. If @var{n} is a logical vector, it is used +## as a mask to selectivelly force the corresponding polynomial +## coefficients to be used or ignored. ## ## The polynomial coefficients are returned in a row vector. ## @@ -35,6 +37,11 @@ ## @item X ## The Vandermonde matrix used to compute the polynomial coefficients. ## +## @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. ## @@ -46,7 +53,9 @@ ## @end table ## ## The second output may be used by @code{polyval} to calculate the -## statistical error limits of the predicted values. +## 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}. ## ## When the third output, @var{mu}, is present the ## coefficients, @var{p}, are associated with a polynomial in @@ -60,6 +69,8 @@ ## Author: KH <Kurt.Hornik@wu-wien.ac.at> ## Created: 13 December 1994 ## Adapted-By: jwe +## Modified on 20120204 by P. Dupuis; added the ability to specify a +## polynomial mask instead of a polynomial degree. function [p, s, mu] = polyfit (x, y, n) @@ -77,8 +88,16 @@ error ("polyfit: X and Y must be vectors of the same size"); endif - if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == fix (n))) - error ("polyfit: N must be a non-negative integer"); + 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; + 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; endif y_is_row_vector = (rows (y) == 1); @@ -92,10 +111,15 @@ v = vander (x, n+1); ## Solve by QR decomposition. - [q, r, k] = qr (v, 0); + [q, r, k] = qr (v(:, polymask), 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) yf = v*p; @@ -104,10 +128,28 @@ else s.yf = yf; endif + s.X = v; - s.R = r; - s.X = v; - s.df = l - n - 1; + ## 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. + try + C = cholinv (r.'*r)(k, k); + catch + C = NaN * ones (m+1, m+1); + end_try_catch + + if (n ~= m) + ## fill matrices if required + 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; + else + s.R = r; + s.C = C; + endif + s.df = l - m - 1; s.normr = norm (yf - y); endif