Function with multiple input parameters to be determined through fitting

Hi, I have a function to be fitted to some experimental data, but this function has multiple fitting parameters (6 parameters). I know that one of the ways to find these fitting parameters is to pass the function, the inital guess to the parameters and the lower and upper bounds into the lsqcurvefit function.

I have already created the function, but it seems like the software doesn't know that the p(1),p(2), p(3)... are parameters. Sorry if the function is very long.

What I wanted was p to be a vector with 6 elements p(1) p(2) p(3) p(4) p(5) p(6) but it seems like the software cant differentiate betwen p and p(1) p(2) p(3) p(4) p(5) p(6). When I try running the code, the error message I got was "Not enough input arguments"

Thank you for helping

function F = EM_SS(p, E_p)
F = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*(int(sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + 126*p(5)^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), E, p(3), Inf, 'ArrayValued', 1)) + p(2)*(2*pi*p(4)^3/2)*1/p(6)*((1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6))); 

end

[Jack]

please refer to this picture for the function to be fitted. everything other than and π is a fitting parameter. is given by: where R is also a fitting parameter. There should be a total pf 6 fitting parameters.

Answer 1

[Torsten]

You have to use "integral" instead of "int" and loop over the elements in E_p:

EM_SS([1 1 1 1 1 1],1)

ans = 3.2940e+03

function F = EM_SS(p, e_p)
  for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
    (integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
    126*p(5)^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
    p(2)*(2*pi*p(4)^3/2)*1/p(6)*(...
    (1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
    (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
    (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
    (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
    (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
    (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
    (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6))); 
  end
end

[Jack]

will give it a try. Thank you!

[Jack]

Hi, I've given it a try, but the software still says that there is insufficient input arguments. Furthermore, the system also says that the variable "F(i)" appears to change size on every loop iteration.

My data are saved as "double" int he workspace

[Torsten]

Before calling lsqcurvefit, call the function "EM_SS" with your initial guesses for the parameters and your xdata vector for "hw" and tell us what the function returns.

When I call EM_SS with ([1 1 1 1 1 1],1), I get a proper output (see above).

I assumed your input "e_p" to the function "EM_SS" is a row vector. If this is not the case, you will have to change the F-vector according to the size of "e_p".

[Jack]

Hi below is the message that I got from the command window when I call EM_SS([1 1 1 1 1 1], 1) in the script that runs the lsqcurvefit.

for context, I have 2 seperate scripts: one runs the lsqcurvefit function (below) and the other runs the EM_SS function

May I ask what does the value 3.29...e+03 represents?

my input data (E_p) to EM_SS is a column vector of size 167 by 1. May I ask how do I change the F_vector according to the size of E_p? I'm sorry, I am quite new to this software

[Jack]

% Data files 
% global T_sample T_substrate Lambda E_p Abs
load Steady_State_Data.mat
h = 4.135667696*10^-15;
c = 3*10^8;
E_p = (Lambda*10^-9/h*c).^-1;
Abs = log10(T_substrate./T_sample);
lb = []; ub = [];
p0 = [0.13 0.1 1.6 0.05 3 1]; % initial guess for parameters (refer to line 13 for the order of parameters) 
% p0 = [A_1 A_2 E_g E_b R g]

% choose between different algorithm (3C1, comment those lines that are not choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt');
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective');
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point');
EM_SS([1 1 1 1 1 1], 1)

% solver
[p, fval, resnorm, residual, exitflag, jacobian] = lsqcurvefit(EM_SS, p0, E_p, Abs, lb, ub); 

plot(E_p, Abs) 

For reference, this is the script that runs the lsqcurefit function. and the EM_SS script is excatly similar to the one that you've editted for me.

T_sample, T_substrate and E_p are arrays of size 167 by 1

[Torsten]

You didn't include your data.

What is the size of your

E_p = (Lambda*10^-9/h*c).^-1

?

Did you notice that the function is defined as

function F = EM_SS(p, e_p)

and not as

function F = EM_SS(p, E_p)

?

And the call to lsqcurvefit must be changed from

[p, fval, resnorm, residual, exitflag, jacobian] = lsqcurvefit(EM_SS, p0, E_p, Abs, lb, ub);

to

[p, fval, resnorm, residual, exitflag, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub);

[Jack]

size of E_p is 167 by 1.

My apologies for missing out the data. Please refer to the attached.

yes I do realise that the functionis defined as

function F = EM_SS(p, e_p)

instead of

function F = EM_SS(p, E_p)

May I ask what is the difference between the 2? What effect does the capital letter have?

[Torsten]

May I ask what is the difference between the 2? What effect does the capital letter have?

I use e_p as the array of hw and E_p as the i-th value from the array.

This code works, but the parameter vector becomes complex in the course of the computation which makes the "integral" function stop:

% Data files 
% global T_sample T_substrate Lambda E_p Abs
load Steady_State_Data.mat
h = 4.135667696*10^-15;
c = 3*10^8;
E_p = (Lambda*10^-9/h*c).^-1;
Abs = log10(T_substrate./T_sample);
lb = []; ub = [];
p0 = [0.13 0.1 1.6 0.05 3 1]; % initial guess for parameters (refer to line 13 for the order of parameters) 
% p0 = [A_1 A_2 E_g E_b R g]

% choose between different algorithm (3C1, comment those lines that are not choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt');
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective');
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point');

% solver
[p, fval, resnorm, residual, exitflag, jacobian] = lsqcurvefit(@EM_SS, p0, E_p, Abs, lb, ub); 

Error using integralCalc (line 33)
Contour endpoints and waypoints must be finite.

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in solution>EM_SS (line 26)
(integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...

Error in lsqcurvefit/objective (line 321)
F = feval(funfcn_x_xdata{3},x,XDATA,varargin{:});

Error in snls (line 344)
newfvec = feval(funfcn{3},xcurr,varargin{:});

Error in lsqncommon (line 204)
snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,F,Jac,caller, ...

Error in lsqcurvefit (line 306)
lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,optimgetFlag,caller,...

plot(E_p, Abs) 

function F = EM_SS(p, e_p)
  for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
    (integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
    126*p(5)^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
    p(2)*(2*pi*p(4)^3/2)*1/p(6)*(...
    (1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
    (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
    (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
    (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
    (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
    (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
    (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6))); 
  end
  F = F(:);
end

[Jack]

Hi I have editted the code and ran it again, this time the error messgae is as follow:

[Jack]

Is there a way to overcome the complexity of the parameter vector and make the integral function run?

[Jack]

Hi I re-ran the code with the new addition of F = F(:) but it seems like the loop does not terminate and just keeps running.

[Torsten]

The code now works technically, but you will get many, many problems with the attempt to fit so many unknown parameters and with the complicated integral. So I'll leave it to you to solve the problem-immanent difficulties and maybe to reduce the number of unknowns.

[Jack]

I see I see. Any idea where I can find resources that will help resolve the possible issues that I may face?

[Torsten]

I think the problems will be too problem-specific to find general ressources. As said: less parameters and good initial guesses for them would be helpful. Did you plot Abs against E_p for the initial guesses you provided to get an impression of the "initial fitting quality" ?

[Jack]

I have tried doing the scatter plot of Abs against E_p but I have been trying for days to plot the curve that will fit the scatter plot.