Force a starting point on exponential graph

2022 9 24
2 个回答
更新时间：2022 9 27
4 次查看（30 天）

0 个投票

Hi,
I'm trying to fit an exponential curve to a process describing a radioactive decay. This is my data:
counts=[2423970,2171372,2065862,1830553,1100899,1037972,914015,752138,684123,606126];
normalized_counts=counts./counts(1);
time=[24.4,40.1,49.2,69.9,137.1,144.4,160.6,185.4,192.7,209.7];
At first I tried using the loveable cftool, which easily managed to fit a "classic" exponential function: y=a*exp(b*x). However, I'm interested in fitting a "pure" exponential [normalized_counts=exp(-a*time)]. When I tried to insert this custom function, I got this error:
Inf computed by model function, fitting cannot continue.
Try using or tightening upper and lower bounds on coefficients.
Then I tried to force an initial condition, so that the exponent will have to pass through the point (24.4, 1) using the code below, but to no avail as well:
x0=24.4;
y0=1;
g = @(p,time)y0*exp(-p*(time-x0));
f = fit(time,normalized_counts,g)
plot(f,time,normalized_counts)
[I got some matrix size errors, tried to add ' (to transpose) to the arrays but that didn't help as well].
Would apprecaite your help!

请先登录，再回答此问题。

回答（2 个）

Davide Masiello
Davide Masiello 2022-9-24
编辑：Davide Masiello 2022-9-25
Ran in:

0 个投票

Ran in:
See if this can help
counts = [2423970,2171372,2065862,1830553,1100899,1037972,914015,752138,684123,606126];
normalized_counts = counts./counts(1);
time = [24.4,40.1,49.2,69.9,137.1,144.4,160.6,185.4,192.7,209.7];
p0 = 1e-3;
g = @(p,x) exp(-p*(x-time(1)));
f = fit(time',normalized_counts',g,'StartPoint',p0)
f =
General model: f(x) = exp(-p*(x-time(1))) Coefficients (with 95% confidence bounds): p = 0.007151 (0.006894, 0.007408)
plot(f,time,normalized_counts)

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Ran Kagan
Ran Kagan 2022-9-24
Thanks for the help!
It seems like having your way dramatically influences the b coefficient; while having the default fit [a*exp(b*x)] results in the coefficient being quite close to the theoretial value I expect, applying your method creates an error of ~20% from that value (I also tried to define p0 to be a smaller value but it seems to leave the coefficient constant).
Davide Masiello
Davide Masiello 2022-9-24
Ran in:
Ran in:
Yes, like below right?
counts = [2423970,2171372,2065862,1830553,1100899,1037972,914015,752138,684123,606126];
normalized_counts = counts./counts(1);
time = [24.4,40.1,49.2,69.9,137.1,144.4,160.6,185.4,192.7,209.7];
x0 = [1 1e-3];
g = @(a,b,x) a*exp(-b*x);
f = fit(time',normalized_counts',g,'StartPoint',x0)
f =
General model: f(x) = a*exp(-b*x) Coefficients (with 95% confidence bounds): a = 1.211 (1.179, 1.243) b = 0.007279 (0.006961, 0.007598)
plot(f,time,normalized_counts)
I didn't consider this because in your code you had specified y0 as a given value rather than an adjustable parameter.
Sam Chak
Sam Chak 2022-9-24
Yes, that is also what I interpreted "pure" exponential function, normalized_counts = exp(-a*time).
It's better for @Ran Kagan to clarify what is truly needed.
Ran Kagan
Ran Kagan 2022-9-25
编辑：Ran Kagan 2022-9-25
@Sam Chak I did explain in the original post that I'm looking for a fit in the form of y=exp(a*x) rather than the classic one provided by cftool: a*exp(b*x). I hope it was clear enough :)
@Davide Masiello What you showed here can be achieved by using the cftool, even without further coding, and yet this is not the exact form I want, although it gives a close value to the theoretial one.
I'm trying to normalize the factor multiplying the exponential form, this is why I'm trying to force the fit to pass on a pre-defined point, but seems like it comes with the cost of propagating a substantial error?
Davide Masiello
Davide Masiello 2022-9-25
@Ran Kagan I understand a bit better now, thanks. I have modified my answer accordingly, see if it helps.
Davide Masiello
Davide Masiello 2022-9-27
My pleasure, if that's the answer you were looking for don't forget to accept it!

请先登录，再进行评论。

Torsten
Torsten 2022-9-25
编辑：Torsten 2022-9-25
Ran in:

0 个投票

Ran in:
You must normalize to time = 0, not time = 24.4.
The initial mass is the key, not the mass when already 24.4 (whatever) have passed.
I think what you should try to do is to introduce another unknown "counts_initial" and fit the curve as
counts_normalized = exp(-p(1)*time)
with
time=[24.4,40.1,49.2,69.9,137.1,144.4,160.6,185.4,192.7,209.7];
counts_normalized = counts/counts_initial
This means that the curve you have to fit is
counts/counts_initial = exp(-p(1)*time)
thus again the general exponential
counts = counts_initial*exp(-p(1)*time)
counts = [2423970,2171372,2065862,1830553,1100899,1037972,914015,752138,684123,606126];
time = [24.4,40.1,49.2,69.9,137.1,144.4,160.6,185.4,192.7,209.7];
p0 = [3e6,1e-3];
g = @(p,x) p(1)*exp(-p(2)*x);
g1 = @(p,x) g(p,x)-counts;
p = lsqnonlin(@(p)g1(p,time),p0);
Local minimum possible. lsqnonlin stopped because the size of the current step is less than the value of the step size tolerance.
hold on
plot(time,counts,'o')
plot([0,time],g(p,[0,time]))
hold off
grid

类别

在 帮助中心 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by