Skip to main content
GIOPHYSICS
Menu
← University Physics V

University Physics V · Radioactive Decay · 15.7

Estimating λ from Poisson Counts

Counts scatter because decay is random, not because the detector is sloppy. This is how to turn a column of scattered integers into a decay constant with an interval you can defend: maximise the likelihood the Poisson process actually implies, and put background and dead time inside the model rather than on the data.

Open the lab ↓Practice & worked examples ↓
01

Build the model

Connect the measurement to the mechanism.

What you actually hold is a list of integers — k₁, k₂, … counts in successive timed bins — each a draw from a Poisson distribution whose mean the model supplies, μᵢ = S e(−λ tᵢ) + bΔ. That Poisson is not an assumption of convenience: it follows from independent decay, binomial survival in each window, thinning by a detection efficiency, and the addition of an independent background, and it is the only honest description at the end of a run, where the model predicts three counts and the detector records zero. Estimating λ therefore means maximising Σ(kᵢ ln μᵢ − μᵢ), a smooth but non-linear function of (λ, S, b) with no closed form, and taking the interval from the curvature of that surface or from a profile scan across λ.

The cost is that nothing can be read off a straight line any more. The gain is threefold: empty bins still speak, efficiency and solid angle cancel out of λ̂ because they enter only through the amplitude S, and every systematic left out of μ — a dead-time loss that tilts the early bins, a background that cannot be subtracted through a logarithm — appears as a bias that more counting makes look worse, not better.

Simple definition
Maximum-likelihood estimation of a decay constant treats each timed bin's count as a Poisson draw whose mean is the model's integrated activity, and chooses the λ that makes the observed set of integers most probable.
Example
A bin whose model mean is μ = 4.0 counts contributes k ln μ − μ = 1 × 1.386 − 4.0 = −2.61 when k = 1, and still contributes −4.00 when k = 0 — precisely the term a log-linear fit has to discard.
Poisson mean of a timed binμᵢ = S e(−λ tᵢ) + b Δ, S = ε N₀ (1 − e(−λΔ))

ε and N₀ enter only through S, so letting the amplitude float makes λ̂ blind to both.

tᵢ bin start (s), Δ bin width (s), b background rate (s⁻¹), ε efficiency, counts dimensionless

Poisson log-likelihoodln L(λ, S, b) = Σᵢ [ kᵢ ln μᵢ − μᵢ ] + const

A bin with kᵢ = 0 still contributes −μᵢ. Logging the data deletes that term.

kᵢ are integer counts; the dropped −Σ ln kᵢ! carries no parameter

Where the maximum sits (b = 0)Σ kᵢ tᵢ / Σ kᵢ = Σ tᵢ e(−λ tᵢ) / Σ e(−λ tᵢ)

λ̂ matches the observed mean count time to the model's — a 1-D root find, not a linear regression.

Equal bins; the amplitude profiles out as Ŝ = Σ kᵢ / Σ e(−λ̂ tᵢ)

Fisher matrix and the covarianceIab = Σᵢ μᵢ⁻¹ (∂μᵢ/∂θₐ)(∂μᵢ/∂θb), Cov ≈ I⁻¹

σλ = √[(I⁻¹)λλ] ≥ 1/√Iλλ: floating S and b inflates the error by 1/√(1 − ρ²).

θ = (λ, S, b). I is a Gram matrix of derivative vectors under the weight 1/μᵢ

What sets the precision on λσλ = 1 / √( K · Varw t ), wᵢ = μᵢ / K

Counts alone are not enough — you need a lever arm in t. Unbinned and complete, Varw t = 1/λ² and σλ = λ/√n.

K expected total counts; Varw t is the variance of arrival time under the model weights (s²)

Dead time, before anything elsenon-paralysable m = n/(1 + n τd) · paralysable m = n e(−n τ_d)

Losses grow with rate, so they flatten the early bins and pull λ̂ low — a bias no extra counting removes.

m observed rate, n true rate (both s⁻¹), τd dead time (s)

01

Why a bin count is exactly Poisson

Start from the boundary condition, not the formula. N₀ identical nuclei decay independently, each surviving to t with probability e(−λt), so the number decaying in the window [tᵢ, tᵢ + Δ] is Binomial(N₀, pᵢ) with pᵢ = e(−λ tᵢ)(1 − e(−λΔ)). For a microcurie of Cs-137, N₀ ≈ 5 × 10¹³ and a 10 s bin has pᵢ ≈ 7 × 10⁻⁹: that is the Poisson limit N₀ → ∞, pᵢ → 0 at fixed N₀ pᵢ = 3.7 × 10⁵ counts, and the binomial variance N₀ p(1 − p) collapses onto the mean. Detection is independent thinning — keep each decay with probability ε — and thinning maps Poisson(μ) to Poisson(εμ). Background arrives as its own independent Poisson stream, and independent Poissons add. So the gross count in a bin is Poisson with mean μᵢ = S e(−λ tᵢ) + bΔ, where S = ε N₀(1 − e(−λΔ)). Two things break the argument, and both return later: dead time, which makes successive counts non-independent, and pile-up, which merges two pulses into one.

02

The likelihood, and what taking logs destroys

Because the kᵢ are independent, ln L(λ, S, b) = Σᵢ [kᵢ ln μᵢ − μᵢ], up to a −Σ ln kᵢ! that carries no parameter and so cannot move the maximum. Maximise that directly and nothing is transformed. The tempting alternative — plot ln kᵢ against tᵢ and fit a straight line — transforms the data instead of the model, and pays four times. Since ln is concave, E[ln k] ≈ ln μ − 1/(2μ), so at μ = 3 counts the plotted point sits at e(−1/6) = 0.85 of the model value, and the bias grows precisely in the late bins that set the slope. Bins with kᵢ = 0 have no logarithm and get dropped, which censors only downward fluctuations and tilts the tail the other way. Var(ln k) ≈ 1/μ is not constant, so an unweighted line lets a 4-count bin speak as loudly as a 4000-count one. And ln(kᵢ − bΔ) is undefined for every bin that fluctuates below background. The likelihood keeps the zeros: a bin with kᵢ = 0 still contributes −μᵢ, a genuine statement that the model predicted counts you did not see.

03

Solving for λ̂: the amplitude profiles out

Set b = 0 for a moment and differentiate. ∂lnL/∂S = Σ (kᵢ/μᵢ − 1) e(−λ tᵢ) = 0 gives Ŝ = Σ kᵢ / Σ e(−λ tᵢ): the amplitude only normalises the model to the total observed counts. Substituting it into ∂lnL/∂λ = 0 leaves one scalar equation, Σ kᵢ tᵢ / Σ kᵢ = Σ tᵢ e(−λ tᵢ) / Σ e(−λ tᵢ) — the observed mean count time must equal the model's mean count time on the same bin grid. That is a moment condition, solvable with Brent's method in a line of SciPy, and it shows why ε and N₀ never appear: they live entirely inside Ŝ. With a background floating there is no such reduction. Minimise −lnL over (λ, ln S, ln b) with scipy.optimize.minimize, using the logs of the positive parameters so the optimiser cannot step to a negative mean, supply the analytic gradient, and confirm the Hessian at the optimum is positive definite instead of trusting that the routine reported success.

04

The interval is a curvature, and the nuisances inflate it

Expand ln L about λ̂. The second-derivative matrix, averaged over data, is the Fisher matrix Iab = Σᵢ μᵢ⁻¹ (∂μᵢ/∂θₐ)(∂μᵢ/∂θb) — a Gram matrix of the model's derivative vectors under the weight 1/μᵢ, hence positive semidefinite, and singular exactly when two parameters are degenerate. The covariance is I⁻¹, and the number you quote is σλ = √[(I⁻¹)λλ], never 1/√Iλλ: the diagonal of the inverse is not the inverse of the diagonal. With one nuisance parameter the two differ by 1/√(1 − ρ²), and for a long run with a floating amplitude ρ² = ½ exactly, so letting S float costs a factor √2 — twice the counts to buy it back. For the interval itself, Wilks' theorem makes 2ΔlnL a χ²₁ variable, so scan λ on a grid, re-maximise the nuisances at each point, and read off where Δ ln L = −½. With n = 9 recorded decays the exact profile interval is [0.70, 1.37] λ̂ while the curvature parabola claims a symmetric [0.67, 1.33]; the skew is real, and only the scan sees it.

05

Background and dead time belong in μ, not in the data

Both corrections precede any statement about λ, and both belong inside the model. A non-paralysable detector with τd = 4.0 µs loses 7.4% of a 20 kHz true rate but only 2.0% of a 5 kHz one, so the loss is a tilt, not a scale factor: left uncorrected it flattens the early bins and biases λ̂ low by about 4% across two half-lives. Inverting the data instead — replacing kᵢ by kᵢ/(1 − mᵢ τd) — gets the mean right and destroys the statistics, because the rescaled counts are no longer integers and their variance is now the squared correction times the mean, while the likelihood still assumes variance equals mean. Fold the live-time fraction into μᵢ instead. Background is the same argument: kᵢ − b̂Δ is a difference of Poissons, Skellam-distributed and free to go negative, and Σ(kₙₑₜ ln μ − μ) is the likelihood of nothing at all. Keep the gross counts, let b float, and add the calibration run's own Poisson term — counts kb with mean b Tb — to the likelihood.

06

λ is efficiency-free; the activity is not

The fit returns Ŝ = ε N₀ (1 − e(−λ̂Δ)), so an absolute activity A₀ = λ̂ N₀ = λ̂ Ŝ / [ε (1 − e(−λ̂Δ))] needs the efficiency, the solid angle, absorption in the source and entrance window, and the branching fraction of the line being counted. On a bench NaI setup those are known to perhaps 5%, while a run of 1.2 × 10⁶ counts pins λ to 0.12%. The same data therefore support a half-life quoted to three digits and an activity quoted to one, and no amount of extra counting shifts that division, because it is systematic rather than statistical. Finish by testing the interval rather than trusting it: draw 10⁴ pseudo-experiments from the fitted μᵢ with a Poisson sampler, refit each, and count how often the Δ ln L = −½ interval covers the input λ. At a few counts per bin the curvature interval undercovers and the profile interval does not, which is how you decide which one to publish.

02

Change one variable at a time

Make the relationship visible.

Interactive model
9
1.50

Set n = 4 and drag the trial value to 1.5, where the dashed parabola crosses −½: the exact curve is still above that line, so the real 68% interval runs further up in λ than the curvature admits. Push n to 60 and the two curves merge — the large-sample limit arriving.

Interactive physics modelThe exact log-likelihood profile Δ ln L(λ) = n [ln u − u + 1] with u = λ/λ̂, for n = 9 recorded decay times (solid), against the parabola −n(u−1)²/2 built from the Fisher curvature (dashed). Zero sits at the top of the frame and the dashed horizontal is Δ ln L = −½. At the trial value u = 1.50 the exact profile reads −0.851 and the parabola −1.125.Δ ln L(λ) = n [ ln u − u + 1 ], u = λ / λ̂, n = 9 decayssolid: exact profile dashed: Fisher parabola −n(u−1)²/2 0 at topΔ ln L = −½00.51.0 λ̂1.5λ / λ̂

Δ lnL EXACT-0.851

Δ lnL PARABOLA-1.125

PARABOLIC σ/λ = 1/√n0.333

WILKS √(2|Δ lnL|)1.30 σ

Live interpretationΔ lnL EXACT: −0.851. Δ lnL PARABOLA: −1.125. PARABOLIC σ/λ = 1/√n: 0.333. WILKS √(2|Δ lnL|): 1.30 σ

03

Catch the common trap

Explain before calculating.

A 40-bin gamma run on an eluted Ba-137m sample has a separately measured background, and several late bins read zero counts. Which procedure estimates λ and its 68% interval correctly?

Choose an answer to test the model.

04

Practice & worked examples

Reason from the model, then test the result.

EasyA gated timer records the individual arrival times of n = 400 counts from an eluted Ba-137m sample over a window long enough that truncation is negligible, and their sum is Σ tⱼ = 8.72 × 10⁴ s. Estimate λ, its standard error, and the half-life with an uncertainty, then compare with the tabulated T½ = 153.0 s.
  1. For a complete unbinned observation, ln L = n ln λ − λ Σ tⱼ, so ∂lnL/∂λ = n/λ − Σ tⱼ = 0 and λ̂ = n / Σ tⱼ: the mean arrival time is the whole sufficient statistic.
  2. λ̂ = 400 / (8.72 × 10⁴ s) = 4.587 × 10⁻³ s⁻¹, so τ̂ = 218.0 s.
  3. Fisher information: I(λ) = −E[∂²lnL/∂λ²] = n/λ², so σλ = λ̂/√n = 4.587 × 10⁻³ / 20 = 2.29 × 10⁻⁴ s⁻¹, a relative error of 1/√400 = 5.0%.
  4. T̂½ = ln2 / λ̂ = 0.6931 / (4.587 × 10⁻³ s⁻¹) = 151.1 s. The half-life is a fixed multiple of 1/λ, so it inherits the same 5.0%: ±7.6 s.
  5. The estimator is biased: E[λ̂] = n λ/(n − 1), so the corrected value is (n − 1)/Σ tⱼ = 4.576 × 10⁻³ s⁻¹, 0.25% lower. That is one twentieth of σλ, and since the ratio √n/(n − 1) falls as n grows this bias never takes over.

Answerλ̂ = (4.59 ± 0.23) × 10⁻³ s⁻¹ and T½ = 151 ± 8 s, consistent with 153.0 s. The 1/(n − 1) bias correction moves λ̂ by 0.25%, one twentieth of the statistical error.

MediumTwo runs on the same nuclide each collect K = 40 000 counts, one over a window of a single half-life and one over five. Using σλ = 1/√(K Varw t) with weights w(t) ∝ e(−λt) on [0, T], find the relative error on λ for each and say what the short run is actually short of.
  1. Put a = λT and E = e(−a). With f₀ = 1 − E, f₁ = 1 − (1 + a)E and f₂ = 2 − (a² + 2a + 2)E, the weighted moments give Varw t = (f₀f₂ − f₁²)/(λ²f₀²), so σλ/λ = f₀ / [√K √(f₀f₂ − f₁²)].
  2. One half-life: a = ln2 = 0.6931, E = 0.5, so f₀ = 0.5000, f₁ = 0.1534, f₂ = 0.06663 and f₀f₂ − f₁² = 0.03331 − 0.02354 = 0.009773, whose root is 0.09886.
  3. σλ/λ = 0.5000 / (200 × 0.09886) = 2.53%.
  4. Five half-lives: a = 3.4657, E = 1/32 = 0.03125, so f₀ = 0.96875, f₁ = 0.86045, f₂ = 1.34554 and f₀f₂ − f₁² = 1.30349 − 0.74037 = 0.56312, whose root is 0.75042; σλ/λ = 0.96875 / (200 × 0.75042) = 0.65%.
  5. The T → ∞ limit gives Varw t = 1/λ² and σλ/λ = 1/√K = 0.50%. Same counts, a factor 3.9 in precision: the short run is starved of lever arm, not of counts, because every arrival is packed into a window over which e(−λt) has only fallen to ½.

Answerσλ/λ = 2.53% over one half-life and 0.65% over five, against an ideal 1/√K = 0.50%. Equal counts, a factor of 3.9 in precision — the missing ingredient is the width of the window.

HardA NaI detector with non-paralysable dead time τd = 4.0 µs watches an eluted Ba-137m source. The observed total rate is 18 529 s⁻¹ in the first bin and 4 913.5 s⁻¹ in a bin 306.0 s later; a separate run gives a background of 12 s⁻¹. The whole run records K = 1.2 × 10⁶ counts over five half-lives. Find λ̂ with and without the dead-time correction, and express the difference in units of the statistical error.
  1. Invert the non-paralysable model, n = m/(1 − m τd): n₁ = 18 529/(1 − 0.074116) = 2.0012 × 10⁴ s⁻¹ and n₂ = 4 913.5/(1 − 0.019654) = 5.012 × 10³ s⁻¹. The correction adds 8.0% at the start of the run and only 2.0% two half-lives later — a tilt, not a common factor.
  2. Subtract the background from the corrected rates: s₁ = 20 012 − 12 = 2.0000 × 10⁴ s⁻¹ and s₂ = 5 012 − 12 = 5.000 × 10³ s⁻¹.
  3. Two-point estimate: λ̂ = ln(s₁/s₂)/Δt = ln 4.000 / 306.0 s = 1.38629/306.0 = 4.530 × 10⁻³ s⁻¹, so T½ = ln2/λ̂ = 153.0 s.
  4. Skip the dead-time step and the same arithmetic on the raw rates gives ln(18 517/4 901.5)/306.0 = ln 3.7778/306.0 = 1.32913/306.0 = 4.344 × 10⁻³ s⁻¹, i.e. T½ = 159.6 s — λ low by 4.1%, T½ high by 4.3%.
  5. Statistical error over five half-lives, using the coefficient from the previous problem: σλ/λ = 1.291/√K = 1.291/1095 = 0.118%.
  6. The uncorrected answer therefore sits 4.12/0.118 = 35 standard deviations from the corrected one. At 10⁶ counts the measurement is limited entirely by the dead-time model, and more counting only makes the discrepancy more significant.

AnswerCorrected: λ̂ = 4.530 × 10⁻³ s⁻¹, T½ = 153.0 s. Uncorrected: λ̂ = 4.344 × 10⁻³ s⁻¹, T½ = 159.6 s — 4.1% low in λ, a 35σ error against the 0.12% statistical uncertainty.