Force a starting point on exponential graph

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!

Answer 1

[Davide Masiello]

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)

[Ran Kagan]

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]

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]

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]

@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]

@Ran Kagan I understand a bit better now, thanks. I have modified my answer accordingly, see if it helps.

[Ran Kagan]

@Davide Masiello Thanks!

[Davide Masiello]

My pleasure, if that's the answer you were looking for don't forget to accept it!

Answer 2

[Torsten]

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