Nonlinear fitting depends too heavily on starting points

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!

Answer 1

[Alan Weiss]

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]

Great, thanks!