[Scilab-users] Fitting correlations to measured data. Write-ups of leastsq is messy ...

[Dang Ngoc Chan, Christophe]

Hello,

> De : Heinz Nabielek
> Envoyé : jeudi 30 janvier 2020 23:25
>
> First: linear least-squares fitting: with measurement data x=[2 7 12]'; y=[2 4.5
> 6.5]'; it is easy, to fit a straight line and plot it:
> M=[ones(x) x]; a=M\y; plot2d(x,M*a); plot(x,y,'r.');
> [...]
> In my case, I wanted to simultaneously fit three straight lines through three
> measurement series with the condition that all 3 lines start at the same point
> on the x-axis. So the following worked well for me:
>
> x=[2 7 12]'; y=[2 4.5  6.5]'; y1=[1 2 3]'; y2=[0 1 1.4]'; par0 = [.5 .15 .09 -5];
> function e=err(par, x,y,y1,y2)
>    e= [y-par(1)*(x+par(4)); y1-par(2)*(x+par(4)); y2-par(3)*(x+par(4))]
> endfunction
>
>  [f,paropt] = leastsq ( list(err,x,y,y1,y2), par0);

For this specific case, you might consider this:

You can use a linear law that goes through zero with a linear regression just by removing the column of ones:

A = x\y

You can thus improve the speed by a factor 8 :


----------

x=[2 7 12]';
y=[2 4.5  6.5]';
y1=[1 2 3]';
y2=[0 1 1.4]';
par0 = [.5 .15 .09 -5];
par01 = -5;
function e=err(par, x,y,y1,y2)
    e= [y-par(1)*(x+par(4)); y1-par(2)*(x+par(4)); y2-par(3)*(x+par(4))]
 endfunction

function e=err1(par, x,y,y1,y2)
    A = (x + par)\y;
    A1 = (x + par)\y1;
    A2 = (x + par)\y2;
    e= [y-A*(x+par); y1-A1*(x+par); y2-A2*(x+par)]
 endfunction

tic();
[f,paropt] = leastsq ( list(err,x,y,y1,y2), par0);
t1=toc();

figure(0);
clf();

plot(x,[y y1 y2],'.');
xx=[5, 15];
plot(xx,paropt(1)*(xx+paropt(4)),'b--');
plot(xx,paropt(2)*(xx+paropt(4)),'g--');
plot(xx,paropt(3)*(xx+paropt(4)),'r--');

tic();
[f,paropt1] = leastsq ( list(err1,x,y,y1,y2), par01);
A = (x+paropt1)\y;
A1 = (x+paropt1)\y1;
A2 = (x+paropt1)\y2;
t2=toc();

figure(1);
clf();

plot(x,[y y1 y2],'.');
xx=[5, 15];
plot(xx,A.*(xx+paropt1),'b--');
plot(xx,A1.*(xx+paropt1),'g--');
plot(xx,A2.*(xx+paropt1),'r--');

disp("t1 = " + string(t1) + " ; t2 = " + string(t2) +...
" ; t1/t2 = " + string(t1/t2));

----------

--
Christophe Dang Ngoc Chan
Mechanical calculation engineer

General
This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error), please notify the sender immediately and destroy this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.