How do I fit an exponential decay curve which accounts for uncertainties?

Question

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:

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

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.

Answer 1

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()