[Scilab-users] the linear least-squares fitting problem

Antoine Monmayrant antoine.monmayrant at laas.fr
Wed Aug 19 12:54:18 CEST 2020


Hello Heinz,

Please see the attached tutorial¹ I wrote for the other physicists in my 
lab on how to use scilab to fit data with a model and calculate the 
confidence interval.
The reference paper mentioned in the comments is worth reading, it 
helped me a lot understanding how the confidence interval relates to the 
covariance matrix.

Hope it helps,

Antoine
¹I made a poor choice for the data/model to fit as sin(x) is not really 
nice (multiple parameters work equally well). I should have used a 
gaussian or lorentzian curve...

Le 01/08/2020 à 16:45, Heinz Nabielek a écrit :
> Dear Scilabers,
>
> for mathematicians, the linear least-squares fitting is a trivial 
> problem, but as an experimental physicist I had problems finding the 
> right correlations for predicting the Confidence Limits of the fitted 
> functional relationship. For a trivial set of experimental data, I 
> wrote a few lines of Scilab code and have 2 questions:
>
> - is the computation of confidence limits correct?
>
> - for large data sets and higher degree polynomials, can we make the 
> computations more efficient and reliable?
>
> Thanks in advance
> Heinz
> _________________________
> x=[0.1269868 0.135477 0.221034 0.45 0.8350086 0.9057919 0.9133759 
> 0.9688678]'; y=[0.0695843 0.0154644 0.0893982 0.5619441 1.3879476 
> 1.5639274 1.6570076 1.9061488]'; 
> plot(x,y,'m.','markersize',14);xgrid(1,1,1); X=[ones(x) x x^2]; // 
> matrix of independent vectors b=X\y; // parameters of linear LSFIT 
> with a 2nd degree polynomial yh = X*b; //Fitted values yh to 
> approximate measured y's e=y-yh; //Errors or residuals SSE=e'*e; //Sum 
> of squared errors ybar=mean(y); R2=1-SSE/sum((y-ybar)^2); [m 
> n]=size(X); MSE = SSE/(m-n-1); //Mean square error se=sqrt(MSE); 
> C=MSE*inv(X'*X); seb=sqrt(diag(C)); CONF=.95; alpha=1-CONF; ta2 = 
> cdft('T',m-n,1-alpha/2,alpha/2); //t-value for alpha/2 
> xx=(0.1:.025:1)'; // cover the whole abscissa range in detail 
> XX=[ones(xx) xx xx^2]; yhh=XX*b; // predicted extended LSFIT curve 
> plot(xx,yhh,'kd-','MarkerSize',8); [mm n]=size(XX); sY = []; sYp = []; 
> //Terms involved in Confidence Interval for Y, Ypred for i=1:mm; sY = 
> [sY; sqrt(XX(i,:)*C*XX(i,:)')]; // standard error in fit sYp = [sYp; 
> se*sqrt(1+XX(i,:)*(C/se)*XX(i,:)')]; // standard error in spread in 
> population distribution end; plot(xx, yhh+ta2*sYp,'rd-'); plot(xx, 
> yhh+ta2*sY,'bd-'); plot(xx, yhh-ta2*sY,'gd-'); plot(xx, 
> yhh-ta2*sYp,'md-'); title('LINEAR LEAST-SQUARES FIT WITH 95% 
> CONFIDENCE BOUNDS ON FIT AND ON THE POPULATION SPREAD','fontsize',2); 
> xlabel('x-values','fontsize',5); ylabel('y-values','fontsize',5); 
> legend('measurement data','LSFIT of 2nd order polynomial','upper 95% 
> CL population','upper 95% CL fit','lower 95% CL fit','lower 95% CL 
> population',2);
>
>
>
>
>
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20200819/7a1ebe31/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: tuto_fit.sce
Type: application/x-scilab-sce
Size: 11713 bytes
Desc: not available
URL: <https://lists.scilab.org/pipermail/users/attachments/20200819/7a1ebe31/attachment.bin>


More information about the users mailing list