Problem by star values lsqcurvefit

1 view (last 30 days)
Zuyu An
Zuyu An on 21 Oct 2019
Commented: Star Strider on 22 Oct 2019
xdata = ...
[0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08];
ydata = ...
[-0.7597 -1.5641 4.331 10.226 10.328 10.43 9.2075 7.9845 6.9538 5.9227 4.7857 3.6488 2.1603 0.67176 0.22867 -0.21442 -0.10787];
fun = @(x,xdata)x(1)+x(2).*sqrt(x(3)./(2.*pi.*(x(4).*(xdata+x(5))).^3)).*exp(-(x(3).*((x(4).*(xdata+x(5)))-x(6)).^2)/((xdata+x(5)).*2.*x(4).*x(6).^2));
x0 = [ -10, 20, 1, 12, 0.005, 1 ];
x = lsqcurvefit(fun,x0,xdata,ydata)
B = fminsearch(@(b)norm(ydata - fun(b,xdata)), x0)
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','MaxIter', 10000);
lb = [];
ub = [];
times = linspace(xdata(1),xdata(end));
plot(xdata,ydata,'ko',times,fun(x,times),'b-')
legend('Data','Fitted exponential')
title('Data and Fitted Curve')
result is :
无标题.png
and with result from B = fminsearch(@(b)norm(ydata - fun(b,xdata)), x0) become i only worse result.
but by x0 = [ -10, 20, 1, 12, -0.1, 1 ] it should be already a good result, because
-10+20*sqrt(1/(2*pi*(12*(x+0.005))^3))*exp(-(1*((12*(x+0.005))-#2)^2)/((x+0.005)*24*2^2)) in Latex is already so:
why can i not use x0 = [ -10, 20, 1, 12, 0.005, 1 ] als start values, and how can i find a start values to fitting it?
Thank you very much first!

Accepted Answer

Star Strider
Star Strider on 21 Oct 2019
I used the genetic algorithm (ga function) to fit them, and it took several attempts before I got a good set of parameters.
The code I used:
xdata = [0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08];
ydata = [-1.5641 4.331 10.226 10.328 10.43 9.2075 7.9845 6.9538 5.9227 4.7857 3.6488 2.1603 0.67176 0.22867 -0.21442 -0.10787];
fun = @(x,xdata)x(1)+x(2).*sqrt(x(3)./((2.*pi.*x(4).*(xdata+x(5))).^3)).*exp(-x(3).*((x(4).*(xdata+x(5))-x(6)).^3)./(2.*x(4).*(xdata+x(5)).*x(6).^2));
ftns = @(x) norm(ydata-fun(x,xdata));
PopSz = 500;
Parms = 6;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-6, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[x,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],[],[],[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(x)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, x(k1))
end
x0 = x(:);
x = lsqcurvefit(fun,x0,xdata,ydata)
xv = linspace(min(xdata), max(xdata), 250);
figure
plot(xdata, ydata, 'p')
hold on
plot(xv, fun(x,xv), '-r')
grid
It can take a while to get good parameter estimates. For example these:
x =
-8.347582247433408
0.838228906949566
0.441192639609259
0.120820987842452
0.004904510624840
-0.022856747683530
were the result of running this ga code for 75 times just now, then ‘fine tuning’ them with lsqcurvefit. It produced a very good fit, with the fitness value of 1.97.
There are no shortcuts for estimating parameters of complicated functions. I wish there were.
  4 Comments

Sign in to comment.

More Answers (1)

Alex Sha
Alex Sha on 22 Oct 2019
Hi, Zuyu An, what's your exact data?
xdata = [0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08];
ydata = [-0.7597 -1.5641 4.331 10.226 10.328 10.43 9.2075 7.9845 6.9538 5.9227 4.7857 3.6488 2.1603 0.67176 0.22867 -0.21442 -0.10787];
or
xdata = [0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.07 0.075 0.08];
ydata = [-1.5641 4.331 10.226 10.328 10.43 9.2075 7.9845 6.9538 5.9227 4.7857 3.6488 2.1603 0.67176 0.22867 -0.21442 -0.10787];
note the first one set data (0,-0.7579) is missed in the later.
Also, your function foumual below seems to be overfit, which lead to multiple solutions
fun = @(x,xdata)x(1)+x(2).*sqrt(x(3)./(2.*pi.*(x(4).*(xdata+x(5))).^3)).*exp(-(x(3).*((x(4).*(xdata+x(5)))-x(6)).^2)/((xdata+x(5)).*2.*x(4).*x(6).^2));
if the above function adjust to the follow by deleting "x(2)", there will be unique solution with same goodness of fit.
fun = @(x,xdata)x(1)+sqrt(x(3)./(2.*pi.*(x(4).*(xdata+x(5))).^3)).*exp(-(x(3).*((x(4).*(xdata+x(5)))-x(6)).^2)/((xdata+x(5)).*2.*x(4).*x(6).^2));
  2 Comments
Zuyu An
Zuyu An on 22 Oct 2019
It should be with first data(0, -0.7579)

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!