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

2 views (last 30 days)
Yi Jiao on 6 Oct 2020
Edited: Matt J on 7 Oct 2020
[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) exp (-k[t]); and Ymax, Y 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!
Yi Jiao on 6 Oct 2020
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) exp (-k[t]); and Y[max], Y 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!

Matt J on 6 Oct 2020
Edited: Matt J on 6 Oct 2020
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
Matt J on 7 Oct 2020
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...

J. Alex Lee on 7 Oct 2020
Edited: J. Alex Lee on 7 Oct 2020
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
J. Alex Lee on 7 Oct 2020
very cool, i was ignorant of the phrase "partitioned least squares", thanks!