Problem fitting the curve with fminsearch

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

Answer 1

[Fabio Freschi]

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));

[Torsten]

@Fabio Freschi

Why a nonlinear fitting tool for a linear problem ?

[Xabier Tamayo]

Yeah that would be a great code to solve It. But they asked me to do with fminsearch. I dont need to use the 5 degree polynomial equation, but i have seen with cftool that It makes a good fit. I also have used a sinusoidal equation but the fit was not good also.

[Fabio Freschi]

@Torsten you are right... lsqnonlin is an automatic choice for me as first attempt with fitting

[Fabio Freschi]

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));

Answer 2

[Torsten]

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.

Answer 3

[Matt J]

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')

[Matt J]

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

Answer 4

[Xabier Tamayo]

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?

[Torsten]

:-)

format long
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   ].'    ; 
IG = [theta_i.^5, theta_i.^4, theta_i.^3, theta_i.^2, theta_i, ones(numel(theta_i),1) ]\x_b

IG = 6×1

-0.000000010602776 -0.000000452771947 0.000010318191427 0.000828248177965 0.014949620439656 0.090661449590186

dF = [theta_i,x_b];
x33 = dF(:,1); %same as theta_i
y33 = dF(:,2); %same as x_b
P = fit2(dF,IG)

P = 6×1

-0.000000010602776 -0.000000452771947 0.000010318191427 0.000828248177965 0.014949620439656 0.090661449590186

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

[Fabio Freschi]

Nice trick to make fminsearch "converge"...

[Torsten]

@Xabier Tamayo

You should take @Fabio Freschi 's solution. It's more ... serious.