From bb81835032931d0295f400518044905dcb36873b Mon Sep 17 00:00:00 2001
From: Matthew Pitkin <matthew.pitkin@ligo.org>
Date: Thu, 2 Aug 2018 23:16:25 +0100
Subject: [PATCH] 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 | 69 +++++++++++---------
1 file changed, 39 insertions(+), 30 deletions(-)
diff --git a/examples/other_examples/radioactive_decay.py b/examples/other_examples/radioactive_decay.py
index 0767632e0..219f7bc58 100644
--- a/examples/other_examples/radioactive_decay.py
+++ b/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()
--
GitLab