uncertainty in polyfit from measurements?

I have some measure point and ive fitted it w/ polyfit. The guesses alright but how can i find out the uncertainty in that coefficients? I used to use Origin for this but it crashes all te time so i decided to switch to matlab.

Answer 1

[John D'Errico]

There are several ways to do this, but I'm not certain what you are looking for. Are you looking for the uncertainty in a prediction of a fitted polynomial? Or are you just asking for the uncertainty in the fitted coefficients? Here is an example.

x = randn(100,1);
y = sin(x) + randn(size(x))/100;

P = polyfitn(x,y,3)
P = 
      ModelTerms: [4x1 double]
    Coefficients: [-0.11492 0.0099578 0.94474 -0.0059965]
    ParameterVar: [2.2444e-06 5.4231e-06 2.7675e-05 1.3305e-05]
    ParameterStd: [0.0014981 0.0023288 0.0052607 0.0036476]
              R2: 0.99788
            RMSE: 0.029182
        VarNames: {'X1'}

Polyfitn returns a standard deviation and variance for each parameter. The ratio of the coefficient to the standard deviation will give you a measure of the significance of each term.

P.Coefficients./P.ParameterStd
ans =
      -76.712        4.276       179.58      -1.6439

My expectation is the linear and cubic terms will have been significant, since a sine wave has a series approximation with only odd order terms. I really want to compare these numbers to a Student t statistic, but I'm feeling lazy now. The error I've added was large enough that the quadratic term is actually potentially significant.

[Orion]

what is polyfitn? this command is not even in product documentation.

[Walter Roberson]

http://www.mathworks.com/matlabcentral/fileexchange/34765-polyfitn

Answer 2

[Walter Roberson]

Asymptotically, the error for a polynomial of degree N fit from N+1 data points, product(x-A[i],i=1..N), will approach +/- sum(eps(A[i]),i=1..N+1) * x^(N-1) -- assuming, that is, perfect representation of x^N.

(If x itself is subject to error in representation then the question of the uncertainty in the coefficients does not arise in the same form.)

That's as x approaches infinity. Higher errors are possible within min(A[i]) to max(A[i]).

Answer 3

[Bruno Luong]

For Gaussian noise, here is a code to show you how to do without stat toolbox. There is a first part that uses brute force method and the second part that shows the method I propose. Both results allow to check that the estimation is correctly done:

x = linspace(0,1);
yext = 1+2*x+3.*x.^2;
sigmay = 0.1; % standard deviation of the noise data

ntest = 1000;
Ptab = zeros(ntest,3);
for n = 1:ntest
    y = yext+sigmay*randn(size(yext));
    P = polyfit(x,y,2);
    Ptab(n,:) = P;
end

% Brute force estimate of uncertainty of P
std(Ptab,1)

% Compute the uncertainty
X = bsxfun(@power,x(:),2:-1:0);
% if not known sigmay can be estimated as std(polyval(P,x)-y) 
stdP = sqrt(diag(inv(X'*X)))*sigmay

Bruno

Answer 4

[Arjen]

Thank you for your answers, but this was not what i meant. Maybe i should have explained it better. I just found the algebraic code:

    D = size(x)*sum(x.^2)-(sum(x)).^2
    S_a =  y_error*sqrt(12/D)
    S_b =  y_error*sqrt(sum(x.^2)/D)

My measure points where 12 points (x,y) This code works, but only when every y point has the same error.

[Bruno Luong]

Yes you should have explained better. The latest attempt is no better. Bye then.

[John D'Errico]

Explain what you want to see, else nobody will bother to try to guess your goal.

Answer 5

[Star Strider]

When this post first appeared, polyparci did not exist. It does now. It may be helpful.