# Nonlinear fitting depends too heavily on starting points

1 view (last 30 days)
Samuele Bolotta on 12 Mar 2021
Commented: Samuele Bolotta on 12 Mar 2021
Hi everyone,
I have a fitting procedure that fits a function to a current, generated this way:
function [EPSC, IPSC, CPSC, t] = generate_current(G_max_chl, G_max_glu, EGlu, EChl, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,tmax)
dt = 0.1; % time step duration (ms)
t = 0:dt:tmax-dt;
% Compute compound current
IPSC = ((G_max_chl) .* ((1 - exp(-t / tau_rise_In)) .* exp(-t / tau_decay_In)) * (Vm - EChl));
EPSC = ((G_max_glu) .* ((1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex)) * (Vm - EGlu));
CPSC = IPSC + EPSC;
end
This is a simulation, but later on the current will come from experimental recordings and therefore we won't know the underlying distribution (as in this case). In this script, I'm testing three different cases to see whether the procedure works well.
%% We have three situations. The first is the most realistic one; the second
% and the third are more and more challenging cases.
for i = 1:3
switch i
case 1
% Realistic values
[~,~,CPSC,t] = generate_current(80,15,0,-70,-50,0.44,15,0.73,3,120);
case 2
% More challenging situation: max values for conductances
% closer to each other
[~,~,CPSC,t] = generate_current(60,40,0,-70,-50,0.44,15,0.73,3,120); % Gs same
case 3
% More challenging situation: max values for conductances
% closer to each other and tau values closer to each other
[~,~,CPSC,t] = generate_current(60,40,0,-70,-50,0.44,8,0.73,3,120); % Gs + Tau same
end
% 5 iterations for each of the three cases
for z = 1:5
% Apply white noise to the CPSC - every time again, in order to have
% different values for the noise each iteration
y = awgn(CPSC,25,'measured');
%% Perform fit
[xData, yData] = prepareCurveData(t, y);
% Set up fittype and options.
ft = fittype( '((G_max_chl) .* ((1 - exp(-t / tau_rise_In)) .* exp(-t / tau_decay_In)) * (Vm - EChl)) + ((G_max_glu) .* ((1 - exp(-t / tau_rise_Ex)) .* exp(-t / tau_decay_Ex)) * (Vm - EGlu))', 'independent', 't', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [-70 0 1 1 -50 0 0 0 0];
opts.StartPoint = [-70 0 50 50 -50 1 1 1 1]; % Starting values
opts.Upper = [-70 0 100 100 -50 20 20 5 5];
%% Fit model to data
switch i
case 1
% Realistic values
if z == 1
[fitresult, gof1(z)] = fit(xData, yData, ft, opts)
% In order to not generate 5 similar plots, we
% generate only one, the first time
%% Plot fit with data
figure( 'Name', 'Fit' );
h = plot( fitresult, xData, yData );
legend( h, 'CPSC at -50mV', 'Fit to CPSC', 'Location', 'NorthEast', 'Interpreter', 'none');
subtitle('Realistic values')
% Label axes
xlabel( 'time', 'Interpreter', 'none' );
ylabel( 'pA', 'Interpreter', 'none' );
grid on
clear y
else
[fitresult, gof1(z)] = fit(xData, yData, ft, opts);
end
case 2
% More challenging situation: max values for conductances
% closer to each other
if z == 1
[fitresult, gof2(z)] = fit(xData, yData, ft, opts)
% In order to not generate 5 similar plots
%% Plot fit with data
figure( 'Name', 'Fit' );
h = plot( fitresult, xData, yData );
legend( h, 'CPSC at -50mV', 'Fit to CPSC', 'Location', 'NorthEast', 'Interpreter', 'none');
subtitle('Maximal values of conductances closer to each other')
% Label axes
xlabel( 'time', 'Interpreter', 'none' );
ylabel( 'pA', 'Interpreter', 'none' );
grid on
clear y
else
[fitresult, gof2(z)] = fit(xData, yData, ft, opts);
end
case 3
% More challenging situation: max values for conductances
% closer to each other and tau values closer to each other
if z == 1
[fitresult, gof3(z)] = fit(xData, yData, ft, opts)
% In order to not generate 5 similar plots
%% Plot fit with data
figure( 'Name', 'Fit' );
h = plot( fitresult, xData, yData );
legend( h, 'CPSC at -50mV', 'Fit to CPSC', 'Location', 'NorthEast', 'Interpreter', 'none');
subtitle('Maximal values of conductances and taus closer to each other')
% Label axes
xlabel( 'time', 'Interpreter', 'none' );
ylabel( 'pA', 'Interpreter', 'none' );
grid on
clear y
else
[fitresult, gof3(z)] = fit(xData, yData, ft, opts);
end
end
end
end
The problem is that the accuracy seems to depend a lot on the starting values I give, which is not desirable. In particular, if these are the starting values:
opts.StartPoint = [-70 0 50 50 -50 1 1 1 1]; % Starting values
The procedure does well. It almost predicts 80 and 15 for the first case, and 60 and 40 for the second and third.
However, if those two 50 become 1 and 1, the procedure doesn't do well anymore. It predicts 69 and 1 (instead of 80 and 15); 60 and 56 (instead of 60 and 40); 28 and 13 (instead of 60 and 40). Weird thing: if I then use these (wrong) predictions as starting values for another iterations, it works perfectly.
What am I doing wrong? Is there a way to improve it?
Thanks!

Alan Weiss on 12 Mar 2021
You might find some relevant information here or here. The point is that fitting problems typically have multiple local optima, and you often need to start from multiple values (meaning run multiple fits) to get the best (global) fit.
Alan Weiss
MATLAB mathematical toolbox documentation
Samuele Bolotta on 12 Mar 2021
Great, thanks!