Re: [Phys-L] Bayesian Inference in Half-Life measurement

[Francois Primeau <fprimeau@uci.edu>]

I'm assuming the problem is to estimate the half-life of a sample given
count data taken at different times.  To be concrete I'm going to assume
that you have counts, c1, c2, and c3 taken at times t1, t2, and t3. Also,
I'm going to assume that the counting was done for 10 minutes, so
dt1=dt2=dt3 = 10*60 seconds.  I'm going to assume that the error in the
data is negligible and that only the Poisson statistics matter.  The Matlab
code and the comments in it explain the solution. The code starts off by
generating synthetic data and then using the synthetic data to estimate the
half-life.

The Bayesian solution (given in Matlab code) is the following:

tau_true = 52*60*60; %  half-life of neutron activated copper in seconds
No = 1e6; % 1 million radioactive atoms at t0 = 0

% all times are in seconds
t1 = 0;        dt1 = 10*60; % count for 10 minutes at t1
t2 = 24*60^2;  dt2 = 10*60; % count for 10 minutes at t2
t3 = 48*60^2;  dt3 = 10*60; % count for 10 minutes at t3

% generate synthetic data
c1 = poissrnd( log(2)*No*exp(-log(2)*t1/tau_true)*dt1/tau_true )
c2 = poissrnd( log(2)*No*exp(-log(2)*t2/tau_true)*dt2/tau_true )
c3 = poissrnd( log(2)*No*exp(-log(2)*t3/tau_true)*dt3/tau_true )

% Now pretend we don't know tau and want to use Bayes' theorem to estimate
it
% Poisson probability function: prob(c) = mu^c*exp(-mu)/c!
% log likelihood function for one of the c's (+ an irrelevant constant that
doesn't depend on tau )
ll = @(tau,c,t,dt) c * log( log(2)*No*exp(-log(2)*t./tau)*dt./tau ) -
log(2)*No*exp(-log(2)*t./tau)*dt./tau;

% The joint probability factors because the probabilities for the counts
*conditioned on tau* are independent
% so the log-likelihood for the full data set is given by the sum of the
individual log-likelihoods
% log-likelihood for the full data set:
LL = @(tau) ll(tau,c1,t1,dt1) + ll(tau,c2,t2,dt2) + ll(tau,c3,t3,dt3);

% I'm playing dumb and assuming the only thing I know is that tau is a
% positive number. In practice, I could use a more informative prior,
% but if I collect even a moderate amount of good data, then the prior won't
% matter much.
% prior: prob(log(tau)) = constant ==> prob(tau) propto 1/tau
lp = @(tau)  1./tau; % log prior (not normalized!)

% log posterior from Bayes' Theorem (not normalized so there is a missing
constant)
LP = @(tau) LL(tau)+lp(tau);

% now make a plot of the posterior
tau = linspace(1.5e5,2.5e5,100); % s
P = exp(LP(tau)-max(LP(tau))); % prevent overflow
Z = trapz(tau,P); % normalization
P = P/Z; %normalized posterior
plot(tau/60,P*60); grid on
xlabel('half life (minutes)')
ylabel('probability density (1/minutes)')
% the plot shows a nice bell shaped posterior
%so the mean and variance will provide a good summary of the posterior

% compute the mean of the posterior
tau_hat = trapz(tau,P.*tau);
% compute the variance of the posterior
sig2 = trapz(tau,P.*(tau-tau_hat).^2);
% print the result
fprintf('tau = %f +/- %f minutes \n',tau_hat/60,sqrt(sig2)/60)



On Tue, Sep 21, 2021 at 10:29 AM Paul Nord <Paul.Nord@valpo.edu> wrote:

> That title doesn't google well.  It seems like there should be a good
> reference for this.  The best I've got is David MacKay's code here:
> http://www.inference.org.uk/mackay/itprnn06/slides/10/mgp00010.html
> Translated into python below.
>
> Two questions:
> 1. If we record the number of counts observed with a geiger tube at
> discrete periods within the decay, is this approach still valid?  Say that
> I've got a sample with a 52 hour half life.  I come back about once a day,
> fire up the good old geiger tube and measure the activity for 10 minutes.
> Can I just use that number of counts as the power of a probability function
> to multiply through here?
>
> 2. Will this give me a good result if I'm extracting multiple decay periods
> from the same data?
>
> The half lives of neutron activated copper is the experiment, actually.
>
> Paul
>
>
> import matplotlib.pyplot as plt
>
> import numpy as np
>
> import math
>
>
> a = 1.0
>
> b = 20.0
>
>
> def p(x,l):
>
>     val = math.exp( -x/l ) / Z(l)/l
>
>     return val
>
>
> def Z(l):
>
>     val = math.exp( -a/l) - math.exp(-b/l)
>
>     return val
>
>
> def like(l):
>
>     return p(3,l) * p(5,l) * p(12,l)
>
>
>
> t1 = np.arange(0.01,1000.0,.1)
>
>
>
> fig,axs = plt.subplots(2)
>
> axs[0].semilogx()
>
> axs[1].semilogx()
>
>
> axs[0].plot(t1,[p(3,x) for x in t1])
>
> axs[0].plot(t1,[p(5,x) for x in t1])
>
> axs[0].plot(t1,[p(12,x) for x in t1])
>
>
> axs[1].plot(t1,[like(x) for x in t1])
>
>
> plt.show()
> _______________________________________________
> Forum for Physics Educators
> Phys-l@mail.phys-l.org
> https://www.phys-l.org/mailman/listinfo/phys-l


-- 
Professor of Earth System Science
Rm 3224 Croul Hall
University of California, Irvine