How to solve unknown coefficients for an exponential fit equation iteratively?

[updated question]

I measured concentrtions [Y] of a gas of interest in an enclosed enviroment over time [t] and I got three data points. They are:

t = [1 7 15]; % min

Y = [557 792 799]; % parts per tillion

The concentration of the gas is suppose to reach an equailibrium level given long eough time, which is Ymax. (in this case, Ymax should be around 800).

And I would like to fit an empirical euqation to it, which is Y(t) = Ymax - (Ymax-Y[0]) exp (-k[t]); and Ymax, Y[0] and k are all unknown. I have seen some papers saying that they solved these three coeffs iteratively, but I just don't know how.

Thank you!

___________________

[Previous Q]

Supoose I have a data set of

x: 0; 7; 15

y: 557; 792; 799;

And I would like to get an exponential fit to the dataset as y(x) = a - (a-y(0))*exp(-kx);

How should I do this iteratively? I just have no idea. Any thoughs would be helpful. Thank you!

[Ameer Hamza]

Why iteratively? You can use built-in functions, such as fit().

[Alan Stevens]

Or, if you don't have the toolbox with fit, look up fminsearch (which should work ok with initial guesses of y0 = 557, a = 800 and k = 1);

[J. Alex Lee]

What do you mean by y(0)...do you mean the value 557?

Say you stepped back and asked if you can fit

If you meant y(0) is unknown, then your problem is equivalent to this one, and you just have the 3 parameter fit with 3 data points. Then it's not a "fitting" problem anymore, but a root-finding one (you get an exact "fit"). You can solve it as a 3D nonlinear algebraic equation, or manipulate it down to a 1D nonlinear algebraic equation.

If you meant by y(0) that you want , then you will have to do a minimization (e.g., least squares problem), but as a result you likely won't be able to satisfy .

[Yi Jiao]

Thank you for your answers. Sorry for the confusion in my writing. Here is a detailed description of my question:

I measured concentrtions [Y] of a gas of interest in an enclosed enviroment over time [t] and I got three data points. They are:

t = [1 7 15]; % min

Y = [557 792 799]; % parts per tillion

The concentration of the gas is suppose to reach an equailibrium level given long eough time, which is Y[max]. (in this case, Ymax should be around 800).

And I would like to fit an empirical euqation to it, which is Y(x) = Y[max] - (Y[max]-Y[0]) exp (-k[t]); and Y[max], Y[0] and k are all unknown. I have seen some papers saying that they solved these three coeffs iteratively, but I just don't know how.

Thank you!

Answer 1

[Matt J]

Here's one way, using fminbnd,

x = [0; 7; 15];
y = [557; 792; 799];


fun=@(k)modelfun(k,x,y);
[k,resnorm]=fminbnd(fun,-10,10);

[~,a]=fun(k);

k,a,resnorm

k = 0.5036

a = 799.1275

resnorm = 9.1306e-04

function [resnorm,a]=modelfun(k,x,y)
 
 x=x(:); y=y(:); y0=y(1);

 e=exp(-k*x);

  a=(1-e)\(y-y0*e);
 
 ymodel=a-(a-y0).*e;
  
 resnorm=norm(y-ymodel);

end

[J. Alex Lee]

nice! You don't lose much in this case by using fminsearch with a starting point of k=0, since, it's not obvious why you'd know to pick the limits -10,10 for fminbnd.

[Walter Roberson]

Yi Jiao thanks Matt J.

[Matt J]

it's not obvious why you'd know to pick the limits -10,10 for fminbnd.

For the given x-values, any value of k outside of the interval [-10,10] would result in absurdly large inputs k*x to the exp() operator. Also, I vaguely wonder if fminbnd is more resilient to local minima than fminsearch...

Answer 2

[J. Alex Lee]

MattJ's answer takes y0 to be known, which Yi Jiao clarified is not the case. Treating y0 unknown and extending MattJ's answer leads to the following, which is basically to do a "2 step" fit finding the pair (y0,ymax) by a linear least squares at some given guess of k. I also implemented the "simpler" 3D search strategy suggested by Alan, and reuse the functions in the 1D strategy for evaluating the model and objective function.

The "iterative" strategy is hidden in the fminsearch algorithm. As others have noted, you can use fitting tools in the optimization toolbox, or write your own Gauss-Newton-like algorithms if you want to do it "correctly", for which you'd need the resfcn.

clc;
close all;
clear;

% data
x = [0; 7; 15];
y = [557; 792; 799];

% initial guess for the 3D fminsearch
p0 = [1 min(y) max(y)];

% execute 3D fminsearch
[p,fval] = fminsearch(@(p)objfcn(p,x,y),p0,optimset('Display','iter'))

% plot results
figure(1); cla; hold on;
plot(x,y,'o','LineWidth',2)
fplot(@(x)modelfcn(p,x),[0,16],'-','LineWidth',2)

% execute 1D fminsearch, similar to MattJ's answer
% but extending to y0 unknown
k0 = 1;
[k,fval] = fminsearch(@(k)objfcn_1(k,x,y),k0,optimset('Display','iter'))

% evalute y0 and ymax
[~,y0,ymax] = objfcn_1(k,x,y);

% plot result, which is the same
fplot(@(x)modelfcn([k;y0;ymax],x),[0,16],'--','LineWidth',2)


%% function for 1D fminsearch on k
%  using linear least squares to find y0,ymax based on k

function [obj,y0,ymax] = objfcn_1(k,x,y)

z = exp(-k*x);

q = [z,ones(size(z))]\y;

ymax = q(2);
y0   = sum(q);

obj = objfcn([k;y0;ymax],x,y);

end


%% functions for 3D fminsearch

function obj = objfcn(p,x,y)
obj = norm(resfcn(p,x,y));
end

function res = resfcn(p,x,y)
res = modelfcn(p,x)./y - 1;
% res = modelfcn(p,x) - y;
end

function y = modelfcn(p,x)

ymax = p(3);
y0   = p(2);
k    = p(1);

y = ymax - (ymax - y0) * exp(-k*x);

end

[Matt J]

Note that the 1D search is already implemented as a general File Exchange tool, fminspleas,

https://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

funlist={1,@(k,x)exp(-k*x)};

[k,lin]=fminspleas(funlist,0,x,y);
Ymax=lin(1);
Y0=Ymax+lin(2);

[J. Alex Lee]

very cool, i was ignorant of the phrase "partitioned least squares", thanks!