Radioactive decay (or radioactivity) is the process by which a nucleus of an unstable atom loses energy by emitting particles such as alpha particles (Helium nucleii) or beta particles (electrons). You can read up more about radioactive decay here.
In this exercise, we look at a (fake) dataset showing the rate of detection of alpha particles close to a radioactive sample, which we will call element X. The dataset used is the data/new_decay_data.txt file. The rate of detection of alpha particles is measured at 100 different times (given by the first column, in days since the start of the experiment), and for each of these times, both the rate of detection and the uncertainty in the rate of detection are given (in the second and third columns, in detections per second).
1 - Plot the rate of detection versus time, including error bars (see the lecture notes for an example of how to do this). Include labels on the x and y axis.
2 - The rate of decay of a radioactive element is given by:
$$R(t) = R_0~e^{-t/\tau}$$where $R_0$ is the rate of decay at time $t=0$, and $\tau$ is the mean lifetime of the element. Define this function (which has two free parameters in addition to the time) as a Python function in your script, and fit it to the data using curve_fit, and overplot the resulting best-fit function over the data using a line. Is the fit good? If not, what do you think is causing the fit to be bad? Print out the best-fit parameters and uncertainties as shown in the lecture notes.
3 - While one would ideally be able to measure the rate of decay in a completely isolated environment, some of the alpha detections are in fact from background radioactive decay, so we need to take this into account by considering a decay rate function that includes a constant background level term:
$$R(t) = R_0~e^{-t/\tau} + R_{\rm bkg}$$Define this function in your script and fit it to the data using curve_fit, overplotting the best-fit function over the data in a new plot. Also draw a horizontal line to show the constant background level. Does this fit look better? Print out the best-fit parameters and uncertainties.
4 - One way to quantify the goodness of a fit is to consider the reduced $\chi^2$ value of the fit, defined as:
$$\chi_{\rm red}^2 = \frac{1}{N - p}\sum_{i=1}^{N}~\left(\frac{y_i - m_i}{\sigma_i}\right)^2$$where $N$ is the number of datapoints, $p$ is the number of parameters for the model, $y_i$ are the data values, $m_i$ are the model values at the same positions, and $\sigma_i$ are the uncertainties on the data $y_i$.
Conventional wisdom has it, that a fit is usually considered to be good if $\chi_{\rm red}^2 \approx 1$, and this is a decent guide.
Compute $\chi_{\rm red}^2$ for both models above, and comment on the values. Do they agree with what you would have concluded from looking at the plots?
# your solution here
A note on the $\chi^2$ distribution: One can and should look at this more precisely: given the $\chi^2$ value
$$\chi^2 = \sum_{i=1}^{N}~\left(\frac{y_i - m_i}{\sigma_i}\right)^2$$and the number of degrees of freedom df (here df=N-p), one can calculate the probability that a chisquare $\chi^2$ greater than the one found could have been obtained:
# example using the cumulative chisquare distribution implemented in scipy.stats
from scipy.stats import chi2
chisquare=20 # enter chisquare value
df=10 # enter dof (is this is low then the probability is lowered)
print (1-chi2.cdf(chisquare,df=df))
You can indeed look at the chisquare probability function with matplotlib, for which the expectation value is $<\chi^2> =$ df.
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
df = 20
x = np.arange(0,2*df,0.1)
# Generate random numbers
r = chi2.rvs(df, size=1000)
plt.hist(r, density=True, histtype='stepfilled', alpha=0.2)
# Display probability density function (pdf) and compare with histogram
plt.plot(x,chi2.pdf(x,df=df),lw=3)
plt.xlabel('$\chi^2$')
plt.ylabel('chi2.pdf')
plt.title('df = '+str(df))
For more information check out the SciPy manual at scipy.stats.chi2