Problem fitting the curve with fminsearch

2 views (last 30 days)
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

Accepted Answer

Fabio Freschi
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);
Local minimum possible. lsqnonlin stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
% plot
figure, hold on
plot(theta_i,x_b);
plot (theta_i,yfit(x));
  4 Comments
Fabio Freschi
Fabio Freschi on 9 Oct 2023
@Torsten you are right... lsqnonlin is an automatic choice for me as first attempt with fitting
Fabio Freschi
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));

Sign in to comment.

More Answers (3)

Torsten
Torsten on 9 Oct 2023
Edited: Torsten on 9 Oct 2023
Using a polynomial of degree 5 to fit your data is a bad idea because the higher the degree, the greater the oscillations between the measurement data.
If you insist on your idea, use "polyfit" instead of your code.

Matt J
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');
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
[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
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

Sign in to comment.


Xabier Tamayo
Xabier Tamayo on 9 Oct 2023
The two ways that you have sent are very good and fast, but unfortunately they force me to use the fminsearch command. Any idea how to do it that way?
  3 Comments
Fabio Freschi
Fabio Freschi on 9 Oct 2023
Edited: Fabio Freschi on 9 Oct 2023
Nice trick to make fminsearch "converge"...
Torsten
Torsten on 9 Oct 2023
You should take @Fabio Freschi 's solution. It's more ... serious.

Sign in to comment.

Categories

Find more on Interpolation in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!