[Scilab-users] Scilab leastsq exponential fitting

[Antoine Monmayrant]

On 07/23/2014 01:36 PM, chloe.kykam wrote:
> I have added the dot to the division but I am still getting errors,
>
>   !--error 21
> Invalid index.
> at line     118 of function numderivative called by :
> at line       2 of function Dfun called by :
> at line       2 of function %opt called by :
> at line      92 of function leastsq called by :
> [f,xopt,gopt]=leastsq(list(myfun,x,y,w),a0)
> at line      16 of exec file called by :
> exec('C:\Users\Kying\Desktop\Stray_light\testing.sce', -1)
>
>
>
> --
> View this message in context: http://mailinglists.scilab.org/Scilab-leastsq-exponential-fitting-tp4030949p4030952.html
> Sent from the Scilab users - Mailing Lists Archives mailing list archive at Nabble.com.
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users
>
Hi,

This does work for me:

/////////////////////////////////////
k=6.63e-34*3e8/1.38e-23
x=[1;2;3;4;5;6;7;8;9;10]
y=[280;320;369.22772;391.25743;414.74257;439.75248;466.06931;493.60396;523.87129;530] 

w=[0;0;1;1;1;1;1;1;1;0]

function y=yth(x,a)
     y=a(1)*exp(-k./x/a(2))
endfunction

a0=[1.0;1.0]

function e=myfun(a,x,y,w)
     e=w.*(yth(x,a)-y)
endfunction

[f,xopt,gopt]=leastsq(list(myfun,x,y,w),a0)

scf();
plot(x,yth(x,a0),'g')//initial guess
plot(x,y,'k.')//experimental
plot(x,yth(x,xopt),'r')//actual fit
/////////////////////////////////////

Be careful however:
1) Your yth funcion is not defined for a(2)=0, you should provide 
boundaries to prevent the optimisation routine to try such a value for 
a(2). See help leastsq.
2) Your initial guess seems to be quite crude, you might want to get a 
better guess to start from.
3) The fit does not seems that good (ie y and yth(x,xopt) differ quite a 
lot) and you should take your values xopt(1) and xopt(2) with a TON of salt!

Moreover, your error seems to be related to the calculation of the 
numerical derivative.
You can avoid that by providing the Jacobian of your fit function with 
two benefits: 1), it's going to be way faster (in case you need to redo 
this fit a lot of time or on more data points), 2) you'll have a way to 
measure the "ton of salt", ie the confidence interval (ie error bars if 
you want) of xopt.
Cheers,

Antoine