# Uncertainty of the area of a Gaussian curve atop a linear background

I have some data from a counting-based spectroscopy experiment. Each data point is an (Energy, Rate) pair. One such data set looks like this:

I choose to fit this data to a Gaussian curve plus a background with a linear slope.

$$I(E) = A e^{-\frac{1}{2}\frac{(E-E_0)^2}{{\sigma}^2}} + m(E-E_0) + b$$

I use Mathematica to perform the fit, using the variance of each data point as the inverse weight. The result of fit is

The uncertainty in the fit parameters are the square root of the diagonal elements of covariance matrix.

Now, the area of the signal peak is $$Area = A\sigma\sqrt{2\pi}$$, and the rate of signal events is $$R = \frac{Area}{W}$$, where $$W$$ is the uniform energy bin width. For the data shown here ($$W = 0.399$$), this produces $$R = 8.103 Hz$$.

But what is the uncertainty of this rate? I could use simple propagation of errors and I get $$\Delta_R = R\sqrt{(\frac{\Delta_A}{A})^2+(\frac{\Delta_\sigma}{\sigma})^2} = 0.848$$

But this doesn't take into account the uncertainty in the background! How do I account for the uncertainty in the background when calculating the area of the Gaussian peak? Or do the variances of the fit parameters for the area - namely $$\Delta_A$$ and $$\Delta_\sigma$$ - already contain contributions owing to the uncertainties in the background parameters?

Many thanks!

• $\Delta_A$ already accounts for the uncertainty in the background. However, your propagation of errors is not quite correct, because it looks likely there is nonzero covariance between $\Delta_A$ and $\Delta_\sigma$ which you ought to incorporate. Another potential issue is how you define or compute the "variance of each data point." Yet another issue concerns how you chose this window of the spectrum for calculation: the estimates can be sensitive to the choice and even might have an inherent bias (interference from tiny but real peaks near the edges).
– whuber
Commented Jul 5, 2023 at 19:40
• Also, I have often wondered whether measurements of photons near the bin edges might induce some negative correlation in the errors, because each photon assigned to the next higher bin is one less photon assigned to the lower bin. Thus, it might be worth doing a quick check of the serial correlation of the residuals to see whether the lag-1 autocorrelation coefficient is significantly negative and, if so, potentially using that to improve the model and the estimate of uncertainty (it should decrease a little).
– whuber
Commented Jul 5, 2023 at 19:43
• @whuber, thank you for the response. The variance in each data point is derived from counting statistics and collection time. The background need not be a pure sloped line; this is a gamma spectrum from Europium-152, and there are >200 gamma lines in the 0 - 2,000 keV range. I'm looking into how to propagate uncertainty when there are off-diagonal elements in the covariance matrix. Thanks for the tip! Commented Jul 5, 2023 at 20:12
• That's the right way to estimate the variances. For low-level counts (say, all less than 100), it would be more rigorous to do the fitting using a Poisson regression based on the counts. That wouldn't require the variances to be specified at all. It doesn't matter what functional form you use for background: regardless what it might be (a spline, a quadratic, an exponential, etc.) only the estimates $\hat A$ and $\hat\sigma$ matter for assessing the uncertainty in the area to second order. For any more detail you would need to analyze the log likelihood in the neighborhood of the estimate.
– whuber
Commented Jul 5, 2023 at 20:46

I believe I've figured this out. There are a few ways to calculate the uncertainty.

NOTE: I've discovered that the results presented here are incorrect because this work deals in "count rates" instead of raw counts. I believe the methods presented are correct, but the conclusions are not. I will update this when I have the time.

## Raw fit method

The raw fit method is what I discuss in the original post. These means 1) fitting a Gaussian with a linear background, 2) calculating the area from $$R=\frac{A\sigma\sqrt{2\pi}}{W}$$, and 3) calculating the uncertainty in $$R$$ by propagating the uncertainty using the variances in the fitted parameters $$A$$ and $$\sigma$$.

The proper way the propagate the uncertainty is $$\Delta_R = R\sqrt{(\frac{\Delta_A}{A})^2 + (\frac{\Delta_\sigma}{\sigma})^2 + 2\frac{\Delta_{A\sigma}^2}{A\sigma}}$$

Notice this is a correction to the original question, because now I've included the covariance between $$A$$ and $$\sigma$$, i.e. $$\Delta_{A\sigma}$$.

Using this approach, I get $$R = 8.103 \pm 1.049$$.

## Gillmore method

Gillmore (Practical Gamma Ray Spectroscopy, pg 110, 2008) discusses estimating the background in an experiment such as mine. He comes up with an uncertainty $$\sigma = \sqrt{A + B(1+\frac{n}{2m})}$$

Here, $$B$$ is the number of background counts - determined by making a linear interpolation of the average value of events over $$m$$ channels just outside of a 'signal range' of width $$n$$ - and $$A$$ is the TOTAL number of counts in the 'signal range' MINUS $$B$$.

Gillmore presents this as and $$B$$ is simply the second term on the right hand side of this expression.

This approach suffers from its presence of free parameters, namely 1) the size of the 'signal range' $$n$$, and 2) the range over which to average the background, $$m$$. For example, here is a set of rather large ranges for the signal and background regions

In this example, the background range (red) is larger than needed, and the signal range (dark green) clearly goes outside of the range of what is reasonable signal.

Here is the result using Gillmore's method for the count rate. I varied the range of data used to count signal events. As the signal range was varied, the background range changed so that the number of background bins on EACH side of the signal was equal to one-half the number of signal bins.

The red line is the result of the area using the raw fit (above). I believe this is a good value, though I don't yet know its uncertainty. Notice when the signal full range is too small, we aren't capturing points in the data that contain the signal, and hence the Gillmore method gives a count rate that is significantly below the raw fit method. At a full range of 12 keV, we get our best estimate using the Gillmore method.

Making the fit larger than 12 keV is not wise, since adding more data points will increase the uncertainty in the count rate.

Here is the uncertainty in those data.

Notice at a signal range of 12 keV, the uncertainty is $$\pm 3.95$$. This is significantly greater the raw fit method, and that is because we are now accounting for the uncertainty in the background.

## Paul's method

I developed an approach that uses a fit and has only one free parameter that can be tuned to give the lowest uncertainty.

As with the raw fit method, we'll fit the data to the function $$I(E) = A e^{−\frac{1}{2}\frac{(E−E_0)^2}{σ^2}}+m(E−E_0)+b$$

Now recognize that the total number of counts $$N$$ is the integral $$\sum_{i=imin}^{imax}{n_i} = N = \int_{Ei}^{Ef}{\frac{I(E)}{W}dE}$$

The left hand side of that equation is simply the sum of the counts in each bin of data! And the right hand side we can write down with some calculus

$$\sum_{i=imin}^{imax}{n_i} = \frac{A\sigma\sqrt{2\pi}}{W} + \frac{\frac{m}{2}(E_f^2-E_i^2)-(mE_0-b)(E_f-E_i)}{W}$$

That first term on the right side is exactly what we are interested in! Let $$R_P = \frac{A\sigma\sqrt{2\pi}}{W}$$. Then

$$R_P = \sum_{i=imin}^{imax}{n_i} - \frac{\frac{m}{2}(E_f^2-E_i^2)-(mE_0-b)(E_f-E_i)}{W}$$

We can calculate this rate. Its uncertainty is $$\Delta_{RP} = \sqrt{\Delta_{\Sigma ni}^2 + (\frac{\partial R_P}{\partial m}\Delta_m)^2 + (\frac{\partial R_P}{\partial E_0}\Delta_{E0})^2 + (\frac{\partial R_P}{\partial b}\Delta_b)^2 + 2(|(\frac{\partial R_P}{\partial m}\frac{\partial R_P}{\partial E_0}\Delta_{mE0}^2)| + |(\frac{\partial R_P}{\partial m}\frac{\partial R_P}{\partial b}\Delta_{mb}^2)| + |(\frac{\partial R_P}{\partial b}\frac{\partial R_P}{\partial E_0}\Delta_{bE0}^2)|)}$$

The factors of $$\Delta_{ij}^2$$ are the elements of the covariance matrix, which is reported by the fitting algorithm. The term $$\Delta_{\Sigma ni}^2$$, the variance in the total number of counts, IS the total number of counts $$N$$. After evaluating, this can be simplified quite a bit. I can go into detail if anyone is interested. The end result is

The last step is to make the fitting range symmetric with full width $$K$$, so that $$E_i = E_0 - K/2$$ and $$E_f = E_0 + K/2$$. Then

Now we have one free parameter, $$K$$. I fit the data over a range of ranges, and calculate the area and uncertainty for each range. Here is the calculated signal count rate:

The trend in the uncertainty looks like this:

We see that having a full fitting range of 22 keV produces the best fit (this nearly matches the best range of 24 keV using the Gillmore method). The resulting fit estimates the signal count rate as $$R_p = (9.77 \pm 3.99) Hz$$. Notice the uncertainty in the count rate is comparable to the result using the Gillmore method. But the signal count rate is larger than both of the other two methods. I'm still investigating why this might be.

I like my (Paul's) method, because there is I only need to minimize the uncertainty by choosing the optimal fit width. With Gillmore's method, I need to set the fit width based upon an assumed 'correct' signal count rate.

It may be most appropriate to calculate the signal count rate using the raw fit method, and then estimate the uncertainty using Paul's method.

Thanks for reading. I hope you found this insightful.