Curve fitting using lsqcurvefit

Hi currently i am trying to fit my measured data to a model. I have followed this example https://au.mathworks.com/matlabcentral/answers/478835-lsqcurvefit-initial-values-stays-the-same#answer_391692 . I am also attaching the code from the example. i want to know what parameters do i need to change to fit my measured data to the model and extract the values from it. Any help is appreicated. I am also attaching my data.

T = readtable('cole.xlsx'); 
x = T{:,1}; %freq
y = T{:,2}; %real
%Transpose
freq = x.';
e_real = y.';

options = optimset('MaxFunEvals',10000);
options=optimset(options,'MaxIter',10000);
guess = 4.5;
UB = guess + 1;
LB = guess - 1;
lb = [];
ub = [];
x0 = [6,guess,2.5,0.6,0.09];
x  = lsqcurvefit(@flsq,x0,freq,e_real,lb,ub,options)

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.

x = 
   4.8828 - 0.0885i   5.0695 - 0.1843i   1.2278 -92.9851i  -1.0163 - 9.3958i  -0.0126 - 0.0739i

e_f = x(1);
e_del = x(2)*1e2;
tau1 = x(3)*1e-12;
alf1 = x(4);
sig = x(5);
yfit = real(flsq(x,freq));
%plot e_real against freq
plot(freq,e_real,'k.',freq,yfit,'b-')
legend('Data','Fitted')
title('Data and Fitted Curve')

function y = flsq(x,freq)
    x(3)=x(3)*1e-12;
    y = x(1) + (x(2)-x(1))./(1 + ((1j*2*pi.*freq*x(3)).^(1-x(4))))+x(5)./(1j*2*pi*freq*8.854e-12);
    
end

[Walter Roberson]

I am not sure what you want done different from what you already have?

[Kabit Kishore]

Hi Walter thank you. I want the fitted curve to go through maximum number of measured data points. Is it possible.

[Image Analyst]

That's a pretty complicated model. Would you rather just fit it to like 3-6 Gaussians instead? Or does that model have some theoretical significance based on the physics of the situation?

[Image Analyst]

@Kabit Kishore, any response to my questions in my last comment above? From Walter's last comment it looks like the model is no good so why use that particular model?

[Kabit Kishore]

Hi @Image Analyst. Yes the model that i am trying to fit has theoretical significance. Curently i am guessing the values which make it pretty dificult to fit the equation to my data. Is there a iterative automatic way which updates the values and provides the best fit. Thank you for your response.

[Walter Roberson]

lsqcurvefit() and fmincon() both do iterative automatic ways to update values to provide best fits.

[Walter Roberson]

If this has physical significance, then are there any constraints on the parameters of the model? Any that have to be real-valued? Any that have to be positive?

Answer 1

[Walter Roberson]

You can get a better fit by either telling lsqcurvefit() a better starting point, or by using a different optimizer.

But... remember that "better fit" mathematically does not necessarily mean "follows the curve more closely"

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/917874/cole.xlsx';
T = readtable(filename); 
x = T{:,1}; %freq
y = T{:,2}; %real
%Transpose
freq = x.';
e_real = y.';

options = optimset('MaxFunEvals',10000);
options=optimset(options,'MaxIter',10000);
guess = 4.5;
UB = guess + 1;
LB = guess - 1;
lb = [];
ub = [];
x0 = [6,guess,2.5,0.6,0.09];
[x, fval]  = lsqcurvefit(@flsq,x0,freq,e_real,lb,ub,options)

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.

x = 
           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

fval = 
          2.50307155416643

e_f = x(1);
e_del = x(2)*1e2;
tau1 = x(3)*1e-12;
alf1 = x(4);
sig = x(5);
yfit = real(flsq(x,freq));
%plot e_real against freq
plot(freq,e_real,'k.',freq,yfit,'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- lsqcurvefit')

Residue = @(x) norm(flsq(x,freq) - e_real);

options = optimset('MaxFunEvals', 1E5);
options = optimset(options, 'MaxIter', 1E5);
options = optimset(options, 'Display', 'none');

for K = 2 : 50
  guess = x0 + randn(1,5);
  [x(K,:), fval(K)]  = fmincon(Residue, guess, [], [], [], [], [], [], [], options);
end

[~, idx] = min(abs(fval));

disp('original -- lsqcurvefit')

original -- lsqcurvefit

disp([fval(1), x(1,:)])

  Columns 1 through 4

           2.50307155416643 +                     0i           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i

  Columns 5 through 6

          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

disp('best found in several tries of fmincon')

best found in several tries of fmincon

disp([fval(idx), x(idx,:)])

          1.82932684842456          4.25784166349341         -638714.651612716          5141.92675473731         -3.02543789002661      0.000633983073703346

e_f = x(idx,1);
e_del = x(idx,2)*1e2;
tau1 = x(idx,3)*1e-12;
alf1 = x(idx,4);
sig = x(idx,5);
yfit = real(flsq(x(idx,:),freq));
%plot e_real against freq
plot(freq, e_real, 'k.', freq, yfit, 'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- fmincon')

function y = flsq(x,freq)
    x(3)=x(3)*1e-12;
    y = x(1) + (x(2)-x(1))./(1 + ((1j*2*pi.*freq*x(3)).^(1-x(4))))+x(5)./(1j*2*pi*freq*8.854e-12);    
end

[Walter Roberson]

Sometimes the model is just not a good fit.

[Walter Roberson]

An improved run

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/917874/cole.xlsx';
T = readtable(filename); 
x = T{:,1}; %freq
y = T{:,2}; %real
%Transpose
freq = x.';
e_real = y.';

options = optimset('MaxFunEvals',10000);
options=optimset(options,'MaxIter',10000);
guess = 4.5;
UB = guess + 1;
LB = guess - 1;
lb = [];
ub = [];
x0 = [6,guess,2.5,0.6,0.09];
[x, fval]  = lsqcurvefit(@flsq,x0,freq,e_real,lb,ub,options)

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.

x = 
           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

fval = 
          2.50307155416643

e_f = x(1);
e_del = x(2)*1e2;
tau1 = x(3)*1e-12;
alf1 = x(4);
sig = x(5);
yfit = real(flsq(x,freq));
%plot e_real against freq
plot(freq,e_real,'k.',freq,yfit,'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- lsqcurvefit')

Residue = @(x) norm(flsq(x,freq) - e_real);

options = optimset('MaxFunEvals', 1E5);
options = optimset(options, 'MaxIter', 1E5);
options = optimset(options, 'Display', 'none');

for K = 2 : 50
  guess = randn(1,5)*100

;
  [x(K,:), fval(K)]  = fmincon(Residue, guess, [], [], [], [], [], [], [], options);
end

[~, idx] = min(abs(fval));

disp('original -- lsqcurvefit')

original -- lsqcurvefit

disp([fval(1), x(1,:)])

  Columns 1 through 4

           2.50307155416643 +                     0i           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i

  Columns 5 through 6

          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

disp('best found in several tries of fmincon')

best found in several tries of fmincon

disp([fval(idx), x(idx,:)])

          1.68283982412135          3.49149418809273           4.1928928642795          147.870041270639          44.9919207152756     -1.52002099046847e-05

e_f = x(idx,1);
e_del = x(idx,2)*1e2;
tau1 = x(idx,3)*1e-12;
alf1 = x(idx,4);
sig = x(idx,5);
yfit = real(flsq(x(idx,:),freq));
%plot e_real against freq
plot(freq, e_real, 'k.', freq, yfit, 'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- fmincon')

function y = flsq(x,freq)
    x(3)=x(3)*1e-12;
    y = x(1) + (x(2)-x(1))./(1 + ((1j*2*pi.*freq*x(3)).^(1-x(4))))+x(5)./(1j*2*pi*freq*8.854e-12);    
end

[Walter Roberson]

Slight improvement. I'm not sure you would be able to do much better -- this does 500 starting points.

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/917874/cole.xlsx';
T = readtable(filename); 
x = T{:,1}; %freq
y = T{:,2}; %real
%Transpose
freq = x.';
e_real = y.';

options = optimset('MaxFunEvals',10000);
options=optimset(options,'MaxIter',10000);
guess = 4.5;
UB = guess + 1;
LB = guess - 1;
lb = [];
ub = [];
x0 = [6,guess,2.5,0.6,0.09];
[x, fval]  = lsqcurvefit(@flsq,x0,freq,e_real,lb,ub,options)

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.

x = 
           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

fval = 
          2.50307155416643

e_f = x(1);
e_del = x(2)*1e2;
tau1 = x(3)*1e-12;
alf1 = x(4);
sig = x(5);
yfit = real(flsq(x,freq));
%plot e_real against freq
plot(freq,e_real,'k.',freq,yfit,'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- lsqcurvefit')

Residue = @(x) norm(flsq(x,freq) - e_real);

options = optimset('MaxFunEvals', 1E5);
options = optimset(options, 'MaxIter', 1E5);
options = optimset(options, 'Display', 'none');

for K = 2 : 500
  guess = rand(1,5)*100;
  [x(K,:), fval(K)]  = fmincon(Residue, guess, [], [], [], [], [], [], [], options);
end

[~, idx] = min(abs(fval));

disp('original -- lsqcurvefit')

original -- lsqcurvefit

disp([fval(1), x(1,:)])

  Columns 1 through 4

           2.50307155416643 +                     0i           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i

  Columns 5 through 6

          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

disp('best found in several tries of fmincon')

best found in several tries of fmincon

disp([fval(idx), x(idx,:)])

          1.66154803509977          3.41081977527717          4.19802178524268          149.532439789817          25.0006532849926      2.49463043656685e-06

e_f = x(idx,1);
e_del = x(idx,2)*1e2;
tau1 = x(idx,3)*1e-12;
alf1 = x(idx,4);
sig = x(idx,5);
yfit = real(flsq(x(idx,:),freq));
%plot e_real against freq
plot(freq, e_real, 'k.', freq, yfit, 'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- fmincon')

function y = flsq(x,freq)
    x(3)=x(3)*1e-12;
    y = x(1) + (x(2)-x(1))./(1 + ((1j*2*pi.*freq*x(3)).^(1-x(4))))+x(5)./(1j*2*pi*freq*8.854e-12);    
end

[Walter Roberson]

ga() is another iterative means for improving fit, but it does not help much. In one of my earlier runs, I got better than the below display, but it was still worse than fmincon() so I am not going to bother to show it.

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/917874/cole.xlsx';
T = readtable(filename); 
x = T{:,1}; %freq
y = T{:,2}; %real
%Transpose
freq = x.';
e_real = y.';

options = optimset('MaxFunEvals',10000);
options=optimset(options,'MaxIter',10000);
guess = 4.5;
UB = guess + 1;
LB = guess - 1;
lb = [];
ub = [];
x0 = [6,guess,2.5,0.6,0.09];
[x, fval]  = lsqcurvefit(@flsq,x0,freq,e_real,lb,ub,options)

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.

x = 
           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

fval = 
          2.50307155416643

e_f = x(1);
e_del = x(2)*1e2;
tau1 = x(3)*1e-12;
alf1 = x(4);
sig = x(5);
yfit = real(flsq(x,freq));
%plot e_real against freq
plot(freq,e_real,'k.',freq,yfit,'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- lsqcurvefit')

Residue = @(x) norm(flsq(x,freq) - e_real);

options = optimset('MaxFunEvals', 1E5);
options = optimset(options, 'MaxIter', 1E5);
options = optimset(options, 'Display', 'none');

tic
for K = 2 : 200
  [x(K,:), fval(K)]  = ga(Residue, 5, [], [], [], [], [], [], [], options);
end
toc

Elapsed time is 38.143890 seconds.

[~, idx] = min(abs(fval));

disp('original -- lsqcurvefit')

original -- lsqcurvefit

disp([fval(1), x(1,:)])

  Columns 1 through 4

           2.50307155416643 +                     0i           4.88278441116842 -    0.0884918454761463i           5.06952748610154 -     0.184297340006755i           1.22778921925455 -      92.9851457164006i

  Columns 5 through 6

          -1.01629693596011 -      9.39581846397004i        -0.0126281060006031 -     0.073855580530327i

disp('best found in several tries of ga')

best found in several tries of ga

disp([fval(idx), x(idx,:)])

           2.0261188723609          4.23046723004257          3.35326686725843          153.911545598562        -0.327809240216987        0.0453959038185919

e_f = x(idx,1);
e_del = x(idx,2)*1e2;
tau1 = x(idx,3)*1e-12;
alf1 = x(idx,4);
sig = x(idx,5);
yfit = real(flsq(x(idx,:),freq));
%plot e_real against freq
plot(freq, e_real, 'k.', freq, yfit, 'b-')
legend('Data','Fitted')
title('Data and Fitted Curve -- ga')

function y = flsq(x,freq)
    x(3)=x(3)*1e-12;
    y = x(1) + (x(2)-x(1))./(1 + ((1j*2*pi.*freq*x(3)).^(1-x(4))))+x(5)./(1j*2*pi*freq*8.854e-12);    
end

[Kabit Kishore]

Thank you @Walter Roberson. Much appreciated.