Fitting procedure using MultiStart - doesn't recognize objective function

3 views (last 30 days)
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!

Accepted Answer

Walter Roberson
Walter Roberson on 15 Mar 2021
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 Comments
Walter Roberson
Walter Roberson on 22 Mar 2021
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
Samuele Bolotta on 22 Mar 2021
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!

Sign in to comment.

More Answers (1)

Alan Weiss
Alan Weiss on 15 Mar 2021
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!