4

I have some radioactive decay data, which has uncertainties in both x and y. The graph itself is all good to go, but I need to plot the exponential decay curve and return the report from the fitting, to find the half-life, and reduced chi^2.

The code for the graph is:

 fig, ax = plt.subplots(figsize=(14, 8))
    ax.errorbar(ts, amps, xerr=2, yerr=sqrt(amps), fmt="ko-", capsize = 5, capthick= 2, elinewidth=3, markersize=5)
    plt.xlabel('Time  /s', fontsize=14)
    plt.ylabel('Counts Recorded in the Previous 15 seconds', fontsize=16)
    plt.title("Decay curve of P-31 by $β^+$ emission", fontsize=16)

The model I am using (which admittedly I'm not confident on my programming here) is:

def expdecay(x, t, A): 
     return A*exp(-x/t)

    decayresult = emodel.fit(amps, x=ts, t=150, A=140)
    ax.plot(ts, decayresult.best_fit, 'r-', label='best fit')
    
    print(decayresult.fit_report())

But I don't think this takes account for the uncertainties, just, plots them on the graph. I would like it to fit the exponential decay curve having taken account for the uncertainties and return the half life (t in this case) and reduced chi^2 with their respective uncertainties.

Aiming for something like the picture below, but accounting for the uncertainties in the fitting:

Aiming for this but accounting for the uncertainties in the fit

Using the weight=1/sqrt(amps) suggestion, and the full data set, I get:

Weighted fit

Which is, I imagine, the best fit (reduced chi^s of 3.89) from this data possible. I was hoping it'd give me t=150s, but hey, that one is on the experiment. Thanks for the help all.

Stef
  • 28,728
  • 2
  • 24
  • 52
Epideme
  • 209
  • 3
  • 11
  • 1
    Scipy curve_fit has the [parameter sigma](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) for weighted fitting. You should also include information which library you used for your current fit - maybe there is an option you are not aware of. – Mr. T Feb 11 '21 at 12:39
  • 1
    you can specify `weights = 1/yerr` where `yerr` is a numpy array with uncertainties, but if you use `yerr = sqrt(amps)` as in your example I guess you won't get much of an improvement – Stef Feb 11 '21 at 12:43
  • I'd seen the weight = 1/yerr thing around but didn't really understand why this might be used or why it would help much. Wouldn't setting it to 1/yerr just mean 'ignoring' big uncertainty points more than small uncertainty points rather than really including it in the fit – Epideme Feb 11 '21 at 14:36
  • I've got numpy, matplotlib, curve_fit from scipy, and lmfit imported if that helps? But I'm happy to use any others if it is a better way to do it – Epideme Feb 11 '21 at 14:38
  • 1
    The use of `lmfit` is not self-explanatory. I do not know this library well enough but maybe @m-newville will help you. However, as Stef said, looking at your posted image, you cannot expect too much from a weighted fit. The original data are quite noisy, and the fit looks reasonably good. – Mr. T Feb 11 '21 at 14:47
  • Yeh, the initial data certainly isn't great, but trying to get the maximum out of it, because there's 17 other data sets the code can be reused on, and to take more after that I have to bid for time on an accelerator. lmfit is a non-linear least-squares minimization fitting. Is what I am describing just done via a 'weighted fit' then? or is there some whole other function specifically for fitting with uncertainties? – Epideme Feb 11 '21 at 14:53
  • 1
    I am not aware of methods other than weighted fits to take uncertainty into consideration during the fitting procedure (that does not mean there are none). – Mr. T Feb 11 '21 at 14:56
  • Okay, thank you for the help then. Is it worth typing this as an answer so I can select it as the accepted answer? – Epideme Feb 11 '21 at 14:57
  • 1
    Nah, too vague for my liking. Maybe somebody else has something more substantial to contribute. – Mr. T Feb 11 '21 at 14:59
  • 1
    yes, @Stef is correct: use `weights=1.0/yerr` is exactly correct. And, yes that does mean giving more importance to data points with very low uncertainty and less importance to data points with high uncertainty. That is exactly what you want - you are more certain about the values with low uncertainty. – M Newville Feb 11 '21 at 15:33

1 Answers1

2

You can specify weights with the weights parameter. To give more weight to values with small uncertainties use for instance 1/uncertainty.
The problem with your uncertainties in the example is, however, that they directly depend on the values of the amplitude (uncertainty=np.sqrt(amps)). If you use such kinds of uncertainties they will just shift your fitted curve downwards. So this approach only makes sense if your uncertainties are real uncertainties obtained from some kind of measurement.

Example:

import matplotlib.pyplot as plt
import numpy as np
import lmfit

ts = np.array([ 15,  32,  51, 106, 123, 142, 160, 177, 196, 213, 232, 249, 269, 286, 323, 340, 359, 375, 394, 466, 484, 520, 539, 645, 681])
amps = np.array([78, 64, 64, 42, 42, 15, 34, 29, 34, 31, 31, 22,  5,  6,  8,  4, 11, 14, 14,  1,  2, 10,  4,  3,  1])
emodel = lmfit.Model(lambda x,t,A: A*np.exp(-x/t))

plt.errorbar(ts, amps, xerr=2, yerr=np.sqrt(amps), fmt="ko-", capsize = 5)
plt.plot(ts, emodel.fit(amps, x=ts, t=150, A=140).best_fit, 'r-', label='best fit')
plt.plot(ts, emodel.fit(amps, x=ts, weights=1/np.sqrt(amps), t=150, A=140).best_fit, 'r--', label='weighted best fit (1/err)')
plt.plot(ts, emodel.fit(amps, x=ts, weights=1/amps, t=150, A=140).best_fit, 'r:', label='weighted best fit (1/err²)')
plt.legend()

enter image description here

Stef
  • 28,728
  • 2
  • 24
  • 52
  • Here the data points are measurements gained from an experiment. I was using yerr = sqrt(amps) since it's a radioactive decay experiment and the uncertainty in a decay measurement is defined as the square root of it's measurement because it's a random process. So it looks like the way to go is the weights method as you suggest. What determines if we should use 1/err, 1/err^2, or any other function that gets larger as the err gets smaller? – Epideme Feb 11 '21 at 16:48
  • Out of interest, how'd you manage to get the data array out of just the graph like that... that's pretty impressive – Epideme Feb 11 '21 at 17:04
  • 1
    I'm not a statistician but I think this can't be decided a priori what are optimal weights. So a common approach is to first fit without weights and then determine the weights dependinging on the deviates delta = abs(y - y(fitted)), e.g. by using w = 1/max(delta, eps) where eps is a lower bound. See for instance https://www.springer.com/de/book/9783658114558, chapter 3.3 or https://www.itl.nist.gov/div898/handbook/pmd/section4/pmd452.htm#uw. I think you should better ask this question on https://stats.stackexchange.com. – Stef Feb 11 '21 at 17:07
  • 1
    I got the data using my favorite online digitizer https://automeris.io/WebPlotDigitizer/ – Stef Feb 11 '21 at 17:07
  • To stats SE it is then! Thank you, it's been helpful. I've posted the results at the end of the question for any looking at this in the future. Think that covers everything on this question. The online digitizer is pretty cool, colour me impressed with that. – Epideme Feb 11 '21 at 17:10
  • One other thing I meant to ask; in this answer the model is changed to `emodel = lmfit.Model(lambda x,t,A: A*np.exp(-x/t))` What is the difference of the change and what role does the lambda take on here? – Epideme Feb 12 '21 at 09:29
  • 1
    There's no difference, it's just for brevity to make the example more compact. Using a regular function like you did is absolutely fine. – Stef Feb 12 '21 at 09:43
  • 1
    @Stef well done! I might suggest printing out the values (and uncertainties!) of the `t` and `amp` parameters too. – M Newville Feb 12 '21 at 12:45
  • Ha, funny you should say that - I've just been trying to do that. Thought it might be nice to have the weighted and unweighted line on the graph and a text box with t, A, chi^2 for each in it. But I keep printing the same text in both boxes at the moment somehow. – Epideme Feb 12 '21 at 15:47