Computing variance of a variable after linear fit using lsqcurvefit function

I just want to confirm whether the follwings are correct or not. I have two data sets and I want to linealy fit these data sets using the "lsqcurvefit" function. The linear fit equation is given by :. Following is the code, which also computes the covariance matrix of the parametrs (a,b):

clear all;
xdata = [-0.41095  -0.410820  -0.41074  -0.41068  -0.41063  -0.41058  -0.41055  -0.41052  -0.41049  -0.41047];
ydata = [0.166666666666667   0.142857142857143   0.125000000000000   0.111111111111111   0.100000000000000   0.090909090909091 0.083333333333333   0.076923076923077   0.071428571428571   0.066666666666667];
fun = @(A,xdata)(A(1)+A(2)*xdata);
A0 = [1,-1];
A = lsqcurvefit(fun,A0,xdata,ydata);
[A,resnorm,residual,exitflag,output,lambda,J]= lsqcurvefit(fun,A0,xdata,ydata);
ACovariance = inv(J.'*J)*var(residual) %computes the covariance matrix of the parameters a,b in the equation y=a+bx

Is the above code correct for computing the covariance matrix of the fitting parameters?

Now I want to compute the variance of "y" at some point let's say at . Is the following equation correct for computing the variance of : . Where, are the variance of a and b , which are equal to i think the diagonal elements of the matrix "ACovariance" in the code, and is the covariance of the parametrs a,b , which is equal to i think offdiagonal elements of "ACovariance". Following is a derivation of the above equation:

Now variance

where, variance .

Is the above way of computing variance of y at a point , is correct? If the above are incorrect, could you please, point out how to find variance of the variable "y" extrapolated to some arbitary point, after linear fitting the given data sets.

[Matt J]

I want to point out that you would never use lsqcurvefit to fit a linear unconstrained model like this. You would just use polyfit.

Answer 1

[Matt J]

It would be,

yvariance = [1,x0]*ACovariance*[1;x0]

[P Rakesh Kumar Dora]

I think your expression of variance matches with my expression in question.

[Matt J]

Does that mean your question has been answered? If so, please Accept-click the answer.

[P Rakesh Kumar Dora]

I am not getting the correct answer.

[Matt J]

Well, provide the Jacobian you've computed, and the x0 you are testing in .mat file attachment.

Also say what you think is the correct answer and why.

[P Rakesh Kumar Dora]

clear all;
ydata = [-0.41095  -0.410820  -0.41074  -0.41068  -0.41063  -0.41058  -0.41055  -0.41052  -0.41049  -0.41047];
xdata = [0.166666666666667   0.142857142857143   0.125000000000000   0.111111111111111   0.100000000000000   0.090909090909091 0.083333333333333   0.076923076923077   0.071428571428571   0.066666666666667];
fun = @(A,xdata)(A(1)+A(2)*xdata);
A0 = [1,-1];
A = lsqcurvefit(fun,A0,xdata,ydata);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

[A,resnorm,residual,exitflag,output,lambda,J]= lsqcurvefit(fun,A0,xdata,ydata);

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

ACovariance = inv(J.'*J)*var(residual);
J

J = 
   (1,1)       1.0000
   (2,1)       1.0000
   (3,1)       1.0000
   (4,1)       1.0000
   (5,1)       1.0000
   (6,1)       1.0000
   (7,1)       1.0000
   (8,1)       1.0000
   (9,1)       1.0000
  (10,1)       1.0000
   (1,2)       0.1667
   (2,2)       0.1429
   (3,2)       0.1250
   (4,2)       0.1111
   (5,2)       0.1000
   (6,2)       0.0909
   (7,2)       0.0833
   (8,2)       0.0769
   (9,2)       0.0714
  (10,2)       0.0667

yvariance = [1 0]*ACovariance*[1;0]

yvariance = 2.8711e-11

ystandarddeviation= sqrt(yvariance)

ystandarddeviation = 5.3583e-06

[P Rakesh Kumar Dora]

Here the point x0 = 0. I expect the standard deviation of y to be around 0.001. I have this number from a reference, I don't know how they compute. I also updated my datasets, in the main question.

[Matt J]

Well, the std deviation of your residuals is orders of magnitude smaller than that,

>> std(residual)

ans =

   4.8895e-06

I don't see how extrapolating to x0=0 ends up magnifying that 1000 times.