nlinfit seems to find only local minimum

Hello!

I am busy fitting my experimental data with this code:

function dCCconv = DiffEqSolverALLSEandSO(param, t, k_1n, k_1p, mu_n, phi, FA, WL, lf, thickness, T); 
RC = 18e-9;
k_2 = param(1);
k_3 = param(2);

filename1 = 'pulse300.dat'
pathname1 = 'D:\My Documents\Experiments\TRMC\PulseShapes\';
pulse = importfile1([pathname1, filename1]);

   assignin('base', 't',t);
   assignin('base', 'RC',RC);

   [P,I] = max(pulse(:,1)); %find the coordinates of the maximum value of t
   [P1, I1] = find(pulse(:,2)>=0);
  pulse_begin_time = pulse(I1(1),1);
  pulse_begin = min(find(t>=pulse_begin_time));
  pulse_end = max(find(t<=P));

[P3, I3] = max(t);
dx = t(2)-t(1);

t_x = t(pulse_begin : pulse_end);
pulses(:,2)=spline(pulse(:,1), pulse(:,2), t_x); % make the pulse and the experimental data have the same time scale
pulses(end+1:I3,2)= 0;
pulses(:,1) = t;
assignin('base', 'pulses',pulses);
pmax = trapz(pulses(:,1), pulses(:,2));

NA = lf*1e-3*WL*1e-9*FA/(6.626e-34*3e8*2.1); % number of absorbed photons per cm^2 (2.1 is in cm^2)
NP = lf*1e-3*WL*1e-9/(6.626e-34*3e8*2.1); % number of incident photons per cm^2

G0 = NA/(pmax*thickness*1e2);

function dx = myode(t, x, pulses, G0, k_1n, k_1p, k_2,k_3, phi, T)  %k_3,phi_n,k_1p,k_2,
    dx = zeros(2,1); 
pulse_actual = interp1(pulses(:,1),pulses(:,2),t );
SE_n = (((T/497))*real(abs(k_1n*1e7))* (t*real(abs(k_1n*1e7)))^((T/497)-1) );
SE_p = ((T/497)*real(abs(k_1p*1e6))* (t*real(abs(k_1p*1e6)))^((T/497)-1) );
dx(1) = real((G0 * pulse_actual*real(abs(phi))) - SE_n * real(x(1)) - abs(k_2*1e-10) * x(1) * x(2))- real(abs(k_3)) * 1e-28 * (x(1)^2 * x(2) + x(1) * x(2)^2) ; 
dx(2) = real((G0 * pulse_actual*real(abs(phi))) - SE_p * real(x(2)) - abs(k_2*1e-10) * x(1) * x(2))- real(abs(k_3)) * 1e-28 * (x(1)^2 * x(2) + x(1) * x(2)^2);  
end

tspan = [t(2): dx : t(end)];
opts = odeset('RelTol',1e-5,'AbsTol',1e-6);
ic = [0;0];
[t,x] = ode45(@(t,x) myode(t, x, pulses, G0, k_1n, k_1p,  k_2, k_3, phi, T), tspan, ic, opts);
CCx = x;

f1x = (real(abs(mu_n)) * CCx(:,1)+ (1/2.06)*real(abs(mu_n))*CCx(:,2))*thickness*1e2 / NP; 
f2x = conv(exp(-t/RC),f1x); %convolution to take into acount of the time response of the cavity cell
tr_cav = sum(exp(-t/RC));
f(1:I3) = 100*f2x(1:I3)/tr_cav; 
dCCconv = f;

end

using nlinfit:

[xfit,resnorm, Jacob, CovB, MSE] = nlinfit(  handles.timecorr, 100*handles.datacorr',@(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature) , handles.x0);

Doing this finds usually a local minimum and I vary the guess parameters until I get a minimum for a chi2 (by hand). Now I tried to fit using lsqcurvefit (just to see if this finds the absolute minimum) using this:

xfit,resnorm] = lsqcurvefit( @(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature), handles.x0, handles.timecorr, 100*handles.datacorr',handles.lb, handles.ub, fitopt ); %

But lsqcurvefit runs infinitely long without returning any result. Is there a better fitting function for my case, or is there something in general wrong? I am not sure how I could improve my fit.

Thanks for your help!

Answer 1

[Matt J]

Make sure you follow guidelines for Optimizing a Simulation or Ordinary Differential Equation. If you continue to have difficulty finding global solutions, you should perhaps consider MultiStart or other tools in the Global Optimization Toolbox.

[Silke]

Hi Matt,

could you please be more specific? Which parameter described in the guideline would you start varying?

Thank you!

[Matt J]

Well, all of the material mentioned there could be relevant to you. However, a direct search method, as described at

https://www.mathworks.com/help/optim/ug/optimizing-a-simulation-or-ordinary-differential-equation.html#btfixp4

would be my first choice. Since you have only 2 unknowns, fminsearch might be sufficient.

[Silke]

Are you suggesting to use fminsearch instead of nlinfit? Or how do I implement fminsearch in my code?

[Matt J]

Yes I am suggesting fminsearch instead of nlinfit.

[Silke]

I have tried using fminsearch instead of nlinfit, but I also do not get a better result.

This image shows a calculation where both variables are set to 2. I would say that this calculation represents the data quite well.

Fitting the data with fminsearch and using 5 and 2 as starting parameters gives this solution:

And this is what nlinfit finds with 5 and 2 as starting parameters.

I do not think that a guess value of 5 and 2 is too far off, when 2 and 2 give good values. Does anyone have any idea how to solve this? Is the problem the ODE solver, or the fitting function?

[Torsten]

The peak is short, the tail is long. As a consequence, the data are fitted such that the tail is approximated very well while the peak is neglected. Either cut off the tail earlier or double/triple the data for the peak.

Best wishes

Torsten.

[Silke]

Just by interpolating more data points in the peak of the measurement data?

[Torsten]

No interpolation required, just double or triple the existing data.

[Matt J]

I do not think that a guess value of 5 and 2 is too far off, when 2 and 2 give good values.

The only way to be sure of that is to use 2 and 2 as the initial guess and see how the iterations modify things.

[Silke]

If I start with 2 and 2 as guess values, the result is 2 and 2 as well. It does not change at all.

[Matt J]

And that confirms that your previous initial guess (5,2) was not sufficiently close.

[Silke]

But how close does my initial guess have to be to find a good solution? If I would already know the parameters for all my experiments, I would not need to perform the fit. I am quite frustrated now that the fitting does not work better at the moment. Even if I start off with 3 and 2, the fit does not find the good parameters.

[Matt J]

It is never apriori clear how close the initial guess has to be, but that is why I recommended the Global Optimization Toolbox for cases (like yours possibly) where it seems overly difficult.

Alternatively, because you only have two unknown parameters, it might be computationally affordable for you to do a sweep of the parameters and plot the objective function as a 2D surface. From this you could see, approximately, where the global minimum lies and get a more informed initial guess.

Answer 2

[Silke]

If I triple the data in the peak, but not in the tail, I will not have a continous dataset anymore, won't I? And in how far will this change my fitting parameters? I am not sure if I fully understand your suggestion.

[Torsten]

The parameters during the fitting process are determined such that the sum of squared differences between experimental and simulated data is minimized.

If you use more data in a certain region of the independent variable (I think it's "time" in your case), the contribution to the sum of squared differences will increase in this region thus forcing the fit to be more exact there.

Best wishes

Torsten.

[Silke]

Okay, thank you for your reply. If I understand you correctly now, you mean that I have more datapoints in the peak, so not doubling the actual value of the peak. So for this I need to interpolate my experimental data, so that I have a smaller time-step in the peak than in the tail. Does the ODE solver work with non-constant time-spacing?

[Torsten]

You don't need to change anything with the ODE solver or do interpolation.

If your data now look like

(t1,y1),(t2,y2),(t3,y3),(t4,y4),(t5,y5),(t6,y6),(t7,y7)

and t3, t4, t5 and t6 lie in the tail, replace the above dataset by

(t1,y1),(t1,y1),(t1,y1),(t2,y2),(t2,y2),(t2,y2),(t3,y3),(t4,y4),(t5,y5),(t6,y6),(t7,y7)

Best wishes

Torsten.

[Torsten]

Or use different weights for your data dependent on the region of the independent variable.

Look at the section

"Nonlinear Regression using observation weights"

under

https://de.mathworks.com/help/stats/nlinfit.html

Best wishes

Torsten.

[Silke]

Well, I have tried now both of your suggestions. First, I have tripled the datapoints in the peak region. This gives an error message in this line of myode:

pulse_actual = interp1(pulses(:,1),pulses(:,2),t );

because for the interpolation the time needs to increase monotonically.

Including a weight did for some reason not change anything in the fit. Even though I used a factor of 10000 for the data at the peak and 1 in the tail. I have no clue why this has no influence.

Do you have any other suggestion on how to solve this?

[Torsten]

I think prescribing the weights is the easiest way to give emphasis to the peak approximation. You should again check why setting weights in your call of "nlinfit" doesn't influence the final result for the parameters being estimated.

Best wishes

Torsten.

[Silke]

How can I check this? I varied them, but they do not seem to have any influence. I included them like that:

 W(1:49) = 5;
        W(50:64) = 1000;
        W(65) = 10000000;
        W(66:80) = 1000;
        W(81:500) = 5;
        W(501:3600) = 0.000001;
 [xfit,resnorm, Jacob, CovB, MSE] = nlinfit( handles.timecorr, 100*handles.datacorr',@(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature, sel2) , handles.x0, 'Weights', W);

[Torsten]

From the documentation:

Weight function for robust fitting, specified as a character vector indicating a weight function, or a function handle. WgtFun is valid only when Robust has value 'on'.

So additionally, you will have to use the "Robust fit Option", I guess.

And maybe you should allocate W as a column vector first:

W = zeros(3600,1);

Best wishes

Torsten.

[Silke]

I did it like that:

W = zeros(3600,1);
        W(1:49) = 5;
        W(50:64) = 1000;
        W(65) = 10000000;
        W(66:80) = 1000;
        W(81:500) = 5;
        W(501:3600) = 0.000001;
         assignin('base', 'W',W); 
         opts = statset('nlinfit');
        opts.RobustWgtFun = [];
  [xfit,resnorm, Jacob, CovB, MSE] = nlinfit( handles.timecorr, 100*handles.datacorr',@(param,timecorr) DiffEqSolverALLSEandSO(param, timecorr, handles.K_1N_GUESS_VALUE, handles.MU_N_GUESS_VALUE, handles.PHI_GUESS_VALUE, handles.K_1P_GUESS_VALUE, FA, WL, lf, thickness, temperature, sel2) , handles.x0,  opts, 'Weights', W');

But it seems not to change anything on the fit. Is this not correct?

[Torsten]

If the estimated parameters remain the same, you can check whether a weight matrix is used at all by looking at the output variable "resnorm". resnorm(i) should be sqrt(w(i))*res(i) where "res" is the output without using the weighting option.

Best wishes

Torsten.