1 ## Author: Paul Kienzle <pkienzle@gmail.com>
2 ## This program is granted to the public domain.
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
11 ## \sum_{i=1}^N (p(x_i) - y_i)^2
16 ## @code{sumsq (p(x(i)) - y(i))},
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}.
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.
25 ## To compute the predicted values of y with uncertainty use
27 ## [y,dy] = polyconf(p,x,s,'ci');
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:
33 ## feval('wpolyfit:plt',x,y,dy,p,s,0.05,'pi')
35 ## Use @var{dy}=[] if uncertainty is unknown.
37 ## You can use a chi^2 test to reject the polynomial fit:
39 ## p = 1-chi2cdf(s.normr^2,s.df);
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.
45 ## You can use an F test to determine if a higher order polynomial
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);
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.
57 ## You can estimate the uncertainty in the polynomial coefficients
60 ## dp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
62 ## but the high degree of covariance amongst them makes this a questionable
65 ## @deftypefnx {Function File} {[@var{p}, @var{s}, @var{mu}] =} wpolyfit (...)
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)}
73 ## @deftypefnx {Function File} wpolyfit (...)
75 ## If no output arguments are requested, then wpolyfit plots the data,
76 ## the fitted line and polynomials defining the standard error range.
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);
86 ## @deftypefnx {Function File} wpolyfit (..., 'origin')
88 ## If 'origin' is specified, then the fitted polynomial will go through
89 ## the origin. This is generally ill-advised. Use with caution.
91 ## Hocking, RR (2003). Methods and Applications of Linear Models.
92 ## New Jersey: John Wiley and Sons, Inc.
95 ## @seealso{polyfit,polyconf}
97 function [p_out, s, mu] = wpolyfit (varargin)
99 ## strip 'origin' of the end
100 args = length(varargin);
101 if args>0 && ischar(varargin{args})
102 origin = varargin{args};
107 ## strip polynomial order off the end
114 ## interpret the remainder as x,y or x,y,dy or [x,y] or [x,y,dy]
127 error("wpolyfit expects vectors x,y,dy or matrix [x,y,dy]");
130 if nc == 3, dy = A(:,3); endif
134 usage ("wpolyfit (x, y [, dy], n [, 'origin'])");
137 if (length(origin) == 0)
139 elseif strcmp(origin,'origin')
142 error ("wpolyfit: expected 'origin' but found '%s'", origin)
145 if any(size (x) != size (y))
146 error ("wpolyfit: x and y must be vectors of the same size");
148 if length(dy)>1 && length(y) != length(dy)
149 error ("wpolyfit: dy must be a vector the same length as y");
152 if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == round (n)))
153 error ("wpolyfit: n must be a nonnegative integer");
157 mu = [mean(x), std(x)];
158 x = (x - mu(1))/mu(2);
163 ## observation matrix
165 ## polynomial through the origin y = ax + bx^2 + cx^3 + ...
166 A = (x(:) * ones(1,n)) .^ (ones(k,1) * (n:-1:1));
168 ## polynomial least squares y = a + bx + cx^2 + dx^3 + ...
169 A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
172 [p,s] = wsolve(A,y(:),dy(:));
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');
186 function plt(x,y,dy,p,s,varargin)
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
199 ## decorate the graph
201 xlabel('abscissa X'); ylabel('data Y');
202 title('Least-squares Polynomial Fit with Error Bounds');
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;");
214 if isscalar(dy), dy = ones(size(y))*dy; end
215 errorbar (x, y, dy, "~;data;");
219 if strcmp(deblank(input('See residuals? [y,n] ','s')),'y')
222 plot(x,y-polyval(p,x),"x;data;");
224 errorbar(x,y-polyval(p,x),dy, '~;data;');
229 xlabel('abscissa X');
230 plot(xf,dyf,'g.;;',xf,-dyf,'g.;;');
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);
241 %! x = linspace(-i,+2i,20);
242 %! noise = ( randn(size(x)) + i*randn(size(x)) )/10;
244 %! y = polyval(P,x) + noise;
250 %! y = polyval (pin, x);
252 %! ## Poisson weights
253 %! # dy = sqrt (abs (y));
254 %! ## Uniform weights in [0.5,1]
255 %! dy = 0.5 + 0.5 * rand (size (y));
257 %! y = y + randn (size (y)) .* dy;
258 %! printf ("Original polynomial: %s\n", polyout (pin, 'x'));
259 %! wpolyfit (x, y, dy, length (pin)-1);