octave_packages/optim-1.2.0/wpolyfit.m

   1 ## Author: Paul Kienzle <pkienzle@gmail.com>
   2 ## This program is granted to the public domain.
   3
   4 ## -*- texinfo -*-
   5 ## @deftypefn {Function File} {[@var{p}, @var{s}] =} wpolyfit (@var{x}, @var{y}, @var{dy}, @var{n})
   6 ## Return the coefficients of a polynomial @var{p}(@var{x}) of degree
   7 ## @var{n} that minimizes
   8 ## @iftex
   9 ## @tex
  10 ## $$
  11 ## \sum_{i=1}^N (p(x_i) - y_i)^2
  12 ## $$
  13 ## @end tex
  14 ## @end iftex
  15 ## @ifinfo
  16 ## @code{sumsq (p(x(i)) - y(i))},
  17 ## @end ifinfo
  18 ## to best fit the data in the least squares sense.  The standard error
  19 ## on the observations @var{y} if present are given in @var{dy}.
  20 ##
  21 ## The returned value @var{p} contains the polynomial coefficients
  22 ## suitable for use in the function polyval.  The structure @var{s} returns
  23 ## information necessary to compute uncertainty in the model.
  24 ##
  25 ## To compute the predicted values of y with uncertainty use
  26 ## @example
  27 ## [y,dy] = polyconf(p,x,s,'ci');
  28 ## @end example
  29 ## You can see the effects of different confidence intervals and
  30 ## prediction intervals by calling the wpolyfit internal plot
  31 ## function with your fit:
  32 ## @example
  33 ## feval('wpolyfit:plt',x,y,dy,p,s,0.05,'pi')
  34 ## @end example
  35 ## Use @var{dy}=[] if uncertainty is unknown.
  36 ##
  37 ## You can use a chi^2 test to reject the polynomial fit:
  38 ## @example
  39 ## p = 1-chi2cdf(s.normr^2,s.df);
  40 ## @end example
  41 ## p is the probability of seeing a chi^2 value higher than that which
  42 ## was observed assuming the data are normally distributed around the fit.
  43 ## If p < 0.01, you can reject the fit at the 1% level.
  44 ##
  45 ## You can use an F test to determine if a higher order polynomial
  46 ## improves the fit:
  47 ## @example
  48 ## [poly1,S1] = wpolyfit(x,y,dy,n);
  49 ## [poly2,S2] = wpolyfit(x,y,dy,n+1);
  50 ## F = (S1.normr^2 - S2.normr^2)/(S1.df-S2.df)/(S2.normr^2/S2.df);
  51 ## p = 1-f_cdf(F,S1.df-S2.df,S2.df);
  52 ## @end example
  53 ## p is the probability of observing the improvement in chi^2 obtained
  54 ## by adding the extra parameter to the fit.  If p < 0.01, you can reject
  55 ## the lower order polynomial at the 1% level.
  56 ##
  57 ## You can estimate the uncertainty in the polynomial coefficients
  58 ## themselves using
  59 ## @example
  60 ## dp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
  61 ## @end example
  62 ## but the high degree of covariance amongst them makes this a questionable
  63 ## operation.
  64 ##
  65 ## @deftypefnx {Function File} {[@var{p}, @var{s}, @var{mu}] =} wpolyfit (...)
  66 ##
  67 ## If an additional output @code{mu = [mean(x),std(x)]} is requested then
  68 ## the @var{x} values are centered and normalized prior to computing the fit.
  69 ## This will give more stable numerical results.  To compute a predicted
  70 ## @var{y} from the returned model use
  71 ## @code{y = polyval(p, (x-mu(1))/mu(2)}
  72 ##
  73 ## @deftypefnx {Function File} wpolyfit (...)
  74 ##
  75 ## If no output arguments are requested, then wpolyfit plots the data,
  76 ## the fitted line and polynomials defining the standard error range.
  77 ##
  78 ## Example
  79 ## @example
  80 ## x = linspace(0,4,20);
  81 ## dy = (1+rand(size(x)))/2;
  82 ## y = polyval([2,3,1],x) + dy.*randn(size(x));
  83 ## wpolyfit(x,y,dy,2);
  84 ## @end example
  85 ##
  86 ## @deftypefnx {Function File} wpolyfit (..., 'origin')
  87 ##
  88 ## If 'origin' is specified, then the fitted polynomial will go through
  89 ## the origin.  This is generally ill-advised.  Use with caution.
  90 ##
  91 ## Hocking, RR (2003). Methods and Applications of Linear Models.
  92 ## New Jersey: John Wiley and Sons, Inc.
  93 ##
  94 ## @end deftypefn
  95 ## @seealso{polyfit,polyconf}
  96
  97 function [p_out, s, mu] = wpolyfit (varargin)
  98
  99   ## strip 'origin' of the end
 100   args = length(varargin);
 101   if args>0 && ischar(varargin{args})
 102     origin = varargin{args};
 103     args--;
 104   else
 105     origin='';
 106   endif
 107   ## strip polynomial order off the end
 108   if args>0
 109     n = varargin{args};
 110     args--;
 111   else
 112     n = [];
 113   end
 114   ## interpret the remainder as x,y or x,y,dy or [x,y] or [x,y,dy]
 115   if args == 3
 116     x = varargin{1};
 117     y = varargin{2};
 118     dy = varargin{3};
 119   elseif args == 2
 120     x = varargin{1};
 121     y = varargin{2};
 122     dy = [];
 123   elseif args == 1
 124     A = varargin{1};
 125     [nr,nc]=size(A);
 126     if all(nc!=[2,3])
 127       error("wpolyfit expects vectors x,y,dy or matrix [x,y,dy]");
 128     endif
 129     dy = [];
 130     if nc == 3, dy = A(:,3); endif
 131     y = A(:,2);
 132     x = A(:,1);
 133   else
 134     usage ("wpolyfit (x, y [, dy], n [, 'origin'])");
 135   end
 136
 137   if (length(origin) == 0)
 138     through_origin = 0;
 139   elseif strcmp(origin,'origin')
 140     through_origin = 1;
 141   else
 142     error ("wpolyfit: expected 'origin' but found '%s'", origin)
 143   endif
 144
 145   if any(size (x) != size (y))
 146     error ("wpolyfit: x and y must be vectors of the same size");
 147   endif
 148   if length(dy)>1 && length(y) != length(dy)
 149     error ("wpolyfit: dy must be a vector the same length as y");
 150   endif
 151
 152   if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == round (n)))
 153     error ("wpolyfit: n must be a nonnegative integer");
 154   endif
 155
 156   if nargout == 3
 157     mu = [mean(x), std(x)];
 158     x = (x - mu(1))/mu(2);
 159   endif
 160
 161   k = length (x);
 162
 163   ## observation matrix
 164   if through_origin
 165     ## polynomial through the origin y = ax + bx^2 + cx^3 + ...
 166     A = (x(:) * ones(1,n)) .^ (ones(k,1) * (n:-1:1));
 167   else
 168     ## polynomial least squares y = a + bx + cx^2 + dx^3 + ...
 169     A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
 170   endif
 171
 172   [p,s] = wsolve(A,y(:),dy(:));
 173
 174   if through_origin
 175     p(n+1) = 0;
 176   endif
 177
 178   if nargout == 0
 179     good_fit = 1-chi2cdf(s.normr^2,s.df);
 180     printf("Polynomial: %s  [ p(chi^2>observed)=%.2f%% ]\n", polyout(p,'x'), good_fit*100);
 181     plt(x,y,dy,p,s,'ci');
 182   else
 183     p_out = p';
 184   endif
 185
 186 function plt(x,y,dy,p,s,varargin)
 187
 188   if iscomplex(p)
 189     # XXX FIXME XXX how to plot complex valued functions?
 190     # Maybe using hue for phase and saturation for magnitude
 191     # e.g., Frank Farris (Santa Cruz University) has this:
 192     # http://www.maa.org/pubs/amm_complements/complex.html
 193     # Could also look at the book
 194     #   Visual Complex Analysis by Tristan Needham, Oxford Univ. Press
 195     # but for now we punt
 196     return
 197   end
 198
 199   ## decorate the graph
 200   grid('on');
 201   xlabel('abscissa X'); ylabel('data Y');
 202   title('Least-squares Polynomial Fit with Error Bounds');
 203
 204   ## draw fit with estimated error bounds
 205   xf = linspace(min(x),max(x),150)';
 206   [yf,dyf] = polyconf(p,xf,s,varargin{:});
 207   plot(xf,yf+dyf,"g.;;", xf,yf-dyf,"g.;;", xf,yf,"g-;fit;");
 208
 209   ## plot the data
 210   hold on;
 211   if (isempty(dy))
 212     plot(x,y,"x;data;");
 213   else
 214     if isscalar(dy), dy = ones(size(y))*dy; end
 215     errorbar (x, y, dy, "~;data;");
 216   endif
 217   hold off;
 218
 219   if strcmp(deblank(input('See residuals? [y,n] ','s')),'y')
 220     clf;
 221     if (isempty(dy))
 222       plot(x,y-polyval(p,x),"x;data;");
 223     else
 224       errorbar(x,y-polyval(p,x),dy, '~;data;');
 225     endif
 226     hold on;
 227     grid on;
 228     ylabel('Residuals');
 229     xlabel('abscissa X');
 230     plot(xf,dyf,'g.;;',xf,-dyf,'g.;;');
 231     hold off;
 232   endif
 233
 234 %!demo % #1
 235 %! x = linspace(0,4,20);
 236 %! dy = (1+rand(size(x)))/2;
 237 %! y = polyval([2,3,1],x) + dy.*randn(size(x));
 238 %! wpolyfit(x,y,dy,2);
 239
 240 %!demo % #2
 241 %! x = linspace(-i,+2i,20);
 242 %! noise = ( randn(size(x)) + i*randn(size(x)) )/10;
 243 %! P = [2-i,3,1+i];
 244 %! y = polyval(P,x) + noise;
 245 %! wpolyfit(x,y,2)
 246
 247 %!demo
 248 %! pin = [3; -1; 2];
 249 %! x = -3:0.1:3;
 250 %! y = polyval (pin, x);
 251 %!
 252 %! ## Poisson weights
 253 %! # dy = sqrt (abs (y));
 254 %! ## Uniform weights in [0.5,1]
 255 %! dy = 0.5 + 0.5 * rand (size (y));
 256 %!
 257 %! y = y + randn (size (y)) .* dy;
 258 %! printf ("Original polynomial: %s\n", polyout (pin, 'x'));
 259 %! wpolyfit (x, y, dy, length (pin)-1);
 260
 261