wpolyfit (optim

2.11 Polynomial fitting suitable for polyconf

Function File: [p, s] = wpolyfit (x, y, dy, n)

Return the coefficients of a polynomial p(x) of degree n that minimizes sumsq (p(x(i)) - y(i)), to best fit the data in the least squares sense. The standard error on the observations y if present are given in dy.

The returned value p contains the polynomial coefficients suitable for use in the function polyval. The structure s returns information necessary to compute uncertainty in the model.

To compute the predicted values of y with uncertainty use

[y,dy] = polyconf(p,x,s,'ci');

You can see the effects of different confidence intervals and prediction intervals by calling the wpolyfit internal plot function with your fit:

feval('wpolyfit:plt',x,y,dy,p,s,0.05,'pi')

Use dy=[] if uncertainty is unknown.

You can use a chi^2 test to reject the polynomial fit:

p = 1-chi2cdf(s.normr^2,s.df);

p is the probability of seeing a chi^2 value higher than that which was observed assuming the data are normally distributed around the fit. If p < 0.01, you can reject the fit at the 1% level.

You can use an F test to determine if a higher order polynomial improves the fit:

[poly1,S1] = wpolyfit(x,y,dy,n);
[poly2,S2] = wpolyfit(x,y,dy,n+1);
F = (S1.normr^2 - S2.normr^2)/(S1.df-S2.df)/(S2.normr^2/S2.df);
p = 1-f_cdf(F,S1.df-S2.df,S2.df);

p is the probability of observing the improvement in chi^2 obtained by adding the extra parameter to the fit. If p < 0.01, you can reject the lower order polynomial at the 1% level.

You can estimate the uncertainty in the polynomial coefficients themselves using

dp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;

but the high degree of covariance amongst them makes this a questionable operation.

Function File: [p, s, mu] = wpolyfit (...): If an additional output mu = [mean(x),std(x)] is requested then the x values are centered and normalized prior to computing the fit. This will give more stable numerical results. To compute a predicted y from the returned model use y = polyval(p, (x-mu(1))/mu(2)

Function File: wpolyfit (...)

If no output arguments are requested, then wpolyfit plots the data, the fitted line and polynomials defining the standard error range.

Example

x = linspace(0,4,20);
dy = (1+rand(size(x)))/2;
y = polyval([2,3,1],x) + dy.*randn(size(x));
wpolyfit(x,y,dy,2);

Function File: wpolyfit (..., 'origin')

If ’origin’ is specified, then the fitted polynomial will go through the origin. This is generally ill-advised. Use with caution.

Hocking, RR (2003). Methods and Applications of Linear Models. New Jersey: John Wiley and Sons, Inc.

See also: polyconf.