Phys 338:  Nonlinear Model Fits

The Importance of Weights

The default behavior of LinearModelFit and NonlinearModelFit is to assume equal uncertainties for all the data points.  But sometimes that isn’t a valid assumption.  Other times, it’s a reasonable assumption for the original data, but some transformation of the data violates that assumption.

In this notebook, we will look at fitting data for the decay of muons from the muon lifetime experiment you may have done in Phys 130.

Clear["Global`*"]
DateString[]
SetDirectory[NotebookDirectory[]]

Theory

The number of decays measured as a function of time is expected to follow an exponential distribution with lifetime τ:

counts[t_, a_, τ_] := a Exp[-t/τ]

This expression assumes there is no background radiation.  For this exercise, an estimate (3.5) of the background counts has already been removed.

Fits

Reading in the Data

The data file contains some comment lines at the top, so just pick the valid data entry lines.

data = Select[Import["netcounts.txt", "Table"], Length[#] == 2 && VectorQ[#, NumericQ] &];

dataplot = ListPlot[data, PlotStyle -> Red, AxesLabel -> {"t (ns)", "Counts"}, AxesStyle -> Medium, ImageSize -> Large]

Unweighted Exponential Fit

This looks more or less like an exponential decay.  Try a simple nonlinear fit.  You need to give plausible guesses for the parameters

fit1 = NonlinearModelFit[data, counts[t, a, τ], {{a, 13000}, {τ, 2000}}, t]
fit1["ParameterConfidenceIntervalTable"]
MatrixForm[fit1["CorrelationMatrix"]]

Note the strong anti-correlation here.  Since we don’t have a very long run (only about 5 τ), increases in ‘a’ can partially be offset by decreases in ‘τ’.

Manipulate[Show[{dataplot, Plot[counts[t, a, τ], {t, 1500, 12000}]}], {{a, 1300}, 500, 3000}, {{τ, 2200}, 1500, 3000}]

Show[{dataplot, Plot[fit1[t], {t, 0, 13000}, PlotRange -> All]} ]

Transformation to an unweighted Linear fit.

On a semilog plot, this should be a straight line.  Take the natural log of the counts.  Note the scatter at the end -- that is significantly affected by our estimate of the background.

lndata = Table[{data[[i, 1]], Log[data[[i, 2]]]}, {i, 1, Length[data]}];
lndataplot =  ListPlot[lndata, PlotStyle -> Red, AxesLabel -> {"t (ns)", "Ln(Counts)"}, AxesStyle -> Medium, ImageSize -> Large]

This may be fit by a straight line.  We could use LinearModelFit and extract the constants, but here we will use NonlinearModelFit so we can use the same command structures for both fits.

lncounts[t_, a_, τ_] := a - t/τ

lnfit1 = NonlinearModelFit[lndata, lncounts[t, a, τ], 
{{a, Log[13000]}, {τ, 2000}}, t]
lnfit1["ParameterConfidenceIntervalTable"]
MatrixForm[lnfit1["CorrelationMatrix"]]

We can show the two fits on the same graph by either taking the exponential of the log fit, or by taking the log of the original exponential fit.

Show[{
lndataplot, 
Plot[{Log[fit1[t]], lnfit1[t]}, {t, 0, 13000}, PlotRange -> All, PlotLegends -> {"Exponential", "Linear"}]
} ]

Show[{
dataplot, 
Plot[{fit1[t], Exp[lnfit1[t]]}, {t, 0, 13000}, PlotRange -> All, PlotLegends -> {"Exponential", "Linear"}]
} ]

The two fits are noticeably different on the log scale.  The Exponential fit matches the early data better, while the linear fit matches the later data better.  The difference has to do with weighting.

Uncertainties and Weighted Fits

Needs["ErrorBarPlots`"]

Unweighted Fits Assuming Equal Uncertainties in the Raw data

Add a third column for uncertainties.  For an initial pass, assume that all the uncertainties are the same.  Arbitrarily pick a number, e.g. 3.5.  (The estimated background in this case was 3.5.)  Add a third column to the data with that assumed uncertainty.

unc = 3.5 (* Arbitrary guess at uncertainty in count *);dataU = Table[{data[[i, 1]], data[[i, 2]], unc}, {i, 1, Length[data]}];

For the semilog data, a constant uncertainty of +- 3 in “n” counts can be changed into an approximate uncertainty of(Log[n+3] - Log[n]).

lndataU = Table[{data[[i, 1]], Log[data[[i, 2]]], Log[data[[i,2]]+unc] - Log[data[[i, 2]]]}, {i, 1, Length[data]}];

dataUplot = ErrorListPlot[dataU, PlotStyle -> Red, AxesLabel -> {"t (ns)", "Counts"}, AxesStyle -> Medium, ImageSize -> Large];

lndataUplot = ErrorListPlot[lndataU, PlotStyle -> Purple, AxesLabel -> {"t (ns)", "Log[Counts]"}, AxesStyle -> Medium, ImageSize -> Large];

Show[{
dataUplot, 
Plot[{fit1[t], Exp[lnfit1[t]]}, {t, 0, 13000}, PlotRange -> All, PlotLegends -> {"Exponential", "Linear"}]
} ]

The uncertainties are not visible on this linear scale.  The exponential fit assumed these equal uncertainties.

Show[{
lndataUplot, 
Plot[{Log[fit1[t]], lnfit1[t]}, {t, 0, 13000}, PlotRange -> All, PlotLegends -> {"Exponential", "Linear"}]
} ]

This plot shows that equal uncertainties on a linear scale (+-3.5 counts) translate into very unequal uncertainties on the semilog scale.  The linear fit assumed equal uncertainties on this graph, and so gave those long-time small counts equal importance.  The curve thus tries to go through some of those points at the expense of the better-constrained high count points at the beginning of the run.

Unequal Uncertainties in the Raw data

For nuclear statistics with large numbers of counts, the uncertainty in “n” counts is given by Sqrt[n].  Although we don’t always have “large” numbers of counts (especially at late times) we will use Sqrt[n] anyway.  Make a separate table of the uncertainties for later use.

unc2 = Table[Sqrt[data[[i, 2]]], {i, Length[data]}];
dataU2 = Table[{data[[i, 1]], data[[i, 2]], unc2[[i]]}, {i, 1, Length[data]}];
lndataU2 = Table[{data[[i, 1]], Log[data[[i, 2]]], Log[data[[i,2]]+unc2[[i]]] - Log[data[[i, 2]]]}, {i, 1, Length[data]}];

dataU2plot = ErrorListPlot[dataU2, PlotStyle -> Red, AxesLabel -> {"t (ns)", "Counts"}, AxesStyle -> Medium, ImageSize -> Large];

lndataU2plot = ErrorListPlot[lndataU2, PlotStyle -> Purple, AxesLabel -> {"t (ns)", "Log[Counts]"}, AxesStyle -> Medium, ImageSize -> Large];

Weighted Fits Assuming Unequal Uncertainties in the Raw data

NonlinearModelFit has a way to incorporate unequal uncertainties.  Each data point is weighted by a weight (Subscript[w, i]) = 1/((Subscript[σ, i])^2) .  Thus terms with small uncertainties are more important than terms with large uncertainties.  For NonlinearModelFit, we generate a vector of weights:

weights = Table[1/((unc2[[i]])^2), {i, 1, Length[data]}];

fitw = NonlinearModelFit[data, counts[t, a, τ], {{a, 13000}, {τ, 2000}}, t, Weights -> weights]
fitw["ParameterConfidenceIntervalTable"]
MatrixForm[fitw["CorrelationMatrix"]]

Show[{
dataU2plot, 
Plot[{fitw[t], fit1[t], Exp[lnfit1[t]]}, {t, 0, 13000}, PlotRange -> All, PlotLegends -> {"Weighted Exponential", "Exponential", "Linear"}]
} ]

The unweighted exponential fit assumes equal uncertainties on this linear scale, but the high count data early on actually has larger uncertainties.

Show[{
lndataU2plot, 
Plot[{Log[fitw[t]], Log[fit1[t]], lnfit1[t]}, {t, 0, 13000}, PlotRange -> All, PlotLegends -> {"Weighted Exponential", "Exponential", "Linear"}]
} ]

Comparison of estimated lifetimes:

TableForm[{
{"Weighted exponential", fitw["BestFitParameters"]},
{"Unweighted exponential", fit1["BestFitParameters"]},
{"Unweighted linear", lnfit1["BestFitParameters"]}
}]

For comparison, here is the standard value:

UnitConvert[ParticleData["Muon", "Lifetime"], "Nanoseconds"]

Summary

By default, both LinearModelFit and NonlinearModelFit assume equal uncertainties in the y-values of the data as given.  This is not always a valid assumption.  If the original data have equal uncertainties but have been transformed in some way (e.g. for a log-log or semilog plot, or square root) then the uncertainties are no longer necessarily equal.  On the other hand, if the original data have unequal uncertainties, then that should be accounted for.  In either case, the quality of the fit as well as the numerical values of the parameters can be affected.  If the uncertainties are all small, this often ends up not mattering much, but if the uncertainties are large, it can become significant.