Make lsqcurvefit go through a defined point?

I am trying to run lsqcurvefit to model some chemical titration data. It uses a custom function I define. However, the fitted result returns a negative x-value at y(0), when I need it to pass through (0,y(0)). I'm attaching an example of my input data. Here is my code:

fun=@(out,xdata) (xdata(:,2).*(out(1)./(1+(xdata(:,1)./out(2))))) + (xdata(:,2).*(xdata(:,5)./(1+(xdata(:,1)./xdata(:,6))))) - (xdata(:,2) + xdata(:,3)).*(xdata(:,1)-xdata(:,4)) + xdata(:,2).*out(3);
out0=[10^-5 10^-4 xdata(1,1)];
options = optimoptions(@lsqcurvefit,'StepTolerance',10^-50,'OptimalityTolerance',10^-15,'FunctionTolerance',10^-15,'MaxFunctionEvaluations',10000);
out_OneGroup = lsqcurvefit(fun,out0,xdata,ydata,[1e-6 1e-5 0],[1e-3 1e-3 0.1],options);
OneGroup_results=feval(fun,out_OneGroup,xdata);

I can take the first value of OneGroup_results (-3.3866e-6), subtract it from the 3rd out_OneGroup value (increasing that value by the offset 3.3866e-6), and re-run the program and all values will be above zero, but the resulting fit it also offset and doesn't work.

Any suggestions greatly appreciated!

[Christopher Hunt]

As a graphical example, here are plots of the fits without (red) and with (cyan) the offset to bring the fit through (x=0,y=y(0))

[Matt J]

The curve model you've shown is not a 1D curve. It is a mapping from . Moreover, xdata(:,7) in your example.mat appears to be unused. Therefore, it is not clear how the plots you've shown are constructed.

[Christopher Hunt]

Sorry, my code is kind of messy. Here is an updated example.mat file, xdata(:,7) was indeed unused. And below is the complete code I am using to generate the fit and plot. I'm not sure what mapping R6-->R means though.

fun = @(out)((xdata(:,2).*(out(1)./(1+(xdata(:,1)./out(2))))) + (xdata(:,2).*(xdata(:,5)./(1+(xdata(:,1)./xdata(:,6))))) - (xdata(:,2) + xdata(:,3)).*(xdata(:,1)-xdata(:,4)) + xdata(:,2).*out(3))-...
            ((xdata(1,2).*(out(1)./(1+(xdata(1,1)./out(2))))) + (xdata(1,2).*(xdata(1,5)./(1+(xdata(1,1)./xdata(1,6))))) - (xdata(1,2) + xdata(1,3)).*(xdata(1,1)-xdata(1,4)) + xdata(1,2).*out(3));
fun1 = @(out) fun(out) - ydata;
out = [10^-5 10^-4 xdata(1,1)];
options = optimoptions(@lsqnonlin,'StepTolerance',10^-50,'OptimalityTolerance',10^-15,'FunctionTolerance',10^-15,'MaxFunctionEvaluations',10000);
out_OneGroup = lsqnonlin(fun1,out,[1e-6 1e-5 0],[1e-3 1e-3 0.1],options);
OneGroup_results=fun(out_OneGroup);

figure
plot(ydata,-1.*log10(xdata(:,1)),'k')
hold on
plot(OneGroup_results,-1.*log10(xdata(:,1)),'r.')
hold on
pK1=-1*log10(out_OneGroup(2));
pK1_diff=abs(-1.*log10(xdata(:,1))-pK1);
idx=find(pK1_diff==(min(pK1_diff)));
hold on
scatter(ydata(idx),pK1,200,'b*')
xlabel('mol NaOH added')
ylabel('pH')
title('One Group')
legend('Obs','Modeled','pKa1')

[Matt J]

It uses a custom function I define. However, the fitted result returns a negative x-value at y(0), when I need it to pass through (0,y(0)).

But in your code above, you are plotting ydata on the x-axis, so it is still not clear in your statement what you consider the "x-values" and whether "y(0)" means ydata(1) or something else.

[Christopher Hunt]

OK, let me start from scratch and try to use terms consistently. I want to fit the data in ydata using the custom function and the input data in xdata. When I fit using the function fun0, the modeled y(0) value is less than the initial value of ydata(0), i.e. it is not 0. However, the general fit of fun0 is OK, excepting the y(0) issue.

But when I use fun1 to force the fit through y(0)==0, the resulting shape of the fit is terrible (i.e. offset from the observations).

Example data attached here, and plot and full code excerpt below.

clear all
load example.mat

out = [10^-5 10^-4 xdata(1,1)]; %initial coefficient guesses

fun0=@(out,xdata) (xdata(:,2).*(out(1)./(1+(xdata(:,1)./out(2))))) + (xdata(:,2).*(xdata(:,5)./(1+(xdata(:,1)./xdata(:,6))))) - (xdata(:,2) + xdata(:,3)).*(xdata(:,1)-xdata(:,4)) + xdata(:,2).*out(3);
options = optimoptions(@lsqcurvefit,'StepTolerance',10^-50,'OptimalityTolerance',10^-15,'FunctionTolerance',10^-15,'MaxFunctionEvaluations',10000);
out_fun0 = lsqcurvefit(fun0,out,xdata,ydata,[1e-6 1e-5 0],[1e-3 1e-4 0.1],options);
y_fun0 =feval(fun0,out_fun0,xdata);

fun1 = @(out,xdata)((xdata(:,2).*(out(1)./(1+(xdata(:,1)./out(2))))) + (xdata(:,2).*(xdata(:,5)./(1+(xdata(:,1)./xdata(:,6))))) - (xdata(:,2) + xdata(:,3)).*(xdata(:,1)-xdata(:,4)) + xdata(:,2).*out(3))-...            
    ((xdata(1,2).*(out(1)./(1+(xdata(1,1)./out(2))))) + (xdata(1,2).*(xdata(1,5)./(1+(xdata(1,1)./xdata(1,6))))) - (xdata(1,2) + xdata(1,3)).*(xdata(1,1)-xdata(1,4)) + xdata(1,2).*out(3));
out_fun1 = lsqcurvefit(fun1,out,xdata,ydata,[1e-6 1e-5 0],[1e-3 1e-4 0.1],options);
y_fun1 =feval(fun1,out_fun1,xdata);

figure
plot(-1.*log10(xdata(:,1)),ydata,'k')
hold on
plot(-1.*log10(xdata(:,1)),y_fun0,'r.')
hold on
plot(-1.*log10(xdata(:,1)),y_fun1,'m.')
ylabel('mol NaOH added')
xlabel('pH')

legend('Observed','fun0','fun1','fun2')

[Torsten]

As said, your model function fun0 is not appropriate for fitting if it does not automatically satisfy fun0(0) = 0.

[Christopher Hunt]

But isn't that the point of multivariate least squares fitting: that a solution can be found using a combination of input data and adjusted model coefficients which satisfies the given function but may not pass through some or most of the points, including the first one? My essential question is whether I can force the routine to make the solution go through one of the given points such that fun0(xdata(1,:))=0. It seems like the answer is yes in a sense is I use fun1(), but then all other fun1() points have an offset. So I am in search of a way to do both: force the solution through fun0(xdata(1,:)) and model the other data points without the offset.

I will also throw out there that the equation underpinning fun0() comes from a couple publications which have used it in the same way to model titration data identical to those I am working with. One of these papers notes a potential need to adjust out(3) to make it result in fun0(out,xdata(1,:))==0 but doesn't say how they adjust it. I think fun1() is trying to do that, but something is not quite right.

[Torsten]

But isn't that the point of multivariate least squares fitting: that a solution can be found using a combination of input data and adjusted model coefficients which satisfies the given function but may not pass through some or most of the points, including the first one?

Of course - including the first one. But you want to exclude the first one because you want the function to pass through this point.

Say you have a model function f0(t) = exp(k*t) for the concentration of a substance over time and you know you start with a concentration of 0 - would you take this model function which starts with a concentration of 1 ? Of course, f1(t) = exp(k*t) - exp(k*0) = exp(k*t) - 1 could be used to satisfy f1(0) = 0. But it is only one possibility to force f0 to pass through (0,0), so you don't know if it is also suited for t>0.

Answer 1

[Torsten]

fun = @(out)((xdata(:,2).*(out(1)./(1+(xdata(:,1)./out(2))))) + (xdata(:,2).*(xdata(:,5)./(1+(xdata(:,1)./xdata(:,6))))) - (xdata(:,2) + xdata(:,3)).*(xdata(:,1)-xdata(:,4)) + xdata(:,2).*out(3))-...
            ((xdata(1,2).*(out(1)./(1+(xdata(1,1)./out(2))))) + (xdata(1,2).*(xdata(1,5)./(1+(xdata(1,1)./xdata(1,6))))) - (xdata(1,2) + xdata(1,3)).*(xdata(1,1)-xdata(1,4)) + xdata(1,2).*out(3));
fun1 = @(out) fun(out) - ydata;
out0 = [10^-5 10^-4 xdata(1,1)];
out_OneGroup = lsqnonlin(fun1,out0);
plot(1:numel(xdata(:,1)),fun(out_OneGroup))
hold on
plot(1:numel(xdata(:,1)),ydata)
hold off

[Christopher Hunt]

Many thanks, that makes sense. The resulting fit does go through (x,y(0)), but the fit is pretty badly offset from the observations. I'm including a plot below. I might be thinking about this wrong, but it seems like it should be possible to drive the fit through (x,y(0)) and have the fit follow the observations closer than this?

[Torsten]

I might be thinking about this wrong, but it seems like it should be possible to drive the fit through (x,y(0)) and have the fit follow the observations closer than this?

Then you are more clever than "lsqnonlin" ...

But seriously: It's kind of strange that your function does not automatically fulfill the constraint that it goes through (x,y(0)). Usually, fitting objects are chosen such that they follow the physics they reflect.

You may try different starting values and see whether the approximation gets better. Look up "MultiStart".

[Christopher Hunt]

Thanks, I gave MultiStart a try, but it gavethe same result coefficients as lsqnonlin, so it seems like a local minimum doesn't seem to be the issue