Problem by star values lsqcurvefit

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 :

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!

Answer 1

[Star Strider]

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.

[Zuyu An]

but why can i not use x0 = [ -10, 20, 1, 12, 0.005, 1 ] als start values, it is already a good one.

und unfortunately your code in my Matlab funktioniert nicht

and when i use

x =

-8.347582247433408

0.838228906949566

0.441192639609259

0.120820987842452

0.004904510624840

-0.022856747683530

als start valuse it is also not fitting.....

i am sorry i have so many questions

[Star Strider]

You need to have the Global Optimization Toolbox to use the ga function, although you can write a relatively simple genetic algorithm yourself if you want to.

You can use any starting values you want to. The problem is that the best fit parameter values ae not necessarily unique, so an algorithm such as ga will search the entire parameter space for the best fit. If your objective function changes, or the data you are fitting change, you will need to search globally for the best parameter set for the new situation.

[Zuyu An]

Ich ändert die Funktion bischen:

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)./(2.*x(4).*(xdata+x(5)).*x(6).^2));

mach 2.*pi.*x(4) in ^3

and now it fitting good with linspace and fminsearch

ganz komisch, but still thank you very much!

[Star Strider]

As always, my pleasure!

Answer 2

[Alex Sha]

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));

[Bjorn Gustavsson]

This is the core of the problem. That the parameterization has the form - that all points on the curve have the same value often makes the optimization-functions to rumble around in the most peculiar paths in the parameter-space. The function-parameterization have another similar possible-problem: is in my experience better to implement as with different meaning of p(4) and p(5). Fixing the first should be the most important.

HTH

[Zuyu An]

It should be with first data(0, -0.7579)