Problem fitting the curve with fminsearch

2 vues (au cours des 30 derniers jours)
Xabier Tamayo
Xabier Tamayo le 9 Oct 2023
Modifié(e) : Matt J le 9 Oct 2023
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

Réponse acceptée

Fabio Freschi
Fabio Freschi le 9 Oct 2023
Modifié(e) : Fabio Freschi le 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 commentaires
Fabio Freschi
Fabio Freschi le 9 Oct 2023
@Torsten you are right... lsqnonlin is an automatic choice for me as first attempt with fitting
Fabio Freschi
Fabio Freschi le 9 Oct 2023
Déplacé(e) : Matt J le 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));

Connectez-vous pour commenter.

Plus de réponses (3)

Torsten
Torsten le 9 Oct 2023
Modifié(e) : Torsten le 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 le 9 Oct 2023
Modifié(e) : Matt J le 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 commentaire
Matt J
Matt J le 9 Oct 2023
Modifié(e) : Matt J le 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

Connectez-vous pour commenter.


Xabier Tamayo
Xabier Tamayo le 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 commentaires
Fabio Freschi
Fabio Freschi le 9 Oct 2023
Modifié(e) : Fabio Freschi le 9 Oct 2023
Nice trick to make fminsearch "converge"...
Torsten
Torsten le 9 Oct 2023
You should take @Fabio Freschi 's solution. It's more ... serious.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Interpolation dans Help Center et File Exchange

Produits


Version

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by