Problem fitting the curve with fminsearch
2 views (last 30 days)
Show older comments
I have been trying to understand how fminsearch works or the result it gives me for at least a while (it doesnt fit my curve). I am trying to fit the curve but when I set the function to be in polynomial or sinusoidal form it does not fit the curve to the experimental like it does with the cftool.
One of the codes I am using is the following:
theta_i = -29:1:30;
theta_i = theta_i';
x_b = 3.85802467750551e-09 1.23456782463727e-07 9.37499560538235e-07 3.95060948021886e-06 1.20562544833058e-05 2.99995500044892e-05 6.48397188022232e-05 0.000126411762446210 0.000227786552702836 0.000385728056932932 0.000621145743629370 0.000959539347420657 0.00143143208410479 0.00207278707531322 0.00292540015345999 0.00403726036144925 0.00546286733815604 0.00726349240257385 0.00950736754606585 0.0122697837643895 0.0156330772872862 0.0196864794049803 0.0245258028990383 0.0302529357565231 0.0369751111287854 0.0448039216911671 0.0538540470230045 0.0642416647639122 0.0760825205564049 0.0894896386199658 0.104570664671209 0.121424846210700 0.140139672234976 0.160787215323602 0.183420243663010 0.208068198426751 0.234733162130625 0.263385974705427 0.293962684097702 0.326361544655273 0.360440796277300 0.396017466794855 0.432867435598924 0.470726974599819 0.509295940232559 0.548242725593590 0.587210994794329 0.625828114407609 0.663715074344481 0.700497560449578 0.735817714176249 0.769346003776893 0.800792550572930 0.829917216722031 0.856537778826176 0.880535591658960 0.901858288833670 0.920519265282326 0.936593924749780 0.950212931632136 ; % I have transposed the data just to copy it here so that it does not occupy much, it is a 60x1 matrix.
% if it were to be plotted it would look like plot(theta_i,x_b), the theta_i vector is the horizontal axis and the x_b the vertical.
IG = ones(1,6);
dF = [theta_i,x_b];
x33 = dF(:,1); %same as theta_i
y33 = dF(:,2); %same as x_b
P = fit2(dF,IG);
a1 = P(1);
a2 = P(2);
a3 = P(3);
a4 = P(4);
a5 = P(5);
a6 = P(6);
yfit = a1.*x33.^5 + a2.*x33.^4 + a3.*x33.^3 + a4.*x33.^2 + a5.*x33 + a6;
plot(theta_i,x_b);
hold on
plot (theta_i,yfit);
hold off
%% And the function is the next one:
function P = fit2(dF,IG)
x33 = dF(:,1);
y33 = dF(:,2);
function [ds] = fit(IG)
a1 = IG(1);
a2 = IG(2);
a3 = IG(3);
a4 = IG(4);
a5 = IG(5);
a6 = IG(6);
f33 = a1.*x33.^5 + a2.*x33.^4 + a3.*x33.^3 + a4.*x33.^2 + a5.*x33 + a6;
ds = (f33-y33).^2;
ds = sum(ds);
end
P = fminsearch (@fit,IG);
end
%maybe there are sombe bugs because y have edited a little bit for this post, but the idea i want to express i think its clear
0 Comments
Accepted Answer
Fabio Freschi
on 9 Oct 2023
Edited: Fabio Freschi
on 9 Oct 2023
I would use lsqnonlin: it's more general than polyfit suggested by @Torsten (that in this case would work perfectly). Qualitatively, it seems that the resulting fitting is good. In addition I have cleaned your code a little.
% your data
theta_i = (-29:1:30).';
x_b = [3.85802467750551e-09 1.23456782463727e-07 9.37499560538235e-07 3.95060948021886e-06 1.20562544833058e-05 2.99995500044892e-05 6.48397188022232e-05 0.000126411762446210 0.000227786552702836 0.000385728056932932 0.000621145743629370 0.000959539347420657 0.00143143208410479 0.00207278707531322 0.00292540015345999 0.00403726036144925 0.00546286733815604 0.00726349240257385 0.00950736754606585 0.0122697837643895 0.0156330772872862 0.0196864794049803 0.0245258028990383 0.0302529357565231 0.0369751111287854 0.0448039216911671 0.0538540470230045 0.0642416647639122 0.0760825205564049 0.0894896386199658 0.104570664671209 0.121424846210700 0.140139672234976 0.160787215323602 0.183420243663010 0.208068198426751 0.234733162130625 0.263385974705427 0.293962684097702 0.326361544655273 0.360440796277300 0.396017466794855 0.432867435598924 0.470726974599819 0.509295940232559 0.548242725593590 0.587210994794329 0.625828114407609 0.663715074344481 0.700497560449578 0.735817714176249 0.769346003776893 0.800792550572930 0.829917216722031 0.856537778826176 0.880535591658960 0.901858288833670 0.920519265282326 0.936593924749780 0.950212931632136].'; % I have transposed the data just to copy it here so that it does not occupy much, it is a 60x1 matrix.
% fitting function
yfit = @(a)a(1)*theta_i.^5+a(2)*theta_i.^4+a(3)*theta_i.^3+a(4)*theta_i.^2+a(5)*theta_i+a(6);
% function to minimize
fun = @(a)yfit(a)-x_b;
% initial values
x0 = ones(6,1);
x = lsqnonlin(fun,x0);
% plot
figure, hold on
plot(theta_i,x_b);
plot (theta_i,yfit(x));
4 Comments
Fabio Freschi
on 9 Oct 2023
Fabio Freschi
on 9 Oct 2023
Moved: Matt J
on 9 Oct 2023
I changed the lsqnonlin with fminsearch. I had to play a little with options because I started from a blind initial guess
theta_i = (-29:1:30).';
x_b = [3.85802467750551e-09 1.23456782463727e-07 9.37499560538235e-07 3.95060948021886e-06 1.20562544833058e-05 2.99995500044892e-05 6.48397188022232e-05 0.000126411762446210 0.000227786552702836 0.000385728056932932 0.000621145743629370 0.000959539347420657 0.00143143208410479 0.00207278707531322 0.00292540015345999 0.00403726036144925 0.00546286733815604 0.00726349240257385 0.00950736754606585 0.0122697837643895 0.0156330772872862 0.0196864794049803 0.0245258028990383 0.0302529357565231 0.0369751111287854 0.0448039216911671 0.0538540470230045 0.0642416647639122 0.0760825205564049 0.0894896386199658 0.104570664671209 0.121424846210700 0.140139672234976 0.160787215323602 0.183420243663010 0.208068198426751 0.234733162130625 0.263385974705427 0.293962684097702 0.326361544655273 0.360440796277300 0.396017466794855 0.432867435598924 0.470726974599819 0.509295940232559 0.548242725593590 0.587210994794329 0.625828114407609 0.663715074344481 0.700497560449578 0.735817714176249 0.769346003776893 0.800792550572930 0.829917216722031 0.856537778826176 0.880535591658960 0.901858288833670 0.920519265282326 0.936593924749780 0.950212931632136].'; % I have transposed the data just to copy it here so that it does not occupy much, it is a 60x1 matrix.
% if it were to be plotted it would look like plot(theta_i,x_b), the theta_i vector is the horizontal axis and the x_b the vertical.
% fitting function
yfit = @(a)a(1)*theta_i.^5+a(2)*theta_i.^4+a(3)*theta_i.^3+a(4)*theta_i.^2+a(5)*theta_i+a(6);
% function to minimize
fun = @(a)sum((yfit(a)-x_b).^2);
% initial values
x0 = ones(6,1);
% options
options = optimset('MaxFunEvals',1e6,'MaxIter',1e6,'TolX',1e-16,'TolFun',1e-16);
% fminsearch
[x,fval,exitflag,output] = fminsearch(fun,x0,options);
% plot
figure, hold on
plot(theta_i,x_b);
plot (theta_i,yfit(x));
More Answers (3)
Matt J
on 9 Oct 2023
Edited: Matt J
on 9 Oct 2023
but when I set the function to be in polynomial or sinusoidal form it does not fit the curve to the experimental like it does with the cftool
Are you sure you wouldn't rather use a Gaussian curve model? It seems much more appropriate, and is easily obtained by downloading gaussfitn:
theta_i = -29:1:30;
x_b = [3.85802467750551e-09 1.23456782463727e-07 9.37499560538235e-07 3.95060948021886e-06 1.20562544833058e-05 2.99995500044892e-05 6.48397188022232e-05 0.000126411762446210 0.000227786552702836 0.000385728056932932 0.000621145743629370 0.000959539347420657 0.00143143208410479 0.00207278707531322 0.00292540015345999 0.00403726036144925 0.00546286733815604 0.00726349240257385 0.00950736754606585 0.0122697837643895 0.0156330772872862 0.0196864794049803 0.0245258028990383 0.0302529357565231 0.0369751111287854 0.0448039216911671 0.0538540470230045 0.0642416647639122 0.0760825205564049 0.0894896386199658 0.104570664671209 0.121424846210700 0.140139672234976 0.160787215323602 0.183420243663010 0.208068198426751 0.234733162130625 0.263385974705427 0.293962684097702 0.326361544655273 0.360440796277300 0.396017466794855 0.432867435598924 0.470726974599819 0.509295940232559 0.548242725593590 0.587210994794329 0.625828114407609 0.663715074344481 0.700497560449578 0.735817714176249 0.769346003776893 0.800792550572930 0.829917216722031 0.856537778826176 0.880535591658960 0.901858288833670 0.920519265282326 0.936593924749780 0.950212931632136 ] ;
params=gaussfitn(theta_i',x_b');
[D,A,mu,sig]=deal(params{:});
x = D + A*exp( -0.5 * (theta_i-mu).^2/sig );
plot(theta_i,x_b,'x',theta_i,x); legend('Data', 'Gauss Fit','Location','SouthEast')
1 Comment
Matt J
on 9 Oct 2023
Edited: Matt J
on 9 Oct 2023
With fminsearch:
theta_i = -29:1:30;
x_b = [3.85802467750551e-09 1.23456782463727e-07 9.37499560538235e-07 3.95060948021886e-06 1.20562544833058e-05 2.99995500044892e-05 6.48397188022232e-05 0.000126411762446210 0.000227786552702836 0.000385728056932932 0.000621145743629370 0.000959539347420657 0.00143143208410479 0.00207278707531322 0.00292540015345999 0.00403726036144925 0.00546286733815604 0.00726349240257385 0.00950736754606585 0.0122697837643895 0.0156330772872862 0.0196864794049803 0.0245258028990383 0.0302529357565231 0.0369751111287854 0.0448039216911671 0.0538540470230045 0.0642416647639122 0.0760825205564049 0.0894896386199658 0.104570664671209 0.121424846210700 0.140139672234976 0.160787215323602 0.183420243663010 0.208068198426751 0.234733162130625 0.263385974705427 0.293962684097702 0.326361544655273 0.360440796277300 0.396017466794855 0.432867435598924 0.470726974599819 0.509295940232559 0.548242725593590 0.587210994794329 0.625828114407609 0.663715074344481 0.700497560449578 0.735817714176249 0.769346003776893 0.800792550572930 0.829917216722031 0.856537778826176 0.880535591658960 0.901858288833670 0.920519265282326 0.936593924749780 0.950212931632136 ] ;
fun=@(p) mdlFcn(p,theta_i, x_b);
p=fminsearch(fun, [15,30]);
[~,DA,x]=fun(p);
plot(theta_i,x_b,'x',theta_i,x); legend('Data', 'Gauss Fit','Location','SouthEast')
function [cost,DA,x]=mdlFcn(p,theta, x_b)
g=exp( -0.5 * (theta(:)-p(1)).^2/p(2) );
G=[g.^0,g];
DA=G\x_b(:);
x=G*DA; %predicted curve
cost=norm(x-x_b(:));
end
Xabier Tamayo
on 9 Oct 2023
3 Comments
Fabio Freschi
on 9 Oct 2023
Edited: Fabio Freschi
on 9 Oct 2023
Nice trick to make fminsearch "converge"...
Torsten
on 9 Oct 2023
See Also
Categories
Find more on Interpolation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!