Fitting procedure using MultiStart - doesn't recognize objective function

Question

0 个投票

I'm trying to fit a function to some data using MultiStart in order to find more reliably the global minimum.

This is the code:

[~,~,CPSC,t] = generate_current(80,15,0,-70,-30,0.44,15,0.73,3,120); % Generate current
% Initial values
gmc = 50;
gmg = 50;
tde = 1;
tdi = 1;
tre = 1;
tri = 1;
G_max_chl = 80;
G_max_glu = 15;
EGlu = -70;
EChl = 0;
Vm = -30;
tau_rise_In = 0.44;
tau_decay_In = 15;
tau_rise_Ex = 0.73;
tau_decay_Ex = 3;
y = awgn(CPSC,25,'measured'); % Add white noise
%% Perform fit
[xdata, ydata] = prepareCurveData(t, y);
lb = [-70 0 1   1   -30 0 0 0 0]; % Lower bound
ub = [-70 0 150 150 -30 20 20 5 5]; % Upper bound
p0 = [-70 0 gmc gmg -30 tde tdi tre tri]; % Starting values
% Set up fittype and options.
fitfcn = @(varargin) ((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));
[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,xdata,ydata,lb,ub);
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,50)
%% Plot fit with data
figure( 'Name', 'Fit' );
h = plot( fitresult1, xdata, ydata );
legend( h, 'CPSC at -30mV', 'Fit to CPSC', 'Location', 'NorthEast', 'Interpreter', 'none');
subtitle('Realistic values')
% Label axes
xlabel( 'time', 'Interpreter', 'none' );
ylabel( 'pA', 'Interpreter', 'none' );
grid on

And this is the function that I call at the beginning:

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

The error that I get is:

Error using lsqcurvefit (line 269)
Function value and YDATA sizes are not equal.
Error in global_m_fit (line 32)
[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,xdata,ydata,lb,ub);

I'm not sure why that is the case.

Thanks!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Star Strider 2021-3-15

See also Nonlinear fitting - No output function with MultiStart

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question

Answer 1

Walter Roberson 2021-3-15

在 MATLAB Online 中打开

1 个投票

fitfcn = @(varargin) ((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));

That ignores the input trial parameters, returning a constant expression. However the curve fitting routines pass in trial coefficients of a smaller size in order to validate the function, and your routine is failing to return something of the same size as the trial coefficients.

It is not clear to me what the point is in using a fitting function that ignores its inputs.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Samuele Bolotta 2021-3-15

在 MATLAB Online 中打开

Thanks!

But still, with this code:

%% Perform fit
[xData, yData] = prepareCurveData(t, y);
lb = [-70 0 1   1   -30 0 0 0 0]; % Lower bound
ub = [-70 0 150 150 -30 20 20 5 5]; % Upper bound
p0 = [-70 0 gmc gmg -30 tde tdi tre tri]; % Starting values
% Set up fittype and options
fitfcn = @(EGlu, EChl, G_max_chl, G_max_glu, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,t) ((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));
[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,t,y,lb,ub);
problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);
ms = MultiStart('PlotFcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,50)

I get this error:

Error in
global_m_fit>@(EGlu,EChl,G_max_chl,G_max_glu,Vm,tau_rise_In,tau_decay_In,tau_rise_Ex,tau_decay_Ex,t)((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))
(line 22)
fitfcn = @(EGlu, EChl, G_max_chl, G_max_glu, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,t) ((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));
Error in lsqcurvefit (line 225)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in global_m_fit (line 24)
[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,t,y,lb,ub);
Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

Walter Roberson 2021-3-21

在 MATLAB Online 中打开

fitfcn = @(EGlu, EChl, G_max_chl, G_max_glu, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,t) ((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));

The function you pass to lsqcurvefit must accept exactly two parameters, the first of which is the trial coefficients, and the second of which is the x locations to calculate the function at given those trial coefficients.

Your p0 has 9 entries, and your fitfcn expects 10 parameters, the last of which is t, which is probably the vector of independent locations.

My guess is that you want

fitfcn_aux = @(EGlu, EChl, G_max_chl, G_max_glu, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,t) ((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));
fitfcn = @(p,t) fitfcn_aux(p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8),p(9),t);

You have some other bugs too; I will post the updated code below.

Walter Roberson 2021-3-21

编辑：Walter Roberson 2021-3-21

在 MATLAB Online 中打开

In the below, the isunix() & ~ismac() is to attempt to detect whether the code is being run in the batch processing facility online -- a facility that displays the answers here in this post, but which does not happen to support the uicontrol('style','push') that gsplotbestf uses to give the user an opportunity to stop the plotting. If you just happen to be using linux yourself then you will want to recode to take the "else" branch

format long g

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

% Initial values

gmc = 50;

gmg = 50;

tde = 1;

tdi = 1;

tre = 1;

tri = 1;

G_max_chl = 80;

G_max_glu = 15;

EGlu = -70;

EChl = 0;

Vm = -30;

tau_rise_In = 0.44;

tau_decay_In = 15;

tau_rise_Ex = 0.73;

tau_decay_Ex = 3;

y = awgn(CPSC,25,'measured'); % Add white noise

%% Perform fit

[xData, yData] = prepareCurveData(t, y);

lb = [-70 0 1 1 -30 0 0 0 0]; % Lower bound

ub = [-70 0 150 150 -30 20 20 5 5]; % Upper bound

p0 = [-70 0 gmc gmg -30 tde tdi tre tri]; % Starting values

% Set up fittype and options

fitfcn_aux = @(EGlu, EChl, G_max_chl, G_max_glu, Vm, tau_rise_In, tau_decay_In, tau_rise_Ex, tau_decay_Ex,t) ((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));

fitfcn = @(p,t) fitfcn_aux(p(1),p(2),p(3),p(4),p(5),p(6),p(7),p(8),p(9),t);

%[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,t,y,lb,ub);

problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...

'lb',lb,'ub',ub,'xdata',xData,'ydata',yData);

if isunix() && ~ismac()

ms = MultiStart('PlotFcns', []);

else

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

end

[xmulti,errormulti] = run(ms,problem,50);

MultiStart completed some of the runs from the start points. 37 out of 50 local solver runs converged with a positive local solver exit flag.

disp(xmulti)

Columns 1 through 7 -70 0 149.999999999992 149.999997227555 -30 0.36130819595096 1.55714549746472 Columns 8 through 9 0.311540436223902 4.99999999999998

disp(errormulti)

138929197.825979

%% Plot fit with data

figure( 'Name', 'Fit' );

fitresult1 = fitfcn(xmulti,xData);

h = plot( xData, yData, xData, fitresult1 );

legend( h, 'CPSC at -30mV', 'Fit to CPSC', 'Location', 'NorthEast', 'Interpreter', 'none');

subtitle('Realistic values')

% Label axes

xlabel( 'time', 'Interpreter', 'none' );

ylabel( 'pA', 'Interpreter', 'none' );

grid on

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

Samuele Bolotta 2021-3-22

编辑：Samuele Bolotta 2021-3-22

在 MATLAB Online 中打开

Hi Walter, thank you for your great answer. I've re-written the code and it works now. Also, I've specified a bit better the parameters (three of them were constant, so there was no point in including them). I was surprised, however, that multistart is still performing the same as the normal fit function. You can see that by running 50 iterations of both.

True values = [80 15 3 8 1.5 0.2]

Both fit function and multistart (avg) = [85 29 3.42 8 2.5 0.2]

Same std for both, also.

Why do you think this is the case?

[EPSC,IPSC,CPSC,t] = generate_current(80,15,0,-70,-30,0.2,8,1.5,3,120);
% Create table for later
True_values = {80,15,3,8,1.5,0.2};
T = cell2table(True_values,'VariableNames',{'Fun_Max conductance Chloride','Fun_Max conductance Glutamate', 'Fun_Tau Decay Ex', 'Fun_Tau decay In', 'Fun_Tau rise Ex', 'Fun_Tau rise In'});
% Initial values for fitting
gmc = 40; gmg = 20; tde = 5; tdi = 6; tre = 1; tri = 1;
ub      = [150 5 20 -30 -70 150 5 20 0];
p0 = [gmc tri tdi -30 -70 gmg tre tde 0]; % Starting values
lb      = [1 0 0 -30 -70 1 0 0 0];
zmax = 30; % number of iterations
start_n = 50; % number of runs
xmulti = zeros(zmax,numel(p0)); % Preallocation
%% Fit model to data
% Objective function
fitfcn = @(p, xdata) ((p(1)) .* ((1 - exp(-xdata / p(1,2))) .* exp(-xdata / p(1,3))) * 40) + ...
    ((p(1,6)) .* ((1 - exp(-xdata / p(1,7))) .* exp(-xdata / p(1,8))) * -30);
for z = 1:zmax
    
    % 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');
        
    % Find the best local fit
    [x,resnorm,~,exitflag,output] = lsqcurvefit(fitfcn, p0, t, y,lb,ub);
        
    %Set up the problem for MultiStart.
    problem = createOptimProblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
            'lb',lb,'ub',ub,'xdata',t,'ydata',y);
        
    % Find a global solution
    ms = MultiStart('PlotFcn',@gsplotbestf);
    [xmulti(z,:),errormulti(z,:)] = run(ms,problem,start_n);
    end
end
% Create one table for each case. Store values for all 20 iterations
T1_ms = table(xmulti(:,1), xmulti(:,6), xmulti(:,8), xmulti(:,3), xmulti(:,7), xmulti(:,2));
T1_ms.Properties.VariableNames = {'Max conductance Chloride','Max conductance Glutamate', 'Tau Decay Ex', 'Tau decay In', 'Tau rise Ex', 'Tau rise In'};
func1 = @(x) mean(x);
func2 = @(x) std(x);
func3 = @(x) var(x);
func4 = @(x) range(x);
func5 = @(x) mad(x);
table_ms1(1,:) = varfun(func1,T1_ms);
table_ms1(2,:) = varfun(func2,T1_ms);
table_ms1(3,:) = varfun(func3,T1_ms);
table_ms1(4,:) = varfun(func4,T1_ms);
table_ms1(5,:) = varfun(func5,T1_ms);
table_ms1 = [T;table_ms1];
table_ms1.Properties.RowNames = {'True value','Mean','Std','Var','Range','Mad'};

Walter Roberson 2021-3-22

在 MATLAB Online 中打开

Look at your code

for z = 1:zmax
    
    % 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');

You assign to y inside the loop. Do you use the changed value of y inside the loop? Yes, you use it on the call to lsqcurvefit() and to createOptimProblem. So the assignment to y is "used".

    % Find the best local fit
    [x,resnorm,~,exitflag,output] = lsqcurvefit(fitfcn, p0, t, y,lb,ub);

You assign to x, resnorm, exitflag, and output inside the loop. Do you use any of those variables for anything else inside the loop? NO. So the very next iteration of z, whatever you had assigned to them will be discarded and replaced with the result for the next value of z. As far as the final results stored in those variables is concerned, you might as well not have done this call to lsqcurvefit() at all, and just do a single call after the for z loop is terminated, if you want to look at the values for debugging purposes. But you have no statements that use even one of those variables after the loop, so why bother making the call to lsqcurvefit() at all?

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

You assign to problem inside the loop? Do you use the changed value of problem inside the loop? Yes, you use it in the call to run()

    % Find a global solution
    ms = MultiStart('PlotFcn',@gsplotbestf);

You assign to ms inside the loop. Do you use the changed value of ms inside the loop? Yes, you use it in the call to run()

[xmulti(z,:),errormulti(z,:)] = run(ms,problem,start_n);

You assign to xmulti(z,:) and errormulti(z,:) . Do you use those inside the loop? No. You do, however, use the stored value for xmulti after the loop, so it is worth saving. You do not currently use errormulti after the loop, so it is not clear it is worth saving.

end

That ends the for z

end

That appears to be a stray end that is going to interfere with the structure of the code.

If it is indicating the end of some other loop, then you potentially have a problem as the output location xmulti(z,:) does not depend upon another other variable and so would be overwritten each iteration of the hypothetical outer loop that this end statement might maybe be there to mark the final point of.

Samuele Bolotta 2021-3-22

Ah yes, I was being naive. I fixed those things you're mentioning. I think the most interesting thing is that if we compare the mean of x (local solver) and the mean of xmulti (multistart, with 50 local solver runs) for 30 iterations, they are almost exactly the same. I'm sure the difference is not statistically significant.

So the problem seems to be related not to the starting points, but to the amount of noise (because with lower levels it works perfectly). I'll try Splitting the Linear and Nonlinear Problems (https://it.mathworks.com/help/optim/ug/nonlinear-data-fitting-problem-based-example.html).

Thanks again for your time!

请先登录，再进行评论。

Answer 2

Alan Weiss 2021-3-15

在 MATLAB Online 中打开

0 个投票

I think that you need to have just one input variable, typically called x, and have each of your other named variables be a component of x. For example,

function y = objfun(x)
EGlu = x(1);
EChl = x(2);
G_max_chl = x(3);
G_max_glu = x(4);
Vm = x(5);
tau_rise_In = x(6);
tau_decay_In = x(7);
tau_rise_Ex = x(8);
tau_decay_Ex = x(9);
t = x(10);
y = ((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));
end

For details, see Writing Scalar Objective Functions. But maybe I got it wrong, and you are using lsqnonlin. Make the appropriate modifications as in Writing Vector and Matrix Objective Functions.

Alan Weiss

MATLAB mathematical toolbox documentation

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Fitting procedure using MultiStart - doesn't recognize objective function

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

类别

产品

标签

Community Treasure Hunt

Fitting procedure using MultiStart - doesn't recognize objective function

1 个评论 显示 -1更早的评论 隐藏 -1更早的评论

采纳的回答

10 个评论 显示 8更早的评论 隐藏 8更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

类别

产品

标签

另请参阅

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论