Statistics error limit on a fitting We are trying to analyze gamma-ray spectra that come from a neutron scattering in the soil sample. Our hope is to find the final error in the fitting parameters of the peak, which should look like a gaussian shape. We expect to have a low count (s) of the actual signal on top of a background (b).
From a simple calculation, error on the signal will be $\sqrt(s + 2b)$ (as in the second half of this post). I am wondering if this is the lowest limit of the error estimation for signal s. If so, when s is small, our error can be larger than the value s. If anyone has a thought on this, whether a confirmation that this is a correct statement and/or a way to calculate errors without being affected by the background value, I would love to hear about it. Thanks!
The calculation and MWE for the signal s that lead to this question is as follows:
Adding two Poissonian results in a new Poissonian
Given $p(i; \mu_1)$ and $p(j; \mu_2)$ we can calculate a new distribution
$$ p(n; \lambda) = \sum^n_{l=0} p(l; \mu_1) p(n-l; \mu_2) $$
which one can simplify to
$$ p(n; \lambda) = \frac{(\mu_1+\mu_2)^k}{k!}\, e^{-(\mu_1+\mu_2)}  $$
So $\lambda=\mu_1+\mu_2$ and the sum is still Poissonian distributed and therefore the error is $\sqrt{\mu_1+\mu_2}$
Difference of two Poissonian
The difference of two Poissonian has a special distribution, the Skellam distribution.
The mean of the Skellam distribution $\mu_1-\mu_2$ and the variance is $\mu_1+\mu_2$.
Application
This has direct application for gamma spectra. If we look at a energy histogram in a gamma measurement, we often have a background signal (b) and a signal (s). Both will be Poissonian distributed. Therefore the sum of both will also be Poissonian with a mean and variance of $s+b$.
If we measure this number and can estimate a background (say from neighboring channels), we can subtract the background out, but this will result in a Skellam distribution and therefore we get:
mean:  $(s+b)-b = s$
variance:  $(s+b)+b = s+2b$  therefore the error is $\sqrt{s+2b}$.
N = 10_000

mu_s = 100
mu_b = 100

measurement = np.random.poisson(mu_s+mu_b, size=N)
background = np.random.poisson(mu_b, size=N)

signal = measurement-background
# need to make sure we get one bin for each integer to avoid binning artifects
plt.hist(signal, bins=120, range=(40, 160))
mu = signal.mean()
sigma = signal.std()
plt.axvline(mu, color="red")
plt.axvline(mu+sigma, color="green")
plt.axvline(mu-sigma, color="green")
plt.xlabel('measured-background')

print(f"measured: {mu} +- {sigma}")
print(f"   sqrt(mu_s+mu_b) = {np.sqrt(mu_s+mu_b)}  <-- understimates")
print(f"   sqrt(mu_s+2*mu_b) = {np.sqrt(mu_s+2*mu_b)}")

Output:
measured: 99.9346 +- 17.360861811557626
   sqrt(mu_s+mu_b) = 14.142135623730951  <-- understimates
   sqrt(mu_s+2*mu_b) = 17.320508075688775


 A: I'm not a particle physisics, however, I believe you are using the wrong tools to subtract the background. What you are doing is the following:

*

*You measure the background distribution and calculate $\lambda_{background}$,

*you use this calculated $\lambda_{background}$ and simulate random numbers from the Poisson distribution. This yields a density distribution, which possesses uncertainty.

*you subtract the simulated density distribution from the measured density distribution. As you have "two random number" in every bin, the uncertainty increases.

I'm pretty sure that step 2 is incorrect. You should not draw random numbers from the Poisson distribution. Instead, you should work with the probability density distribution. This distribution is a well-defined mathematical function. By using this function, we do not obtain an uncertainty. Therefore, subtracting the background from the measured signal does not increase the uncertainty. Does this make sense?
Here my simulated data, to explain how I would do it:

The graph displays the background data
$X_1 \sim Pois(\lambda=70)$, the (in the experiment unknown) signal
$X_2 \sim Pois(\lambda=100)$, and the combination of the background data and the signal,
$X_{total} = X_1 + X_2$. First, I fit the background data and obtain
Fitting of the distribution ' pois ' by maximum likelihood 
Parameters : 
       estimate Std. Error
lambda  70.2215  0.1873785

Next, I fit the total data and obtain
Fitting of the distribution ' pois ' by maximum likelihood 
Parameters : 
       estimate Std. Error
lambda   170.41  0.2918989

Therefore, the point estimate of the signal data is
$E[X_2] = E[X_{total}] - E[X_1] = 170.41 - 70.22 = 100.19$, and its uncertainty is given by
$\sqrt{Var[X_{total}] + Var[X_1]} = \ldots = 0.34$.
A: You have not said how you are estimating $b$. The implication from your $\sqrt{s + 2b}$ result is that you measure source plus background in one channel and background in one other channel. In which case I think your result is correct (at least in a frequentist interpretation, with no prior information about $s$ and $b$).
I would have thought it more likely that your information about $b$ comes from measurements in several ($N$) channels, in which case, the variance in this would not be $b$, it would be $b/N$.
A better (Bayesian) method would be to forward-model the observed source plus background and background channels, calculating the (log) likelihood for each combination of $s$ and $b$ producing the observed data. To this you would add the log of the prior probabilities for $s$ and $b$ -maybe flat, but not allowing either $s$ or $b$ to be less than zero (your current method and error bar does allow $s<0$). The anti-log of the resultants would give the posterior probability distributions for $s$ and $b$.
A: This problem is widely-known and called the 'on-off problem' and occasionally within astrophysics as the Li & Ma problem owing to the paper:
T. P. Li and Y. Q. Ma. Analysis methods for results in gamma-ray astronomy.
Astrophys. J., 272:317–324, 1983
A good starting point for reading more on this might be
A comparison of limit setting methods for the on–off problem
Wolfgang A. Rolke
Nucl.Instrum.Meth.A 806 (2016) 318-324
1505.07027 [physics.data-an]
That paper reviews limit setting in a frequentist context. There are other papers with Bayesian treatments and that focus on discoveries. You can browse this search to see some of the other physics literature on the topic (NB sometimes it is written on-off and sometimes on/off):
https://inspirehep.net/literature?sort=mostrecent&size=25&page=1&q=ft%20%22on-off%20problem%22%20or%20ft%20%22on%2Foff%20problem%22
