Nonlinear fitting - No output function with MultiStart

I'm fitting a function to some nonlinear data. However, the output that I get is:

  MultiStart with properties:

          UseParallel: 0
              Display: 'final'
    FunctionTolerance: 1.0000e-06
              MaxTime: Inf
            OutputFcn: []
              PlotFcn: []
     StartPointsToRun: 'all'
           XTolerance: 1.0000e-06

No plot and not output function.

This is the code:

% Parameters
p(1) = G_max_chl;
p(2) = G_max_glu;
p(3) = tau_rise_In;
p(4) = tau_decay_In;
p(5) = tau_rise_Ex;
p(6) = tau_decay_Ex;

% Create the objective function

fitfcn = @(p, Vm, EGlu, EChl, dt, tmax) ((p(1)) .* ((1 - exp(-(0:dt:tmax-dt) / p(3))) .* exp(-(0:dt:tmax-dt) / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-(0:dt:tmax-dt) / p(5))) .* exp(-(0:dt:tmax-dt) / p(6))) * (Vm - EGlu));

% Create the training data

rng default % For reproducibility
N = 200; % Number of data points
preal = [80,15,0.44,0.73,15,3]; % Real coefficients
Vm = -30;
Eglu = 0;
EChl = -70;
dt = 0.1;
tmax = 120;
ydata = fitfcn(preal,xdata); % Response data with noise
ydata = awgn(ydata,35,'measured');

% Set bounds and initial point.
lb = [0,0,0,0,0,0];
ub = [150,150,5,5,20,20];
p0 = 5*ones(1,6); % Arbitrary initial points

% Find the best local fit
[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,xdata,ydata,lb,ub)

%Set up the problem for MultiStart.
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);

% Find a global solution.
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,200)

I am not sure what I'm doing wrong.

Answer 1

[Star Strider]

Unfortunately, this:

fitfcn = @(p, Vm, EGlu, EChl, dt, tmax) ((p(1)) .* ((1 - exp(-(0:dt:tmax-dt) / p(3))) .* exp(-(0:dt:tmax-dt) / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-(0:dt:tmax-dt) / p(5))) .* exp(-(0:dt:tmax-dt) / p(6))) * (Vm - EGlu));

is not going to work with lsqcurvefit, or likely any other optimisation function.

The problems (that I can easily see) are:

fitfcn = @(p, Vm, EGlu, EChl, dt, tmax) ((p(1)) .* ((1 - exp(-(0:dt:tmax-dt) / p(3))) .* exp(-(0:dt:tmax-dt) / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-(0:dt:tmax-dt) / p(5))) .* exp(-(0:dt:tmax-dt) / p(6))) * (Vm - EGlu));
                                                                ↑ ← 't  MUST BE AN INDEPENDENT VARIABLE ↑                 ↑ ←  ↑ ← EXTRA PARAMETERS               ↑ ←  t                           ↑← 't                 ↑ ←  ↑ ← EXTRA PARAMETERS         

So it should be:

 fitfcn = @(p, t, Vm, EGlu, EChl) ((p(1)) .* ((1 - exp(-t / p(3))) .* exp(-t / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-t / p(5))) .* exp(-t / p(6))) * (Vm - EGlu));

and the lsqcurvefit call:

 B = lsqcurvefit(@(p,t)fitfcn(p, t, Vm, EGlu, EChl), B0, t, ydata)

with ‘Vm’, ‘EGlu’ and ‘EChl’ existing in the calling function workspace.

That should come close to working, although it may require some revision.

I cannot test it so I am posting it as UNTESTED CODE.

[Samuele Bolotta]

Thank you! I attach the whole code so that you can test it as well.

% Simulated data
[~,~,CPSC,t] = generate_current(80,15,0,-70,-30,0.44,15,0.73,3,120);

% % Parameters
G_max_chl = 80;                          % Maximal conductance of Chl
G_max_glu = 15;                          % Maximal conductance of Glu
tau_rise_Ex = 0.73;                      % Tau rise for E
tau_rise_In = 0.44;                      % Tau rise for I
tau_decay_Ex = 3;                        % Tau decay for E
tau_decay_In = 15;                       % Tau decay for I

p(1) = G_max_chl;
p(2) = G_max_glu;
p(3) = tau_rise_In;
p(4) = tau_decay_In;
p(5) = tau_rise_Ex;
p(6) = tau_decay_Ex;
dt = 0.1;
tmax = 120;
t = 1:dt:tmax;
Vm = -30;
Eglu = 0;
EChl = -70;

% Create the objective function
fitfcn = @(p, t, Vm, EGlu, EChl) ((p(1)) .* ((1 - exp(-t / p(3))) .* exp(-t / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-t / p(5))) .* exp(-t / p(6))) * (Vm - EGlu));

% Create the training data
rng default % For reproducibility
N = 200; % Number of data points
preal = [80,15,0.44,0.73,15,3]; % Real coefficients

xdata = CPSC;
ydata = fitfcn(preal,xdata); % Response data with noise

% Set bounds and initial point.
lb = [0,0,0,0,0,0];
ub = [150,150,5,5,20,20];
p0 = 5*ones(1,6); % Arbitrary initial points

% Find the best local fit
[xfitted,errorfitted] = lsqcurvefit(@(p,t)fitfcn(p, t, Vm, EGlu, EChl), B0, t, ydata)
% [xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,xdata,ydata,lb,ub)

%Set up the problem for MultiStart.
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);

% Find a global solution.
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,200)

Now I get this error:

Not enough input arguments.

Error in MultiStart>@(p,t,Vm,EGlu,EChl)((p(1)).*((1-exp(-t/p(3))).*exp(-t/p(4)))*(Vm-EChl))+((p(2)).*((1-exp(-t/p(5))).*exp(-t/p(6)))*(Vm-EGlu)) (line 19)
fitfcn = @(p, t, Vm, EGlu, EChl) ((p(1)) .* ((1 - exp(-t / p(3))) .* exp(-t / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-t / p(5))) .* exp(-t / p(6))) * (Vm - EGlu));

Error in MultiStart (line 31)
ydata = fitfcn(preal,xdata); % Response data with noise

I don't understand why, since I defined every input value I think.

[Star Strider]

As always, my pleasure!

Problem:

Unrecognized function or variable 'generate_current'.

That aside, It’s necessary to use ‘fitfcn’ will all the arguments:

ydata = fitfcn(preal,xdata, Vm, EGlu, EChl); % Response data with noise

[Samuele Bolotta]

Ah sorry, I didn't include the function.

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

So, I realized I was specifying the xdata and the ydata wrongly. In theory, xdata is time and ydata is CPSC. Both 1x1200. I want to fit fitfcn to this data. Now the code works better, but there's an error at the end:

Error in MultiStart (line 42)
ms = MultiStart('PlotFcns',@gsplotbestf);

This is the code:

% Simulated data
[~,~,CPSC,t] = generate_current(80,15,0,-70,-30,0.44,15,0.73,3,120);

% % Parameters
G_max_chl = 80;                          % Maximal conductance of Chl
G_max_glu = 15;                          % Maximal conductance of Glu
tau_rise_Ex = 0.73;                      % Tau rise for E
tau_rise_In = 0.44;                      % Tau rise for I
tau_decay_Ex = 3;                        % Tau decay for E
tau_decay_In = 15;                       % Tau decay for I

p = [G_max_chl,  G_max_glu, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex];

dt = 0.1;
tmax = 121;
t = 1:0.1:tmax-dt;
Vm = -30;
EGlu = 0;
EChl = -70;

% Create the objective function
fitfcn = @(p, t, Vm, EGlu, EChl) ((p(1)) .* ((1 - exp(-t / p(3))) .* exp(-t / p(4))) * (Vm - EChl)) + ((p(2)) .* ((1 - exp(-t / p(5))) .* exp(-t / p(6))) * (Vm - EGlu));

preal = [80,15,0.44,0.73,15,3]; % Real coefficients

xdata = t;
ydata = CPSC; % Response data with noise

% Set bounds and initial point.
lb = [0,0,0,0,0,0];
ub = [150,150,5,5,20,20];
p0 = [50,50,1,1,1,1]; % Arbitrary initial points

% Find the best local fit
[xfitted,errorfitted] = lsqcurvefit(@(p,t)fitfcn(p, t, Vm, EGlu, EChl), p0, xdata, ydata,lb,ub)

%Set up the problem for MultiStart.
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);

% Find a global solution.
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,200)

So it looks like I need better data, a better model or better starting values. What do you think?

Thanks!

[Star Strider]

The error appears to be here:

ms = MultiStart('PlotFcns',@gsplotbestf);

that will likely be solved by eliminating the ‘s’:

ms = MultiStart('PlotFcn',@gsplotbestf);

I usually use ga (genetic algorithm) of this sort of problem because I have more experience with it. I have used GlobalSearch as well a few times, however not MultiStart.

If the data are generated by your actual objective function to which you just added noise, the Golbal Optimization Toolbox functions will likely converge on the correct parameter set providing that the noise does not completely obscure the underlying signal. Providing the correct magnitudes of the initial parameter set is likely important (or likewise using the bound limits to force them into a specific range of magnitudes), however the initial parameter estimates provided to global optimisation functions are likely less important because the algorithms will eventually find the best set.

[Samuele Bolotta]

Thanks for your explanation!

Unfortunately, even removing the 's' it still gives me the same error. Doing this, I can find a local but not the global minimum. I'll try with GA, but would like to understand what I'm doing wrong here first.

% Find the best local fit
[xfitted,errorfitted] = lsqcurvefit(@(p,t)fitfcn(p, t, Vm, EGlu, EChl), p0, xdata, ydata,lb,ub)

%Set up the problem for MultiStart.
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);

% Find a global solution
ms = MultiStart('PlotFcn',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,200)