radioactive_decay.py: complete example

 - complete the radioactive decay example to test the Poisson
   likelihood
 - refs Monash/tupak!132

examples/other_examples/radioactive_decay.py

@@ -2,7 +2,7 @@
 """
 An example of how to use tupak to perform paramater estimation for
 non-gravitational wave data. In this case, fitting the half-life and
-initial radionucleotide number for Carbon 14. 
+initial radionucleotide number for Lead 214. 
 """
 from __future__ import division
 import tupak
@@ -17,56 +17,65 @@ label = 'radioactive_decay'
 outdir = 'outdir'
 tupak.utils.check_directory_exists_and_if_not_mkdir(outdir)
 
-NA = 6.02214078e23 # Avagadro's number
+# generate a set of counts per minute for Polonium 214 with a half-life of 20 mins
+halflife = 20
+N0 = 1.e-19 # initial number of radionucleotides in moles
 
-# generate a set of counts per minute for Carbon 14 with a half-life of
-# 5730 years
-halflife = 5730.
-N0mol = 1. # initial number of radionucleotides in moles
-N0 = N0mol * 
+def decayrate(deltat, halflife, N0):
+    """
+    Get the decay rate of a radioactive substance in a range of time intervals
+    (in minutes). halflife is in mins. N0 is in moles.
+    """
 
+    times = np.cumsum(deltat) # get cumulative times
+    times = np.insert(times, 0, 0.)
 
-# First, we define our "signal model", in this case a simple linear function
-def model(time, m, c):
-    return time * m + c
+    ln2 = np.log(2.)
+    NA = 6.02214078e23 # Avagadro's number
 
+    rates = (N0*NA*(np.exp(-ln2*(times[0:-1]/halflife)) 
+             - np.exp(-ln2*(times[1:]/halflife)))/deltat)
+
+    return rates
 
-# Now we define the injection parameters which we make simulated data with
-injection_parameters = dict(m=0.5, c=0.2)
 
-# For this example, we'll use standard Gaussian noise
-sigma = 1
+# Now we define the injection parameters which we make simulated data with
+injection_parameters = dict(halflife=halflife, N0=N0)
 
 # These lines of code generate the fake data. Note the ** just unpacks the
 # contents of the injection_parameters when calling the model function.
-sampling_frequency = 10
-time_duration = 10
-time = np.arange(0, time_duration, 1/sampling_frequency)
-N = len(time)
-data = model(time, **injection_parameters) + np.random.normal(0, sigma, N)
+sampling_frequency = 1.
+time_duration = 30.
+time = np.arange(0, time_duration, 1./sampling_frequency)
+deltat = np.diff(time)
+
+rates = decayrate(deltat, **injection_parameters)
+# get radioactive counts
+counts = np.array([np.random.poisson(rate) for rate in rates])
 
 # We quickly plot the data to check it looks sensible
 fig, ax = plt.subplots()
-ax.plot(time, data, 'o', label='data')
-ax.plot(time, model(time, **injection_parameters), '--r', label='signal')
+ax.plot(time[0:-1], counts, 'o', label='data')
+ax.plot(time[0:-1], decayrate(deltat, **injection_parameters), '--r', label='signal')
 ax.set_xlabel('time')
-ax.set_ylabel('y')
+ax.set_ylabel('counts')
 ax.legend()
 fig.savefig('{}/{}_data.png'.format(outdir, label))
 
 # Now lets instantiate a version of the Poisson Likelihood, giving it
-# the time, data and signal model
-likelihood = PoissonLikelihood(time, counts, model)
+# the time intervals, counts and rate model
+likelihood = PoissonLikelihood(deltat, counts, decayrate)
 
-# From hereon, the syntax is exactly equivalent to other tupak examples
-# We make a prior
+# Make the prior
 priors = {}
-priors['m'] = tupak.core.prior.Uniform(0, 5, 'm')
-priors['c'] = tupak.core.prior.Uniform(-2, 2, 'c')
+priors['halflife'] = tupak.core.prior.LogUniform(1e-5, 1e5, 'halflife',
+                                                 latex_label='$t_{1/2}$')
+priors['N0'] = tupak.core.prior.LogUniform(1e-25, 1e-10, 'N0',
+                                           latex_label='$N_0$')
 
 # And run sampler
 result = tupak.run_sampler(
     likelihood=likelihood, priors=priors, sampler='dynesty', npoints=500,
-    walks=10, injection_parameters=injection_parameters, outdir=outdir,
-    label=label)
+    nlive=1000, walks=10, injection_parameters=injection_parameters,
+    outdir=outdir, label=label)
 result.plot_corner()