Fit statistics¶
Introduction¶
This page describes common fit statistics used in gamma-ray astronomy. Results were tested against results from the Sherpa and XSpec X-ray analysis packages.
All functions compute per-bin statistics. If you want the summed statistics for
all bins, call sum on the output array yourself. Here’s an example for the
cash statistic:
>>> from gammapy.stats import cash
>>> data = [3, 5, 9]
>>> model = [3.3, 6.8, 9.2]
>>> cash(data, model)
array([ -0.56353481, -5.56922612, -21.54566271])
>>> cash(data, model).sum()
-27.678423645645118
Gaussian data¶
TODO
Poisson data¶
TODO
Poisson data with background measurement¶
If you not only have a measurement of counts \(n_{\mathrm{on}}\) in the signal region, but also a measurement \(n_{\mathrm{off}}\) in a background region you can write down the likelihood formula as
where \(\mu_{\mathrm{sig}}\) is the number of expected counts in the signal regions, and \(\mu_{\mathrm{bkg}}\) is the number of expected counts in the background region, as defined in the Introduction. By taking two time the negative log likelihood and neglecting model independent and thus constant terms, we define the WStat.
In the most general case, where \(\mu_{\mathrm{src}}\) and \(\mu_{\mathrm{bkg}}\) are free the minimum of \(W\) is at
Profile Likelihood¶
Most of the times you probably won’t have a model in order to get \(\mu_{\mathrm{bkg}}\). The strategy in this case is to treat \(\mu_{\mathrm{bkg}}\) as so-called nuisance parameter, i.e. a free parameter that is of no physical interest. Of course you don’t want an additional free parameter for each bin during a fit. Therefore one calculates an estimator for \(\mu_{\mathrm{bkg}}\) by analytically minimizing the likelihood function. This is called ‘profile likelihood’.
This yields a quadratic equation for \(\mu_{\mathrm{bkg}}\)
with the solution
where
Goodness of fit¶
The best-fit value of the WStat as defined now contains no information about the goodness of the fit. We consider the likelihood of the data \(n_{\mathrm{on}}\) and \(n_{\mathrm{off}}\) under the expectation of \(n_{\mathrm{on}}\) and \(n_{\mathrm{off}}\),
and add twice the log likelihood
to WStat. In doing so, we are computing the likelihood ratio:
Intuitively, this log-likelihood ratio should asymptotically behave like a
chi-square with m-n degrees of freedom, where m is the number of
measurements and n the number of model parameters.
Final result¶
Special cases¶
The above formula is undefined if \(n_{\mathrm{on}}\) or \(n_{\mathrm{off}}\) are equal to zero, because of the \(n\log{{n}}\) terms, that were introduced by adding the goodness of fit terms. These cases are treated as follows.
If \(n_{\mathrm{on}} = 0\) the likelihood formulae read
and
WStat is derived by taking 2 times the negative log likelihood and adding the goodness of fit term as ever
Note that this is the limit of the original Wstat formula for \(n_{\mathrm{on}} \rightarrow 0\).
The analytical result for \(\mu_{\mathrm{bkg}}\) in this case reads:
When inserting this into the WStat we find the simplified expression.
If \(n_{\mathrm{off}} = 0\) Wstat becomes
and
For \(\mu_{\mathrm{sig}} > n_{\mathrm{on}} (\frac{\alpha}{1 + \alpha})\), \(\mu_{\mathrm{bkg}}\) becomes negative which is unphysical.
Therefore we distinct two cases. The physical one where
\(\mu_{\mathrm{sig}} < n_{\mathrm{on}} (\frac{\alpha}{1 + \alpha})\).
is straightforward and gives
For the unphysical case, we set \(\mu_{\mathrm{bkg}}=0\) and arrive at
Example¶
The following table gives an overview over values that WStat takes in different scenarios
>>> from gammapy.stats import wstat
>>> from astropy.table import Table
>>> table = Table()
>>> table['mu_sig'] = [0.1, 0.1, 1.4, 0.2, 0.1, 5.2, 6.2, 4.1, 6.4, 4.9, 10.2,
... 16.9, 102.5]
>>> table['n_on'] = [0, 0, 0, 0, 0, 5, 5, 5, 5, 5, 10, 20, 100]
>>> table['n_off'] = [0, 1, 1, 10 , 10, 0, 5, 5, 20, 40, 2, 70, 10]
>>> table['alpha'] = [0.01, 0.01, 0.5, 0.1 , 0.2, 0.2, 0.2, 0.01, 0.4, 0.4,
... 0.2, 0.1, 0.6]
>>> table['wstat'] = wstat(n_on=table['n_on'],
... n_off=table['n_off'],
... alpha=table['alpha'],
... mu_sig=table['mu_sig'])
>>> table['wstat'].format = '.3f'
>>> table.pprint()
mu_sig n_on n_off alpha wstat
------ ---- ----- ----- ------
0.1 0 0 0.01 0.200
0.1 0 1 0.01 0.220
1.4 0 1 0.5 3.611
0.2 0 10 0.1 2.306
0.1 0 10 0.2 3.846
5.2 5 0 0.2 0.008
6.2 5 5 0.2 0.736
4.1 5 5 0.01 0.163
6.4 5 20 0.4 7.125
4.9 5 40 0.4 14.578
10.2 10 2 0.2 0.034
16.9 20 70 0.1 0.656
102.5 100 10 0.6 0.663
Notes¶
All above formulae are equivalent to what is given on the XSpec manual statistics page with the substitutions: