Estimating λ from Poisson counts

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⁵ decays, 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.

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.

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.

Expand ln L about λ̂. Minus the second-derivative matrix, averaged over data, is the Fisher matrix I_ab = Σᵢ μᵢ⁻¹ (∂μᵢ/∂θₐ)(∂μᵢ/∂θ_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Δ ln L 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 profile interval is [0.70, 1.37] λ̂ while the curvature parabola claims a symmetric [0.67, 1.33]. Both cover close to 68% at that n, so the gap is not a verdict on either; what separates them is that the parabola depends on the variable it is drawn in — built in τ = 1/λ it maps back to [0.75, 1.50] λ̂ — while the Δ ln L = −½ interval is the same whether you scan λ, τ or T½.

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 k_b with mean b T_b — to the likelihood.

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 over five half-lives pins λ to 0.12%. The same data therefore support a half-life quoted to four digits and an activity quoted to two, 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 λ. Coverage is set by the total count, not the count per bin, because the information adds across bins: 400 counts spread over 100 bins cover 68% with either interval. Only when the total itself is small do the two part company, and then the pseudo-experiments, not a formula, decide which one to publish.