1
$\begingroup$

Problem:

I'm trying to plot an exponential decay curve on a set of integer data, and print the various parameters and their uncertainties on the plot itself. Having adapted code from another question, this worked successfully with non-integer values with uncertainties in the data set. Here the data set is integer only and does not contain uncertainties and my attempt at adapting it has been unsuccessful.

The target output is something like this: Non poisson fit

Code: (contains data)

dataHist5 = {{18, 74}, {36, 64}, {54, 62}, {72, 54}, {90, 47}, {108, 
    39}, {126, 40}, {144, 35}, {162, 34}, {180, 29}, {198, 34}, {216, 
    30}, {234, 31}, {252, 22}, {270, 14}, {288, 14}, {306, 13}, {324, 
    25}, {342, 18}, {360, 11}, {378, 13}, {396, 16}, {414, 13}, {432, 
    12}, {450, 10}, {468, 12}, {486, 11}, {504, 13}, {522, 9}, {540, 
    8}, {558, 7}, {576, 5}, {594, 5}, {612, 4}, {630, 5}, {648, 
    2}, {666, 1}, {684, 3}, {702, 2}, {720, 2}, {738, 1}, {756, 
    1}, {774, 1}, {792, 0}, {810, 0}, {828, 1}, {846, 2}, {864, 
    0}, {882, 1}, {900, 1}};
fitData = dataHist5;

Clear[A, k]
uncertainties = sqrt[dataHist5[[2]]];
fit = NonlinearModelFit[fitData, A Exp[-k t], {A, k}, t, 
  Weights -> 1/dataHist5[[2]]]

{A, k} = {A, k} /. fit["BestFitParameters"];
{\[Sigma]A, \[Sigma]k} = fit["ParameterErrors"];
hLife = Log[2]/Around[k, \[Sigma]k];
halfLife = hLife[[1]];
seA = Around[A, \[Sigma]A];
sehalfLife = hLife[[2]];

Show[Plot[fit[x], {x, 0, 900}, PlotRange -> All, 
  PlotTheme -> "Detailed", PlotStyle -> Red, Axes -> False, 
  Frame -> {{True, False}, {True, False}}, 
  FrameLabel -> {"Time /s", 
    "Counts Recorded in the Previous 15 seconds"}, 
  ImageSize -> Large], ListPlot[dataHist5, ImageSize -> Large], 
 Graphics[Inset[
   Framed[Column[{Style["Run 0", Bold], 
      Row[{"Data Points = ", Length[dataHist5], "/50"}], 
      Row[{Subscript[t, Style["1/2", FontSize -> 10]], " = ", 
        PlusMinus[NumberForm[halfLife, 4], 
         NumberForm[sehalfLife, 3]]}], 
      Row[{"A = " PlusMinus[Round[A], Round[\[Sigma]A]]}],
      Row[{"\[Lambda] = " PlusMinus[NumberForm[k, 3], 
          NumberForm[\[Sigma]k, 2]]}], 
      Row[{Superscript[\[Chi], 2], "= ", 
        NumberForm[fit["ANOVATableSumsOfSquares"][[2]], 4]}], 
      Row[{"Reduced " Superscript[\[Chi], 2], "= ", 
        NumberForm[fit["ANOVATableMeanSquares"][[2]], 3]}]}], 
    Background -> White, RoundingRadius -> 5], {Right, Top}, 
   Scaled[{1.1, 1.2}]]], 
 PlotLabel -> 
  Style["Decay Curve of Phosphorus-30 by \[Beta]+ Emission", Bold]]

Current Output:

poisson fit

Among many errors are 'NonlinearModelFit::wtsln: The number of weights 2 specified by Weights->{0.0277778,0.015625} is not the same as the number of data points 50.' and 'Set::shape: Lists {A,k} and {A,k}/. NonlinearModelFit[<<1>>][BestFitParameters] are not the same shape.'(though there are other errors as well)

Attempted Solution (Following JimB's answer):

The ParameterTable elements are seemingly being pulled, but not the DevianceTable ones, and I'm not super sure I've worked out the uncertainty in A and t correctly. Also the axes labels aren't appearing.

glm = GeneralizedLinearModelFit[dataHist5, t, t, 
  ExponentialFamily -> "Poisson"]

halfLife = (Log[E, 2]/glm["ParameterTableEntries"][[2, 1]])
sehalfLife = (((glm["ParameterTableEntries"][[2, 2]])/(glm[
    "ParameterTableEntries"][[2, 1]])))*halfLife
k = glm["ParameterTableEntries"][[2, 1]]
\[Sigma]k = glm["ParameterTableEntries"][[2, 2]]
edp = glm["DevianceTableEntries"][[4, 2]]
redp = (glm["DevianceTableEntries"][[4, 2]])/(glm[
     "DevianceTableEntries"][[3, 2]])
A = E^(glm["ParameterTableEntries"][[1, 1]])
\[Sigma]A = (A*(glm["ParameterTableEntries"][[1, 1]])/E)

Show[ListPlot[dataHist5], 
 Plot[glm[t], {t, 0, 900}, PlotRange -> All, PlotTheme -> "Detailed", 
  PlotStyle -> Red, Axes -> False, 
  Frame -> {{True, False}, {True, False}}, 
  FrameLabel -> {"Time /s", 
    "Counts Recorded in the Previous 15 seconds"}, 
  ImageSize -> Large], 
 Graphics[Inset[
   Framed[Column[{Style["Run 0", Bold], 
      Row[{"Data Points = ", Length[dataHist5], "/50"}], 
      Row[{Subscript[t, Style["1/2", FontSize -> 10]], " = ", 
        PlusMinus[NumberForm[halfLife, 4], 
         NumberForm[sehalfLife, 3]]}], 
      Row[{"A = " PlusMinus[Round[A], Round[\[Sigma]A]]}],
      Row[{"\[Lambda] = " PlusMinus[NumberForm[k, 3], 
          NumberForm[\[Sigma]k, 2]]}], 
      Row[{Superscript[\[Chi], 2], "= ", NumberForm[edp, 4]}], 
      Row[{"Reduced " Superscript[\[Chi], 2], "= ", 
        NumberForm[redp, 3]}]}], Background -> White, 
    RoundingRadius -> 5], {Right, Top}, Scaled[{1.1, 1.2}]]], 
 PlotLabel -> 
  Style["Decay Curve of Phosphorus-30 by \[Beta]+ Emission", Bold]]

enter image description here

$\endgroup$
9
  • 2
    $\begingroup$ I think maybe you meant Weights -> 1/dataHist5[[All, 2]]? Also, you have a lower case s in Sqrt’ for your uncertainties and you have [[2]]` there instead of [[All, 2]]. $\endgroup$ – MassDefect Mar 31 at 11:46
  • $\begingroup$ Yep that turns out to be it. Also I need to remove the points with a zero in it, otherwise I get a lot of 1/0 complexiinfinity type errors. Would you want to post it as an answer so it can be accepted / upvoted? $\endgroup$ – Epideme Mar 31 at 11:53
  • 1
    $\begingroup$ If you specify Weights->1/(eps+dataHist5[[All, 2]]) with eps a small number, you do not need to delete points. $\endgroup$ – Daniel Huber Mar 31 at 19:46
  • $\begingroup$ Ah, good point, thank you $\endgroup$ – Epideme Apr 1 at 17:25
  • $\begingroup$ Might throwing out points with large weight not be a good thing? Throwing out points with infinite weight even worse? Fortunately in this case you don't need to use even the approximated weights that you have. $\endgroup$ – JimB Apr 3 at 4:23
2
$\begingroup$

If the "integer data" consists of counts that might be expected to follow a Poisson distribution given the predictor value, then a Poisson regression should be considered. That can be accomplished using GeneralizedLinearModelFit.

glm = GeneralizedLinearModelFit[dataHist5, t, t, ExponentialFamily -> "Poisson"]
Show[ListPlot[dataHist5], Plot[glm[t], {t, 0, 900}]]

Poisson regression fit

One of the characteristics of a lack-of-fit is that the variability about the fitted curve is larger (or at least different) from what one would expect when there is Poisson variability. A parameter that characterizes that is the "dispersion" parameter. It should be around 1.0 if the observed variability is what is expected with a Poisson distribution.

GeneralizedLinearModelFit has an option for "EstimatedDispersion" but it appears to always have the value 1.0 so I don't know what it is doing. A more reasonable estimator is residual deviance divided by the residual degrees of freedom:

glm["DevianceTable"]

Deviance table

The estimated dispersion parameter is 42.4657/48 = 0.884702 which is a bit less than one suggesting that there is no additional variation in the data to be explained after the Poisson model is considered.

So what model is actually being fit?

$$Y|t \sim \text{Poisson}(e^{a+b t})$$

We have the coefficients of the fit:

glm["ParameterTable"]

Parameter table

So the estimate of the intercept ($a$) is found in the row labeled 1 and the estimate of the slope ($b$) is found in the row labeled t.

In short in this case there is no need to use weights.

$\endgroup$
6
  • $\begingroup$ Hi, as always thank you for the answer. I attempted to recreate, and got the graph and the tables, but didn't manage to print anything to the graph. I think mostly because I'm pulling details from the tables wrong, but not sure how. Also not totally sure on my values for the uncertainty in the intercept either. I'm trying to make a graph something like the one at the top but with the poisson distribution fitting, with the important values for the graph visible in the top corner $\endgroup$ – Epideme Apr 7 at 16:21
  • $\begingroup$ Also meant to ask here - how come the fitting doesn't produce an integer intercept, would that not be expected? Here I presume it is A +/- sigma(A) = e^(estimate,1) +/- e^(standard error,1) $\endgroup$ – Epideme Apr 8 at 9:22
  • $\begingroup$ There is an intercept: $a$. But if you wanted a model that was $Y|t\sim \text{Poisson}(\text{intercept}+e^{a+bt})$, then you wouldn't have a generalized linear model. You need to specify the exact model you want. Then one chooses the function to estimate the parameters. Specifying the model includes both the "fixed effects" ($\text{intercept}+e^{a+bt}$) and the "random effects" or error structure ($\sim \text{Poisson}()$). $\endgroup$ – JimB Apr 8 at 16:00
  • $\begingroup$ Here I am trying to find the half-life and goodness of fit, and the intercept. It would seem to me, deriving from the model Y = Poisson (e^(a+bt)), that the intercept would be e^4.33202 +/- e^0.0521607 ? $\endgroup$ – Epideme Apr 8 at 16:05
  • $\begingroup$ That's not how it works. $e^{4.33202}$ is the estimate of the Poisson mean when $t=0$. Is that a quantity of interest? If that is a quantity of interest, then an approximate 95% confidence interval for that mean would be $e^{4.33202\pm 1.96*0.0521607}$. $\endgroup$ – JimB Apr 8 at 16:19

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