Error estimation for linear fit parameters accounting for known x- and/or y-uncertainty in fitted data

Hello everybody!

It seems like an easy task, but I've been looking around for quite some while now and couldn't find a proper answer yet... so I hope you guys can help me. :) My problem is the following: I got a set of data (x,y(x)) with corresponding uncertainties (dx, dy). What I would like to do is a weighted linear fit to the data and extrapolate to x=0 to obtain y(0) and dy(0). The weighted linear fit is no problem, it can easily be done by:

lin_fit=fit(x, y, 'poly1', 'weights', weights);

where

weights=sum(dy)/dy.^2

are the normalized weights for individual data points according to their uncertainties. Then I calculate 95% confidece intervals from the output of confint and here the problem arises: the error estimations only depend on the scatter of the data points, but completely ignore the given uncertainties. To give an explicit example: create two vectors x and y which give just a straight line. Fit a linear function an you will get zero error in the estimated fit parameters. Now add error bars but no scatter to the data and the estimated errors for the fit parameters are still zero, independently of how large our assigned error bars are for the individual data points. Can somebody please tell me how I can produce error estimations for the fit parameters that also take into account the error bars assigned to the individual data points rather than "just" using them as a weight for the fit?

Many thanks in advance, Axel.

PS: Concerning x-errors I do not even know how to include them into the fit... please help!

Answer 1

[Pavel Nikitin]

Unfortunately, it is still a problem. I have used three built-in functions ('nlinfit', 'fitlm' and 'fit') for different number of points Nx (from 2 to 10). The experimental error was assumed proportional to 'y' (0.1*y). The uncertainty for the parameters was about 10^-15 -- 10^(-8) for Nx>2, and 0 for Nx=2. Here is some code.

Can someone explain such strange behavior of Matlab?

Nx = 2 : 10;
for n = 1 : length(Nx)
    x = linspace(0,10,Nx(n));  y = x + 3;   w = 0.1*y;
    modelFun = @(p,x) p(1) + p(2)*x;
    init_cond = [2.8,1];
    [p,R,J,~,~] = nlinfit(x,y,modelFun,init_cond,'Weights',1./w.^2);
    ci = nlparci(p,R,'Jacobian',J);
    a(n) = (ci(1,2)+ci(1,1))/2;
    b(n) = (ci(2,2)+ci(2,1))/2;
    err_a(n) = abs(ci(1,2)-ci(1,1))/2;
    err_b(n) = abs(ci(2,2)-ci(2,1))/2;
end
plot(Nx,err_a,'-',  Nx,err_b,'--')

another example

Nx = 2 : 10;
for n = 1 : length(Nx)
    x = linspace(0,10,Nx(n));  y = x + 3;   w = 0.1*y;
    mdl = fitlm(x,y,'Weights',1./w.^2);
    coeff = mdl.Coefficients.Estimate;
    a(n) = coeff(1); b(n) = coeff(2);
    err = mdl.Coefficients.SE;
    err_a(n) = err(1); err_b(n) = err(2);                
end
plot(Nx,err_a,'-',  Nx,err_b,'--')

and one more

Nx = 2 : 10;
for n = 2 : length(Nx)
    x = linspace(0,10,Nx(n));  y = x + 3;   w = 0.1*y;
    [res,gof,output] = fit(x',y',fittype('a*x+b'),'StartPoint',[2,0.8],'Weights',1./w.^2);
    coeff = coeffvalues(res);
    a(n) = coeff(1); b(n) = coeff(2);
    err = confint(res);
    err_a(n) = abs(err(1,1) - err(2,1)) / 2;
    err_b(n) = abs(err(1,2) - err(2,2)) / 2;
end
plot(Nx,err_a,'-',  Nx,err_b,'--')