Non-linear data fitting to a model
1 view (last 30 days)
Show older comments
I am trying to fit some data using the following code:
xdata = [0.0600 0.1250 0.2540 0.4050 0.5450];
ydata = [0.0800 0.1000 0.1200 0.1400 0.1600];
objfcn = @(x)(x(1)*x(2))./((x(2)*(1-xdata)./xdata)+(1-xdata)*(x(2)-1))-ydata;
x0 = rand(2,1);
opts = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','Display','iter');
[x,resnorm]=lsqnonlin(objfcn,x0,[],[],opts);
M = x(1)
B = x(2)
I am getting results but sometimes M value is coming negative which is not desirable. Also, the predicted ydata deviates largely from the original ydata values. If someone could be help to analyze the small code mentioned above and what is the wrong with results.
Thanks,
0 Comments
Accepted Answer
Matt J
on 3 Nov 2014
Edited: Matt J
on 3 Nov 2014
Random guesses for the initial point x0 are a bad idea. I generate a better initial guess below, by re-arranging the model equation as linear. Also, if it is undesirable for M to be negative, then you should use the lb argument to bound it. This means you will then have to use the trust-region algorithm instead of Levenberg-Marquardt, since the latter doesn't support bound constraints.
xdata=xdata(:); ydata=ydata(:);
objfcn = @(x)(x(1)*x(2))./((x(2)*(1-xdata)./xdata)+(1-xdata)*(x(2)-1))-ydata;
expr1=(1-xdata)./xdata;
expr2=(1-xdata);
A=[ones(size(xdata)), -(expr1+expr2).*ydata, ];
b=ydata.*expr2;
s=A\b;
res=A*s-b
x0 = [s(1)/s(2),s(2)],
opts = optimoptions(@lsqnonlin,'Display','iter');
[xf,resnorm]=lsqnonlin(objfcn,x0,[0,-inf],[],opts);
This leads to
K>> yfit= objfcn(xf)+ydata
yfit =
0.0827
0.0903
0.1067
0.1342
0.1757
I tend to think this is the best fit that will be possible given such a small number of xdata, ydata samples.
2 Comments
Matt J
on 3 Nov 2014
Edited: Matt J
on 3 Nov 2014
Generating the initial guess is a problem-specific art and not a science, I'm afraid. A frequent strategy, though, is to ignore measurement noise and then re-arrange or transform the model equation to something simpler. When we ignore noise in your original model, it can be re-arranged into the equation
(x(1)*x(2))=((x(2)*(1-xdata)./xdata)+(1-xdata)*(x(2)-1)).*ydata
Making the change of variables s(1)=x(1)*x(2), s(2)=x(2) this becomes linear in s
s(1)-((s(2)*(1-xdata)./xdata)+(1-xdata)*(s(2)-1)).*ydata = 0
and can be solved using standard linear solvers, which is what I did for you in my code above.
Another way to generate an initial guess, because you only have 2 unknown variables, is to make a surf plot of the resnorm function and inspect it for the approximate location of the minimum. Or maybe instead of plotting, just sample the resnorm function on some coarse sampling grid and use min() to locate an approximate global minimum.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!