Asked by Kelly McGuire
on 16 Jul 2019

My data requires two model functions for the curve fitting. I use one function for the first part of the data, and the second function for the last portion. I have the two model functions in different for loops. The second function picks up in the correct row vector element where the first function leaves off. Both model functions use PredCurrFun as the function name. It appears that the output from the second model function is overwriting the outputs from the first model function. How can I use two model functions in the same LSQCURVEFIT routing? Do I need to call one PredCurrFun and the other PredCurrFun2, and LSQCURVEFIT can use both? I tried setting n = 0 before the first for loop, and then use PredCurrFun(n) hoping that would write the output values in the correct order and position, but I get "Conversion to double from function_handle is not possible". In the code below, that's why there is n = 0 and n = n +1 in the for loops. How do I keep the second function from overwriting the first?

%Initialize Variables

%A = ;

Co = concentration; %Initial concentration in uM, make concentration into row vector i.e. 10 10 10 100 100 100 500 500 500 etc., in same concentration order as current row vector

Diff = 2; %Diffusion coefficient in um^2/ms

T = traceLength; %Trace length row vector

xdata = time;

ydata = current; %Current traces concatenated end to end, beginning with washin traces (lowest concentration to highest), and then washout traces added to end of washin traces (lowest conc. to highest)

%t = time2; %For plotting oberseved data, put time traces in their own column and import as numeric matrix

%y = current2; %For plotting observed data, put current traces in their own column and import as numeric matrix

%Initialize PredCurrFun vector with zeros

PredCurrFun = zeros(1,1271810);

%Plot of actual data

%plot(t1,y,'--b');

%xlabel('Time');

%ylabel('Current');

%Model function used by LSQCURVEFIT - %p(1) is in uM^(-1) ms^(-1), p(2) is

%in ms^(-1), p(3) is in um, total current p(4) is in nA, and p(5) is a

%unitless fraction between 0 and 1.

%Trace Length vector here

[k,m] = size(Co);

%Fitting Wash-in

n = 1;

for i = 1:m

t = 0;

fprintf('t: %d\n', t);

for j = T(i)+1:T(i+1) %Ex. T(1)+1 = 1 to T(1+1) = T(2) = 386610, T(2)+1 = 386611 to T(2+1) = T(3) = 556681, etc...

PredCurrFun = @(p,t) p(5).*p(4) + (1-p(5)).*p(4) .* ((1-p(2)./(p(1).*(Co-(4.*Co)./pi).*exp(-p(1).*(Co-(4.*Co./pi).*exp(-(Diff.*pi^2.*t./(4.*p(2)^2))) + p(2)).*t) ...

+ (p(2) ./ (p(1).*(Co-(4.*Co./pi).*exp(-Diff.*pi^2.*t./(4.*p(3)^2)) + p(2)))))));

t = t + 1;

n = n + 1;

end

end

%Fitting Wash-out

for k = 4:6

t = 0;

fprintf('t: %d\n', t);

for m = T(k)+1:T(k+1) %Ex. T(4)+1 = 945665+1 = 945666 to T(1+4) = T(5) = 973519

PredCurrFun = @(p,t) p(4).*(p(5)+(1-p(5)).*(1-(0.9997).*exp(-p(2).*t)));

t = t + 1;

n = n + 1;

end

end

%Parameter upper, lower, and starting values

lb = [0,0,0,-250,0]; %k1 is in uM^(-1) ms^(-1), k2 is in ms^(-1), L is in um, and total current, I, is in nA

ub = [1E-5,1E-5,300,-50,1];

startingVals = [3E-7,3E-7,10,-140,0.5];

%Non-linear least squares fit of data using model function

format long

options=optimset('disp','iter','LargeScale','off','TolFun',0.000001,'MaxIter',100000,'MaxFunEvals',1000000);

[p,resnorm] = lsqcurvefit(PredCurrFun, startingVals, xdata, ydata, lb, ub, options)

%Chisquare Calculation

ChiSquareRed = resnorm ./ size(ydata)

%Plot LSQCURVEFIT Result

%The number in parentheses next to Co is the row number of the

%concentration column number. For example, Co(1) refers to the first row, which in

%this case was 10 uM, Co(4) refers to the fourth row which was 100 uM, and

%so on.

xgrid = linspace(0,1271810,1271810);

line(xgrid, PredCurrFun(p,xgrid), 'Color', 'r');

Answer by Star Strider
on 16 Jul 2019

‘Do I need to call one PredCurrFun and the other PredCurrFun2, and LSQCURVEFIT can use both?’

Yes, in separate calls and other calculations.

If I understand your Question correctly (no promises), they are two different functions applicable (apparently) to two different segments of your ‘time’ vector. I would use two separate lsqcurvefit calls, one for each segment of your study. Then calculate the fitted functions and plot them with respect to the data they are fitted to. You can then plot them serially, in the same plot call.

Star Strider
on 16 Jul 2019

The lsqcurvefit function would try to estimate ‘p(1)’ and ‘p(2)’, unless you changed the second funciton to have them defined already (so they would not change) and only have the others in the parameter vector.

Something like this:

Do the first lsqcurvefit call, then:

p1prev = p(1);

p2prev = p(2);

then define in the second function:

p1 = p1prev;

p2 = p2prev;

p4 = p(1); % ‘p(4)’ Is The First Parameter In The New Parameter Vector

p5 = p(2); % ... Similarly

PredCurrFun = @(p,t) p4.*(p5+(1-p5).*(1-(0.9997).*exp(-p2.*t)));

That preserves the first parameter estimates, then allows lsqcurvefit to estimate the others. (Alternatively, in the second lsqcurvefit call, you could simply constrain ‘p(1)’ and ‘p(2)’ to be exactly those values, however I would not advise that.)

Kelly McGuire
on 16 Jul 2019

Great, I will give that a try, thanks.

Star Strider
on 16 Jul 2019

My pleasure.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.