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 tau = -1/slope dtau = tau**2 * sigma intercept= popt Nfit = exp(intercept) dNfit = Nfit * sigma 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]") ylabel("ln($N$)") # 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>): <operations> 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
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
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.