Exercise 2: Nonlinear fitting methods I

PHY224H1F Exercise 2: Nonlinear fitting methods I
Radiation and nonlinear circuits
We continue with the curve fitting from the past exercise, but extend it to nonlinear models.
Radioactive decay is an example of an exponential relation. We will measure the intensity from
a radioactive isotope, and use curve fitting tools to find the half-life.
Background knowledge for exercise 2
Python lists, arrays, numpy, SciPy, PyPlot, curve_fit()
Error Analysis Chi squared (χ
2
), goodness of the fit
Physics Nuclear decay. Review this from your first year text.
1 Introduction
You’ve probabably learned about exponential decay of radioactive materials in the past. The number of
emissions (during a sampling period) depends on the number of atoms of that isotope. Because an emission
event corresponds to a death of an isotope, it follows that the intensity of radiation decays exponentially,
I(t) = I0e
− t
τ , (1)
where I(t) is the intensity of radiation at time t, and τ is the mean lifetime of the isotope. A more intuitive
parameter is the half-life, t1/2
, defined so that
I(t) = I0
1
2
t
t1/2
. (2)
Our interest for this experiment is calculating the half-life of an isotope by measuring the intensity (emissions/sample) of radiation.
The sample you will measure today is a metastable Barium (Ba-137m). It is prepared by flowing an
acid through a generator containing Cesium (Cs-137). Inside the generator, the Cs-137 beta decays (with a
half-life of 30.1 years), into Ba-137m. The acid washes out the Ba-137m which decays by gamma emission
(662 keV) with a half-life of 2.6 minutes. We will measure the rate of emissions using a Geiger counter.
2 Analysis with logarithms
Nonlinear relations become more complicated to analyze in general. There are analytical formulae to find
the parameters in a linear regression, but typically not for a nonlinear regression.
One method for nonlinear regression is transforming the data into a linear model and using linear regression methods. Logarithms are a particularly useful tool for analyzing exponential relations. Consider the
general exponential relation,
y = y0e
ax
. (3)
Taking the logarithm of both sides,
log(y) = ax + log(y0). (4)
Thus the exponential relation for y became a linear relation for log(y). Using a semi-log plot, with a
logarithmic scale for the y-axis turns an exponential into a straight line, where the slope of the line is the
growth/decay rate.
By first transforming the relation, linear regression can be used to find the parameter a. The model
function is the same, f(x) = ax + b, but now, the data to match is zi = log(yi).
Page 1
PHY224H1F Exercise 2: Nonlinear fitting methods I
The transformation of the data also requires an error propagation to get more accurate results from the
regression,
σzi =
 
 
 
 
σyi
yi
 
 
 
 
. (5)
This propagation formula is the standard one for logarithms.
Note that when you are plotting the predictions from your model, you will need to convert the log(y)
back into y so that it is in the correct units and can be compared with the experimental data.
3 Nonlinear regression
Traditionally, the method described in Section 2 was used for analyzing exponential relations, because a
formula could be used to calculate the parameters. However, there is a problem with performing a linear
regression on zi = axi
.
The maximum likelihood of a χ
2
linear model assumes that the measurement errors are normally distributed. When the axes are rescaled, this assumption is potentially invalidated. A related side-effect is that
certain ranges of y are given more weight in the regression.
With the use of computers and nonlinear optimization, it is now possible to perform a least-squares
optimization using a wide range of nonlinear models. Recall that curve_fit() uses least squares to minimize
χ
2 =
X
N
i=1

yi − y(xi)
σi
2
, (6)
where yi
is the dependent data you measured, xi
is the independent data, σi
is the measurement error in the
dependent data. For an exponential model, we simply use y(x) = aebx as the model function.
Computer-aided nonlinear regression brings with it other things to consider like whether the model
function is suitable, if the initial guess for parameters is close enough to the true optimum values, and still
if the measurement errors fit a normal distribution.
4 The experiment
The radioactive sample used in this experiment will be prepared for you by a lab technologist, as a demonstration.
You will not acquire the data on your computer. Instead you’ll find the radioactive decay text file
at: /Desktop/2nd Year lab Files/Exercises/Cesium Radioactive Decay.txt. The instrument used
to acquire the data is called a Geiger counter.
You will also be provided with a text file with data for background radiation: /Desktop/2nd Year lab
Files/Exercises/Radioactive Background 20min 20secdwell.txt
The radioactive energy is quite low and safe. We are using a single demonstration Geiger
counter because of the complexity of setting up 30 identical experiments before the sample
has completely decayed. You still need to document the experiment in the lab–book as if you
were conducting it.
4.1 Uncertainty in the Geiger counter
The Geiger counter follows counting statistics for determining the uncertainty. This is much like counting
the number of heads in a series of coin-flip trials (the √
N rule). The rate is the number of events N divided
by the sample time ∆t. The best estimate of the count rate is given by,
R =
N
∆t
±
√
N
∆t
. (7)
Page 2
PHY224H1F Exercise 2: Nonlinear fitting methods I
There is also background radiation to consider; you should have heard the counter beeping with no sample
present. The counts from the two sources of radiation add together for our final reading,
Ns = NT − Nb . (8)
The uncertainty can be calculated via standard error propagation, so
σs =
q
σ
2
T + σ
2
b
=
p
NT + Nb .
Finally, because the background radiation cannot be measured at the same time, we assume it has constant
statistics, and use its mean in the calculations (Nb = N¯
b). The mean background radiation must be
subtracted from each data point before analysis.
5 The Python program
We will use our programs to compare the transformation method (Section 2) and the non-linear leastsquares method (Section 3). The best parameters will be found with both methods, and used to plot best
fit curves over the data. For this experiment, we have the luxury of comparing our calculated half-time to
the theoretical value.
The program should be organized as follows:
• Import the required functions and modules (numpy, scipy.optimize.curve_fit(), matplotlib.pyplot).
• Define the model functions (f(x, a, b) = ax + b, g(x, a, b) = b eax)
• Load the data using loadtxt().
• Subtract the mean background radiation from the data.
• Calculate the standard deviation for each data point.
• Convert the count data into rates (divide by ∆t).
• Perform the linear regression on (xi
, log(yi)) using f as the model function.
• Perform the nonlinear regression on (xi
, yi) using g as the model function.
• Output both of the half-life values you calculated.
• Plot the errorbars, both curves of best fit (linear and non-linear), and the theoretical curve. Include
a legend for the different curves.
• Make a second plot of the same data, but using a logarithmic y axis. One way of doing this is to
call the pylab.semilogy function. Another way is to call pylab.yscale(‘log’) after you make the
plot.
An alternate form of the non–linear g function can be used where the exp is replaced by 1
2
or 2. You can
also treat the linearized model f as a natural logarithm (ln), or a base–2 logarithm (log2
). These transforms
might change the way you process your data for input into the function, or change the way you interpret
the errors from the curve_fit.
Write the program, and run it using the data gathered in the experiment. Save all plots and the parameters
you determined.
Which regression method gave a half-life closer to the expected half-life of 2.6 minutes?
Can you see the difference on the plots?
6 Analyzing the quality of the fit
There are two ways that we can assess the quality of the fit of our model: variance of the calculated
parameters, and the reduced chi-squared statistic.
Page 3
PHY224H1F Exercise 2: Nonlinear fitting methods I
6.1 Variance of parameters
The variance of the parameters is returned by curve_fit() as the diagonal entries in the covariance matrix,
p_cov. Recall that the error in measurements is understood as the standard deviation,
t1/2 = t˜1/2 ± σt1/2
, (9)
where σt1/2 =
p
Var(t1/2
). The variance of the parameter a is p_cov[0, 0].
To get σt1/2
from σa, you will need to use the error propagation formula for a reciprocal,
σ

1
a

=
σa
a
2
. (10)
Modify your program from Section 5 to calculate the standard deviation of the half-life.
What values did you find? Does the nominal half-life (2.6 minutes) fall in the range?
6.2 Reduced χ
2
Recall that the χ
2 distribution gives a measure of how unlikely a measurement is. The more a model deviates
from the measurements, the higher the value of χ
2
. But if χ
2
is too small, it’s also an indication of a problem:
usually that there weren’t enough samples.
This best (most likely) value of χ
2 depends on the number of degrees of freedom of the data, ν = N − n,
where N is the number of observations and n is the number of parameters. This dependence is removed
with the reduced chi-squared distribution,
χ
2
red =
1
ν
X
N
i=1

yi − y(xi)
σi
2
. (11)
For χ
2
red, the best value is 1:
χ
2
red  1 indicates that the model is a poor fit;
χ
2
red > 1 indicates an incomplete fit or underestimated error variance;
χ
2
red < 1 indicates an over-fit model (not enough data, or the model somehow is fitting the noise).
Add a function to your program to calculate χ
2
red. What values were computed? How would
you interpret your results?
Prepared by John Ladan and Ruxandra Serbanescu, 2016. Edited by Christopher Lee, 2017.
Page 4