Introduction

Naturally occuring copper is composed principally of two isotopes, Cu63 and Cu65. These isotopes are stable and do not emit radiation. Through thermal neutron bombardment we form isotopes with an additional neutron. After some time we will have produced radioactive Cu63 and Cu65 in differing quantities. Using a simple Geiger Tube, a counter, and a clock we will attempt to determine the half lives of the two isotopes. Numerous methods of analysis could be used. Here I present a Bayesian method using Markov Chain Monte Carlo (MCMC). In essence, this calculates the probability that the parameters of our model would have produced the observed data and searches for the most probable model.

Files for this exercise can be found here:https://github.com/paulnord/bayesian_analysis

The Model

Even with good shielding around our detector and sample, there is some level of background radiation. For the purposes of this model we will assume that the number of background count rate is constant. The total number of observed background counts will simply be a function of the background rate and time.

The total number of observed counts at time t will be a function of five parameters. These include the initial number of atoms for each species a and b, the half lives for each species, and the background count rate.

We model this using the pystan code as follows.

import pystan

import pandas as pd

import pickle

import multiprocessing

multiprocessing.set_start_method("fork")

model = """

data {

int<lower=0> N;

vector[N] a;

int d[N];

}

parameters {

real<lower=0> A0; // Initial number of atoms a

real<lower=0.0, upper=20> t0; // Half Life a

real<lower=0> A1; // Initial number of atoms b

real<lower=400, upper=1000> t1; // Half Life b

real<lower=0, upper=20> background;

}

model {

for(n in 1:N) {

d[n] ~ poisson( A0 + A1 - A0*(0.5^((a[n])/t0)) - A1*(0.5^((a[n])/t1)) + background*a[n] );

}

}

"""

def main():

df = pd.read_csv("my_data.csv")

x = df["t"].to_numpy()

y = df["count"].to_numpy()

# Put our data in a dictionary

data = {'N': len(x), 'a': x, 'd': y}

# Compile the model

sm = pystan.StanModel(model_code=model)

# Train the model and generate samples

fit = sm.sampling(data=data, iter=2000, chains=5, warmup=1000, thin=1, seed=91801, n_jobs=4, control=dict(adapt_delta=0.9))

pickle.dump( sm, open( "model.p", "wb" ) )

pickle.dump( fit, open( "fit.p", "wb" ) )

pickle.dump( data, open( "data.p", "wb" ) )

## Diagnostics #################################################################

print(fit)

if __name__ == '__main__':

main()

Note that the number observed at time t, denoted in the data as d[n], is predicted to be a Poisson distribution around the precise values generated by the mathematical model.

Data

t,count

0.5,87

1,164

1.5,259

2,331

3.5,494

4,551

4.5,604

5,659

5.5,714

6,766

6.5,825

7,875

7.5,924

8,971

8.5,1007

9,1054

9.5,1100

10,1143

10.5,1179

11,1221

11.5,1260

12,1305

12.5,1346

13,1388

13.5,1418

14,1456

14.5,1490

192.3833333,11184

196.8333333,11407

307.15,16679

312.1833333,16919

916.6666667,40583

918.6666667,40660

1086.633333,46063

1089.583333,46138

1271.666667,51484

1273.666667,51538

1446.666667,56114

1448.666667,56166

1629.1,60640

1631.85,60697

1747.716667,63394

1748.716667,63414

Results

A0 738.2109 +/- 37.703469

t0 4.674207 +/- 0.411689 minutes

A1 48328.88 +/- 6742.159141

t1 765.6829 +/- 74.729073 minutes

background 13.99500 +/- 2.212040 counts/minute

Expected values for the two decay rates are well within the uncertainties described here. The plots above show an overlay of 1000 models selected from a sample of the models used in this analysis. They are plotted in a light color so that they form a fuzzy edge around the plotted line described by the mean of the parameters. All of these models generate curves which are very close to the data. While the errors seem very large, they are actually a better representation of the true uncertainty in applying this model to this data. Many least-squares fitting functions will give uncertainties which give too much confidence in the model predictions.

Thanks

I'd like to thank students William Bakke and Anand Agrawal for sharing the data from their lab exercise.