PHY113 Python Lab (Exercise 6)

Here’s the solution:

# Exercise 6: Fit linearised Radioactive Decay
# Transform data and parameters to linear form: Y = A + B*X
X = t         # transform t data for fitting (trivial)
Y = log(N)    # transform N data for fitting
dY = dN/N     # transform uncertainties for fitting

def linearfit(x,m,c):
    return m*x+c

popt,pcov = curve_fit(linearfit, X, Y, sigma=dY)
sigma = sqrt(diag(pcov))

slope = popt[0]
tau = -1/slope
dtau = tau**2 * sigma[0]

intercept= popt[1]
Nfit = exp(intercept)
dNfit = Nfit * sigma[1]

print Nfit,dNfit
print tau,dtau

xfit = linspace(0, 500, 100) 
yfit = slope*xfit + intercept

# Plot Data
errorbar(X, Y, fmt='o',label='Data',yerr=dY)
plot(xfit,yfit, label='Weighted Linear Fit')

xlabel("time [days]")

# Add Legend
legend(loc='best') # best location to avoid hiding data
# Save figure as png file with 300dpi
savefig('DecayLinearisedFit.png', dpi=300, format='png', bbox_inches='tight')

So what’s going on ?

After creating variables for convenience, we define a new function called linearfit. Notice the syntax :

def function_name(<arguments>):
    return <output>

In our case, we don’t need this function to do anything else than evaluate a simple expression, so we can just set it to return a number. Next we’re calling the curve_fit method we loaded from scipy.optimize (see here). This is an algorithm that uses non-linear least squares to fit our function linearfit to the data arrays X and Y. It assumes that Y is a function of X with some parameters to adjust and (possibly) an uncertainty sigma (offset). It returns an array, popt, which contains the optimised parameters of linearfit (m and c) such as to minimize

which is exactly what a $\chi^2$ fit is. The algorithm also returns a 2d-array pcov, which you can think of as a matrix. The diagonal contains the variance of the popt, such that we can extract the one standard deviation errors by taking the square root of those elements. Note that the way we define sigma makes it a 1d-array.

The rest of the exercise should by now be familiar and quite straightforward.

Note : if you’re wondering how curve_fit knows what the parameters and the x-variable are when we only pass linearfit as argument, you’re starting to ask yourself the right questions. Kudos. In fact, curve_fit assumes that the function you pass as argument is of the form function(x,*parameters), which is exactly why we ordered the arguments of linearfit in its definition the way we did.

Solutions to exercise 7

comments powered by Disqus