Error propagation in finding full width at half maximum

I have an equation given by:

f[x_,M_]:= 96.4529 E^(M x (-0.216331 - 38.7396 M x)) + 31.0508 E^(M x (0.306405 - 18.585 M x)) + 4.36041 E^(M x (3.95974 - 7.37814 M x)) + 2.00366 E^(M x (-1.54639 - 3.79704 M x)) + 119.8 E^(M x (-0.0235058 - 0.0245919 M x));


where $$x$$ is a variable, $$M$$ is the mean value from a measurement and $$\delta M$$ is the error in $$M$$. For example, let $$M \pm \delta M$$ be equal to $$15.584 \pm 0.045$$.

The approximate error $$\delta f$$ in $$f$$ can be given, in the simplified form assuming indipendence of $$x$$ and $$M$$, as:

$$\delta f = \sqrt{\left(\frac{\partial f}{\partial M}\right)^2(\delta M)^2 + \left(\frac{\partial f}{\partial x}\right)^2(\delta x)^2}$$

As $$x$$ is not a measured value, I take $$\delta x = 0$$ :

DfM[x_, M_] := Evaluate@ D[f[x, M], M];
Dfx[x_, M_] := Evaluate@ D[f[x, M], x];
Deltaf[x_, M_, deltaM_] := Evaluate@ Simplify[Sqrt[(DfM[x, M])^2 deltaM^2
+ (Dfx[x, M])^2 deltax^2]] /. {deltax -> 0};


This way I can find the error $$\delta f$$ in $$f$$ that is a consequence of the error $$\delta M$$ in $$M$$.

But my aim is to find the full width at half maximum (FWHM) i.e. $$FWHM \pm \delta FWHM = (x_2-x_1) \pm \sqrt{(\delta x_1)^2 + (\delta x_2)^2}$$ along with the error $$\delta FWHM$$ that propagates in doing so. Here, $$x_1$$ and $$x_2$$ are the corresponding $$x$$-axis values for the half maximum HM. The errors in $$x_1$$ and $$x_2$$ are $$\delta x_1$$ and $$\delta x_2$$, respectively.

To find the half of the maximum (i.e. HM) , I find the maxima and minima of f[x,M] as:

fMin = First[NMinimize[{f[x, 15], -0.08 < x < 0.04}, x]];
fMax = First[NMaximize[f[x, 15], x]];
HM = fMax + (fMin - fMax)/2;


There will be some error in HM i.e. $$\delta HM$$ but I am unable to figure out how this error is propagating via the NMaximize and NMinimize functions. Additionally, to find the $$x_1$$ and $$x_2$$ I need to solve the function f[x,M] for the value HM but this function does not seem to have an inverse and therefore Solve fails if I use:

Solve[f[x, 15] == HM, x]


Is there a way I could find $$FWHM$$ and its error $$\delta FWHM$$?.

EDIT:

Link to the data of which fun[x] is a nonlinear fit: Data

fun[x_] := 96.4529 E^((-0.216331 - 38.7396 x) x) + 31.0508 E^((0.306405 - 18.585 x) x) +  4.36041 E^((3.95974 - 7.37814 x) x) +  2.00366 E^((-1.54639 - 3.79704 x) x) +  119.8 E^((-0.0235058 - 0.0245919 x) x);


the function f[x_,M_] given above is derived using fun[x_] as:

f[x_,M_] := fun[x*M];

• I think you're almost there: Replace Solve[f[x, 15] == HM, x] with x1 = x /. FindRoot[f[x, 15] == HM, {x, -0.1}] and x2 = x /. FindRoot[f[x, 15] == HM, {x, -x1}]. (I'll work on the error part later if no one else does it sooner.) – JimB Oct 31 '18 at 19:43
• Why the apparent arbitrary restriction for finding the minimum? -0.08 < x < 0.04 Otherwise with minimum would seem to be zero. – JimB Oct 31 '18 at 19:59
• The function f[x,M] is a fitted function from a data that is collected in the domain -0.08 < x < 0.04 and has the background amplitude closer to fMin rather than a 0 amplitude. – jsid Oct 31 '18 at 20:14
• Because this is all based on data, you might be better off using a bootstrap approach rather than a propagation of error (Delta method to us statisticians) approach. Is it possible to make the data available? – JimB Oct 31 '18 at 21:46
• You might want: x /. Solve[f[x, 15] == HM, x, Reals]. – Daniel Lichtblau Oct 31 '18 at 23:25

Here is a bootstrap approach to estimate a standard error associated with the estimate of $$x_2-x_1$$.

data = {{-1.5, 117.955}, {-1.452, 118.277}, {-1.404, 119.012}, {-1.356, 119.277},
{-1.308, 120.204}, {-1.26, 120.545}, {-1.212, 120.866}, {-1.164, 120.712},
{-1.116, 120.37}, {-1.068, 120.523}, {-1.02, 120.848}, {-0.972, 120.798},
{-0.924, 120.98}, {-0.876, 121.441}, {-0.828, 121.496}, {-0.78, 121.163},
{-0.732, 120.515}, {-0.684, 121.504}, {-0.636, 121.139}, {-0.588, 119.96},
{-0.54, 119.393}, {-0.492, 120.752}, {-0.444, 122.253}, {-0.396, 123.639},
{-0.348, 126.277}, {-0.3, 131.765}, {-0.252, 139.704}, {-0.204, 154.437},
{-0.156, 180.332}, {-0.108, 214.055}, {-0.06, 242.769}, {-0.012, 255.787},
{0.036, 249.115}, {0.084, 224.131}, {0.132, 194.485}, {0.18, 171.191},
{0.228, 155.534}, {0.276, 143.708}, {0.324, 134.348}, {0.372, 129.606},
{0.42, 127.123}, {0.468, 124.993}, {0.516, 122.399}, {0.564, 120.834},
{0.612, 120.375}, {0.66, 119.378}, {0.708, 118.306}, {0.756, 116.853},
{0.804, 115.839}};


Now the function to be fit:

f[x_, a1_, b1_, c1_, a2_, b2_, c2_, a3_, b3_, c3_, a4_, b4_, c4_] :=
a1 E^(x (b1 + c1 x)) +
a2 E^(x (b2 + c2 x)) +
a3 E^(x (b3 + c3 x)) +
a4 E^(x (b4 + c4 x));

(* Initial parameter values *)
inits = {
{a1, 55}, {b1, 1}, {c1, -20},
{a2, 20}, {b2, 2}, {c2, -1},
{a3, 12}, {b3, -1.5}, {c3, -53},
{a4, 120}, {b4, -0.5}, {c4 , -0.2}};

findx1x2[data_] :=
Module[{nlm, fmax, fmin, HM, x1, x2, xfmax, xfmin, x2Init},
nlm = NonlinearModelFit[data, f[x, a1 , b1 , c1 , a2 , b2 , c2 , a3 , b3, c3, a4 , b4, c4],
inits, x, MaxIterations -> 1000];
fmax = FindMaximum[nlm[x], {x, Select[data, #[[2]] == Max[data[[All, 2]]] &][[1, 1]]}];
fmin = FindMinimum[{nlm[x], -1.5 < x < 0.8}, {x, Select[data, #[[2]] == Min[data[[All, 2]]] &][[1, 1]]}];
xfmax = x /. fmax[[2]];
xfmin = x /. fmin[[2]];
HM = (fmin[[1]] + fmax[[1]])/2;
x1 = x /. FindRoot[nlm[x] == HM, {x, (x /. fmax[[2]])}];
x2Init = If[x1 < xfmax, xfmax + (xfmax - x1), xfmax + (xfmax - x1)];
x2 = x /. FindRoot[nlm[x] == HM, {x, x2Init}];
{Min[x1, x2], Max[x1, x2], HM, fmin[[1]], fmax[[1]]}]


The bootstrap process (random selection of residuals):

(* Get predicted responses and fit residuals *)
nlm = NonlinearModelFit[data,
f[x, a1 , b1 , c1 , a2 , b2 , c2 , a3 , b3, c3, a4 , b4, c4],
inits, x, MaxIterations -> 1000];
predictedResponse = nlm["PredictedResponse"];
fitResiduals = nlm["FitResiduals"];

(* Array to hold bootstrap estimates *)
nboot = 1000;
x1x2 = ConstantArray[{0, 0, 0, 0, 0}, nboot];

(* Perform bootstraps *)
n = Length[data];
SeedRandom[12345];
Do[
boot = Transpose[{data[[All, 1]], predictedResponse + RandomChoice[fitResiduals, n]}];
x1x2[[iboot]] = findx1x2[boot],
{iboot, nboot}]

(* Summarize results *)
mle = findx1x2[data]
(* Estimate of x2 - x1 *)
estimate = mle[[2]] - mle[[1]]
(* 0.29482173864942085 *)
(* Standard error of x2 - x1 *)
stderr = StandardDeviation[x1x2[[All, 2]] - x1x2[[All, 1]]]
(* 0.0014648255351604935 *)

(* Histogram of bootstrap estimates *)
Histogram[x1x2[[All, 2]] - x1x2[[All, 1]]]


Edit

If the quantity of interest is $$(x_2-x_1)/M$$, then applying the Propagation of Error (Delta) method results in an estimate of the variance as

$$\frac{\sigma_M^2 (x_2-x_1)^2}{M^4}+\frac{\sigma_{x_2-x_1}^2}{M^2}$$

assuming that the estimate of $$x_2-x_1$$ is statistically independent from the estimate of $$M$$. So we have

results = {x1 -> mle[[1]], x2 -> mle[[2]], σx1x2 -> stderr, M -> 15.584, σM -> 0.045}
estimate = (-x1 + x2)/M /. results
(* 0.0189182 *)
seEstimate = (σx1x2^2/M^2 + σM^2 (x2 - x1)^2/M^4)^0.5 /. results
(* 0.000108717 *)

• This certainly looks like a better approach to me however, how do I introduce the factor $M$ here?. Do I divide the $x$-axis values of the data by $M=15.584$ and replace the $x$ variable in your nonlinear model function f with {x -> x*M}?. – jsid Nov 1 '18 at 18:08
• Injecting $M$ into the mix would only affect the estimated coefficients but NOT any predictions of the dependent variable or inferences about $x_1$, $x_2$, or $x_2-x_1$. Or maybe I'm not understanding what you want. Is it $M(x_2-x_1)$ or $(x_2-x_1)/M$ the quantity of interest? – JimB Nov 1 '18 at 18:30
• What I did was that I first fit the data and then change the fitted equation fun[x_] by introducing $M$ as f[x_,M_] := fun[x*M]. The factor $M$ basically reduces the width of fitted function while keeping the same amplitude. Since your answer works directly on the NonlinearModelFit, consider the following manipulation: data = data/. {a_Real, b_Real} :> {a/15.584, b}. The fitted parameters would change accordingly and $(x_2-x_1)$ will become small. The $FWHM$ should then be $(x_2-x_1)/M$ I guess?. – jsid Nov 1 '18 at 19:59

This is an extended comment rather than an answer.

While I still think that a bootstrap approach would be more accurate than using the Propagation of Error/Delta Method, there are two issues you might want to consider first.

1. The parameter correlation matrix for your current model suggests that using 5 kernels is way overfitting (with $$a_i e^{M x(b_i+c_i M x)}$$ being a kernel). Fitting a model with just 4 kernels has a much smaller AIC (Akaike Information Criterion) and fits almost perfectly with no errors.

2. If you are after predictions or making inferences about $$x_1$$ and $$x_2$$, $$M$$ is completely irrelevant because $$M$$ and $$x$$ always appear together as a product with an associated coefficient.

More specifically the kernel $$a e^{M x(b+c M x)}$$ results in the same prediction as $$a e^{x(B+C x)}$$. (If $$M$$ appeared away from $$x$$, then that would be a different story.)