Poisson distribution & exponential decay

[cjurban]

Homework Statement

t(s) = 1 15 30 45 60 75 90 105 120 135
N(counts) = 106 80 98 75 74 73 49 38 37 22
Consider a decaying radioactive source whose activity is measured at intervals of 15 seconds. the total counts during each period are given. What is the lifetime τ of the source? What are the uncertainties on N_0 and τ?

Homework Equations

N(t) = N_0*exp(-t/τ)

The Attempt at a Solution

ln[N(t)]=ln[N_0]+(-t/τ) gives an affine function, so if I plot it I can get N_0 and τ from the graph. However this doesn't feel satisfactory to me, I'd prefer to do it mathematically, but I can't figure it out. I don't know how to do a linear fit for a Poisson distributed set of data. I think the error on the N(t) is given by (σ/√n) where n is the number of data points and σ=√avg[N(t)]. Any help?

 

[mfb]

There are formulas how to get the uncertainties for a linear fit, if you know the uncertainties of the individual values. Alternatively, use a computer, fit programs now those formulas.
The individual values follow (approximately) a Poisson distribution.

Strictly speaking, this direct approach will lead to slightly wrong values. This can be seen with datapoints that are zero (not present here, but they are just the most extreme case) - the calculated uncertainty would be zero, which is obviously wrong. There are better methods to treat this properly, but I think that is beyond the scope of this task.

 

[haruspex]

Homework Statement

t(s) = 1 15 30 45 60 75 90 105 120 135
N(counts) = 106 80 98 75 74 73 49 38 37 22
Consider a decaying radioactive source whose activity is measured at intervals of 15 seconds. the total counts during each period are given.

I'm a bit confused. There are ten timestamps, starting at (effectively) 0, so that should be nine intervals, but there are ten counts.

 

[haruspex]

Glossing over the counts versus intervals problem, I tried attacking it from first principles using MLE. It gets nasty.
Suppose there are N potential decay events initially, all with parameter lambda. The prob that a given decay occurs in (t_i, t_i+1) = $e^{-\lambda t_i}-e^{-\lambda t_{i+1}}$.
The prob that k_i occur in each such interval is $\left(\stackrel{N}{k_1, k_2, ... k_r}\right)\prod \left(e^{-\lambda t_{i}}-e^{-\lambda t_{i+1}}\right)^{k_i}$
Taking logs and differentiating wrt lambda:
$\Sigma k_i \frac{t_{i+1}e^{-\lambda t_{i+1}}-t_{i}e^{-\lambda t_{i}}}{e^{-\lambda t_{i}}-e^{-\lambda t_{i+1}}}=0$
Solving that for lambda will not be fun.

 

[mfb]

Why do you want to solve it for lambda? I would use a numerical approximation.

 

[TSny]

Glossing over the counts versus intervals problem, I tried attacking it from first principles using MLE. It gets nasty.
Suppose there are N potential decay events initially, all with parameter lambda. The prob that a given decay occurs in (t_i, t_i+1) = $e^{-\lambda t_i}-e^{-\lambda t_{i+1}}$.
The prob that k_i occur in each such interval is $\left(\stackrel{N}{k_1, k_2, ... k_r}\right)\prod \left(e^{-\lambda t_{i}}-e^{-\lambda t_{i+1}}\right)^{k_i}$
Taking logs and differentiating wrt lambda:
$\Sigma k_i \frac{t_{i+1}e^{-\lambda t_{i+1}}-t_{i}e^{-\lambda t_{i}}}{e^{-\lambda t_{i}}-e^{-\lambda t_{i+1}}}=0$
Solving that for lambda will not be fun.

Interesting calculation.

Since in this problem, $t_{i+1} = t_i + \Delta t$ where $\Delta t$ is fixed at 15 s, it doesn't seem too hard to solve your equation for $\lambda$. [Unless I'm making a stupid mistake.]

 

[haruspex]

Why do you want to solve it for lambda? /QUOTE]
Because that is the MLE for 1/τ in the OP. I note that the question asks for the uncertainty in the estimate, but a first step, surely, is establishing how to make such an estimate.

 

[haruspex]

Interesting calculation.

Since in this problem, $t_{i+1} = t_i + \Delta t$ where $\Delta t$ is fixed at 15 s, it doesn't seem too hard to solve your equation for $\lambda$. [Unless I'm making a stupid mistake.]

You're right - I looked at that too briefly. If the time intervals are all T, looks like it reduces to $e^{-\lambda T} = \frac{x}{x+1}$ where $x = \frac{\Sigma i k_i}{\Sigma k_ i}$. Is that what you get?
Not sure how to go about determining the uncertainty though.

 

[mfb]

I don't think this has to be done analytically, and it is a standard fit problem. Find lambda where the likelihood is maximal, find lambda+, lambda- where the likelihood gets worse by a factor of [I would have to check that].

 

[TSny]

If the time intervals are all T, looks like it reduces to $e^{-\lambda T} = \frac{x}{x+1}$ where $x = \frac{\Sigma i k_i}{\Sigma k_ i}$. Is that what you get?

Yes, that's what I get.

 

[haruspex]

I don't think this has to be done analytically, and it is a standard fit problem.

I don't understand what you mean by a standard fit problem in this context. The decays are counts in an interval, not instantaneous decay rates. Doesn't that complicate things a little? Moreover, the counts in one interval are not independent of those in the next.

 

[mfb]

The decays are counts in an interval, not instantaneous decay rates.

The number of decays in an interval is a binomial distribution where N is the number of atoms at the start and p follows from lambda (actually, we can estimate p and find lambda afterwards - p is all we need).

Moreover, the counts in one interval are not independent of those in the next.

Good point. Well, it looks easy to integrate that in the likelihood. Just calculate it interval by interval. The total number of atoms is a second free parameter, of course.
Once you have the likelihood function, use a standard fit algorithm to find the maximum and the uncertainty on the parameters.