Problem fitting the curve with fminsearch
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
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);
% plot
figure, hold on
plot(theta_i,x_b);
plot (theta_i,yfit(x));
4 commentaires
Fabio Freschi
le 9 Oct 2023
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));
Plus de réponses (3)
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');
[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
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
Xabier Tamayo
le 9 Oct 2023
3 commentaires
Fabio Freschi
le 9 Oct 2023
Modifié(e) : Fabio Freschi
le 9 Oct 2023
Nice trick to make fminsearch "converge"...
Torsten
le 9 Oct 2023
Voir également
Catégories
En savoir plus sur Interpolation dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!