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
[email protected]
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
_______________________________________________
users mailing list
[email protected]
http://lists.scilab.org/mailman/listinfo/users
