4
$\begingroup$

I have some data from a radioactive decay experiment that I'm trying to fit an exponential decay curve on that will take account of the uncertainties on the data, and then ideally return the data for half-life and such (time taken for the y-value to decrease by half)

The uncertainties are entered in using 'around'. Currently my code is (using a cut down data set):

dataHist5 = {{Around[16.5, 1.5], 
   Around[77.8, 8.8]}, {Around[34.5, 1.5], 
   Around[60.5, 8.0]}, {Around[52.5, 1.5], 
   Around[63.8, 8.0]}, {Around[106.5, 1.5], 
   Around[42.4, 6.5]}, {Around[124.5, 1.5], 
   Around[41.7, 6.5]}, {Around[142.5, 1.5], 
   Around[14.6, 3.8]}, {Around[160.5, 1.5], 
   Around[33.9, 5.8]}, {Around[178.5, 1.5], 
   Around[29.4, 5.4]}, {Around[196.5, 1.5], 
   Around[33.5, 5.8]}, {Around[214.5, 1.5], 
   Around[30.9, 5.6]}, {Around[232.5, 1.5], 
   Around[31.1, 5.8]}, {Around[250.5, 1.5], 
   Around[21.5, 4.6]}, {Around[268.5, 1.5], 
   Around[4.3, 2.1]}, {Around[286.5, 1.5], 
   Around[6.4, 2.5]}, {Around[322.5, 1.5], 
   Around[7.5, 2.7]}, {Around[340.5, 1.5], 
   Around[4.5, 2.1]}, {Around[358.5, 1.5], 
   Around[11., 3.3]}, {Around[376.5, 1.5], 
   Around[14.0, 3.7]}, {Around[394.5, 1.5], 
   Around[14.0, 3.7]}, {Around[466.5, 1.5], 
   Around[0.6, 0.7]}, {Around[502.5, 1.5], 
   Around[2.2, 1.5]}, {Around[520.5, 1.5], 
   Around[9.4, 3.1]}, {Around[538.5, 1.5], 
   Around[4.1, 2.0]}, {Around[646.5, 1.5], 
   Around[2.2, 1.5]}, {Around[682.5, 1.5], Around[0.6, 0.7]}}
ListPlot[dataHist5, PlotTheme -> "Detailed"]
model = a Exp[-kt];
fit = FindFit[dataHist5, model, {a, k}, t]

which returns: decay graph

And

FindFit::fitm: Unable to solve for the fit parameters; the design matrix is nonrectangular, non-numerical, or could not be inverted.

for the fitting. However, I think even had I got it to fit, it would just fit the points rather than the uncertainties around the points as well. Anyone know of a good way to implement this and have mathematica report the fitting back? Thanks for the help

EDIT: Trying to recreate JimB's answer to get Chi^2:

data = Transpose[{dataHist5[[All, 1, 1]], dataHist5[[All, 2, 1]]}];
ListLogPlot[data]
logData = data;
logData[[All, 2]] = Log[data[[All, 2]]]/data[[All, 1]];
nlm = NonlinearModelFit[logData, loga/t - k, {loga, k}, t];
nlm["BestFitParameters"]
(*{loga\[Rule]4.47236,k\[Rule]0.00675936}*)

Show[ListPlot[data], Plot[Exp[t nlm[t]], {t, 1, 700}]]
ListPlot[Transpose[{logData[[All, 1]], nlm["FitResiduals"]}], 
 PlotRange -> All]
nlm["ANOVATable"]
Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ Exponential least-squares fits of your data assume that your error bars are symmetric (Gaussian) instead of a more realistic distribution like a gamma distribution. It would be better to take the logarithm of your data points, convert the error bars carefully to this logarithmic scale, and then do a linear fit. $\endgroup$
    – Roman
    Commented Feb 11, 2021 at 18:33
  • 1
    $\begingroup$ The uncertainty values for the response variable appear to be just the square root of the response. Why is that assumed? Certainly with count data that has a Poisson distribution an estimate of the standard deviation is the square root of the counts. But your response variable is not an integer and therefore doesn't have a Poisson distribution. Also, I think more targeted help would be available if you gave the complete dataset shown in the figure. $\endgroup$
    – JimB
    Commented Feb 12, 2021 at 4:41
  • $\begingroup$ Completed data set added to the question $\endgroup$
    – Epideme
    Commented Feb 12, 2021 at 16:47
  • $\begingroup$ Added the 'Finding X^2' section, and corrected the full data set to reflect some typos I made previously $\endgroup$
    – Epideme
    Commented Mar 3, 2021 at 9:45

3 Answers 3

Reset to default
5
$\begingroup$

Note: While I was putting this answer together @Roman commented to perform a regression after taking the log of the response variable. I've implemented that suggestion below with a modification to take care of the unequal uncertainties across the values of the response variables.

If the predictor variable has substantial uncertainty, then you should NOT perform a regular regression but rather consider error-in-variables regression. Fortunately it appears that the predictor uncertainty is much smaller than the uncertainty in the response variable (both in an absolute and relative sense) so we can avoid the error-in-variables approach.

The uncertainty in the response variable appears to be simply the square root of the response variable. I can only assume that the OP is applying what one reasonably does with a Poisson distribution (variance = mean or standard deviation = square root of the mean). But the response variable is not even a count so that shouldn't be assumed and I'm going to ignore those particular uncertainties because they are unjustified.

Using the OP's data a plot of the data on a log scale shows a negative linear trend but with increasing variability as the predictor variable gets larger. That change in variability needs to be accounted for in the model as both the functional form of the model and the form of the error needs to be modeled appropriately.

data = Transpose[{dataHist5[[All, 1, 1]], dataHist5[[All, 2, 1]]}];
ListLogPlot[data]

Plot of data

Fortunately, the proposed model ($y=a e^{-kt}$) is linear on the log scale:

$$E(\log{y})=\log{a}-k t$$

where $E$ is the expectation operator. To account for the increasing variance with increasing values of the predictor variable we consider the following model with an error structure:

$$\log{y}=\log{a}-k t+t \epsilon$$ where $\epsilon \sim N(0,\sigma^2)$. Note that we can divide both sides by $t$ to get a linear model with a constant variance:

$$\log{y}/t=\log{a}/t-k+\epsilon$$

We can now modify the dataset and run a relatively simple linear model:

logData = data;
logData[[All, 2]] = Log[data[[All, 2]]]/data[[All, 1]];
nlm = NonlinearModelFit[logData, loga/t - k , {loga, k}, t];
nlm["BestFitParameters"]
(* {loga -> 4.47236, k -> 0.00675936} *)

Show[ListPlot[data], Plot[Exp[t nlm[t]], {t, 1, 700}]]

Data and fit on log scale

Alternatively, we can get the exact same fit using the "Weights" option in either LinearModelFit or NonlinearModelFit:

logData = data;
logData[[All, 2]] = Log[data[[All, 2]]];
lm = LinearModelFit[logData, -t, t, Weights -> 1/data[[All, 1]]^2];
lm["BestFitParameters"]
(* {4.47236, 0.00675936} *)

The residuals look not totally symmetric across zero but relatively uniform across values of the predictor variable:

ListPlot[Transpose[{logData[[All, 1]], nlm["FitResiduals"]}], PlotRange -> All]

Predictors vs residuals

Obtaining standard errors and an estimate of half-life

The standard errors of the parameters are found using

nlm["ParameterTable"]

Parameter Table

There are several (equivalent) ways of obtaining an estimate of the half-life along with an associated standard error. A "black box" approach is the following:

hLife = Log[2]/Around[k, σk]
(* Around[k^(-1) Log[2], Abs[k^(-2) σk] Log[2]] *)

halfLife = hLife[[1]] /. nlm["BestFitParameters"]
(* 102.546 *)
sehalfLife = hLife[[2]] /. nlm["BestFitParameters"] /. σk -> nlm["ParameterErrors"][[2]]
(* 8.22627 *)

This corresponds to an approximate 95% confidence interval of (84,116). This approach uses (internally) the Delta Method (as we statisticians call it and the rest of the world calls it Propagation of Error).

Alternatively, a more direct approach (that obtains exactly the same results) is parameterizing the model in terms of the half-life:

logData = data;
logData[[All, 2]] = Log[data[[All, 2]]]/data[[All, 1]];
nlm = NonlinearModelFit[logData, loga/t - Log[2]/halfLife, {loga, halfLife}, t];
nlm["ParameterTable"]

Parameter Table for halfLife model formulation

Improve this answer
$\endgroup$
10
  • $\begingroup$ Thank you, this is a really helpful answer. In terms of the data, it came from a radioactive decay experiment, it's measured by a coincidence detector, with each channel corresponded to some angle of a circular detector, the counts are plotted on a histogram, and a gaussian plotted every 18seconds. The area under each gaussian corresponding to the number of counts, each gaussian being the next 18second's result. Particularly anomolous results excluded. Because it came from this applied model, we got some non-integer values. Thus the lack of poisson statistics. Using sqrt(counts) just came from $\endgroup$
    – Epideme
    Commented Feb 12, 2021 at 16:34
  • $\begingroup$ that property of radioactive decay, though admittedly that should be applied with a 'real' count and not a fitted one. How did you return the whole data set by the way? Interesting that it had some to higher significant figure than I have written down. Thank you for all the help with the regression stuff, because that's saved me a ton of sitting there going through things, and the explanation is very helpful. Is there a way to return epsilon? It also appears that both the weighted linear model and the modified nlmreturn similiarly, is there any reason to prefer one of them? $\endgroup$
    – Epideme
    Commented Feb 12, 2021 at 16:39
  • $\begingroup$ The value of k returned seems to suggest the half life value of t=99.9s, which, is much lower than expected (~150s), which is much lower than I've got plotting this in other programs, is this expected? and is there a way to easily see the uncertainty on t(1/2)? $\endgroup$
    – Epideme
    Commented Feb 12, 2021 at 16:43
  • $\begingroup$ "...is this expected?" Yogi Berra said “In theory, there is no difference between practice and theory. In practice, there is.” I'll add in how to get uncertainty (in terms of standard errors) for k and t(1/2). $\endgroup$
    – JimB
    Commented Feb 12, 2021 at 17:47
  • $\begingroup$ "weighted linear model" vs. "modified nlm". While the fits are identical (not just similar) I prefer the "modified nlm" as it makes the error structure explicit. To say "I'll just use weights" is too vague and many times indicates that the user doesn't understand what the weights are doing. $\endgroup$
    – JimB
    Commented Feb 12, 2021 at 17:50
4
$\begingroup$

Try this:

dataHist5 = {{Around[19.5, 1.5], 
    Around[77.8, 8.8]}, {Around[37.5, 1.5], 
    Around[63.5, 8.0]}, {Around[55.5, 1.5], 
    Around[63.8, 8.0]}, {Around[109.5, 1.5], 
    Around[42.4, 6.5]}, {Around[127.5, 1.5], 
    Around[41.7, 6.5]}, {Around[145.5, 1.5], Around[14.6, 3.8]}};
dataHist5A = dataHist5 /. Around[x_, y_] -> x

Then you fit the dataHist5A, rather than dataHist5. Here please take care that Mma understands kt as a single variable. The expression you want to be in the exponent is k*t:

model = a Exp[-k*t];
fit = FindFit[dataHist5A, model, {a, k}, t]

(*  {a -> 91.5711, k -> 0.00798938}  *)

Now, let us show the fit on the background of the data. This time we will use already dataHist5:

Show[{
  ListPlot[dataHist5],
  Plot[model /. fit, {t, 0, 700}, PlotStyle -> Red]
  }]

enter image description here

Have fun!

Improve this answer
$\endgroup$
5
  • $\begingroup$ Oh excellent, thank you! It was that Hist5A line including the around that I was missing. It's odd that I needed k*t since i got kt from here: reference.wolfram.com/language/ref/FindFit.html Do you know how I would return half-life from here (the time it takes to deplete the y-value by half? $\endgroup$
    – Epideme
    Commented Feb 11, 2021 at 16:42
  • $\begingroup$ Your use of kt doesn't work because there is a space between the k and the t (which Mathematica interprets as multiplication). Alexei just explicitly included the multiplication sign. Half life for a first-order reaction is Ln[2]/k. $\endgroup$
    – bobthechemist
    Commented Feb 11, 2021 at 17:00
  • $\begingroup$ Is there anyway to return the uncertainties on a and k? (and by extension the half-life?) $\endgroup$
    – Epideme
    Commented Feb 12, 2021 at 16:45
  • $\begingroup$ Do you mean the errors of the fitting itself? $\endgroup$
    – Alexei Boulbitch
    Commented Feb 12, 2021 at 19:27
  • $\begingroup$ Yeh. Ideally I need to output a figure with it's uncertainties for a, k, and t. $\endgroup$
    – Epideme
    Commented Mar 1, 2021 at 11:05
1
$\begingroup$

I wrote a new ResourceFunction to automatically handle Around values in LinearModelFit and NonlinearModelFit and convert them to weights. It handles error bars in both the x and y directions. Just wrap ResourceFunction["FitWithErrors"] around the whole fit expression:

fit = ResourceFunction["FitWithErrors"] @ NonlinearModelFit[dataHist5, model, {a, k}, t];

fit // Normal

Show[
 Plot[fit[t], {t, 0, 700}],
 ListPlot[dataHist5, PlotTheme -> "Detailed"]
 ]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.