Source code for gammapy.stats.fit_statistics

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Common fit statistics used in gamma-ray astronomy.

see :ref:fit-statistics
"""
import numpy as np
from gammapy.stats.fit_statistics_cython import TRUNCATION_VALUE

__all__ = ["cash", "cstat", "wstat", "get_wstat_mu_bkg", "get_wstat_gof_terms"]

[docs]def cash(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
r"""Cash statistic, for Poisson data.

The Cash statistic is defined as:

.. math::
C = 2 \left( \mu_{on} - n_{on} \log \mu_{on} \right)

and :math:C = 0 where :math:\mu <= 0.
For more information see :ref:fit-statistics

Parameters
----------
n_on : array_like
Observed counts
mu_on : array_like
Expected counts
truncation_value : array_like
Minimum value use for mu_on
mu_on = truncation_value where n_on <= truncation_value
Default is 1e-25.

Returns
-------
stat : ndarray
Statistic per bin

References
----------
* Sherpa statistics page section on the Cash statistic
<http://cxc.cfa.harvard.edu/sherpa/statistics/#cash>_
* Sherpa help page on the Cash statistic
<http://cxc.harvard.edu/sherpa/ahelp/cash.html>_
* Cash 1979, ApJ 228, 939
<https://ui.adsabs.harvard.edu/abs/1979ApJ...228..939C>_
"""
n_on = np.asanyarray(n_on)
mu_on = np.asanyarray(mu_on)
truncation_value = np.asanyarray(truncation_value)
if np.any(truncation_value) <= 0:
raise ValueError("Cash statistic truncation value must be positive.")

mu_on = np.where(mu_on <= truncation_value, truncation_value, mu_on)

# suppress zero division warnings, they are corrected below
with np.errstate(divide="ignore", invalid="ignore"):
stat = 2 * (mu_on - n_on * np.log(mu_on))
return stat

[docs]def cstat(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
r"""C statistic, for Poisson data.

The C statistic is defined as

.. math::
C = 2 \left[ \mu_{on} - n_{on} + n_{on}
(\log(n_{on}) - log(\mu_{on}) \right]

and :math:C = 0 where :math:\mu_{on} <= 0.

truncation_value handles the case where n_on or mu_on is 0 or less and
the log cannot be taken.
For more information see :ref:fit-statistics

Parameters
----------
n_on : array_like
Observed counts
mu_on : array_like
Expected counts
truncation_value : array_like
n_on = truncation_value where n_on <= truncation_value.
mu_on = truncation_value where n_on <= truncation_value
Default is 1e-25.

Returns
-------
stat : ndarray
Statistic per bin

References
----------
* Sherpa stats page section on the C statistic
<http://cxc.cfa.harvard.edu/sherpa/statistics/#cstat>_
* Sherpa help page on the C statistic
<http://cxc.harvard.edu/sherpa/ahelp/cash.html>_
* Cash 1979, ApJ 228, 939
<https://ui.adsabs.harvard.edu/abs/1979ApJ...228..939C>_
"""
n_on = np.asanyarray(n_on, dtype=np.float64)
mu_on = np.asanyarray(mu_on, dtype=np.float64)
truncation_value = np.asanyarray(truncation_value, dtype=np.float64)

if np.any(truncation_value) <= 0:
raise ValueError("Cstat statistic truncation value must be positive.")

n_on = np.where(n_on <= truncation_value, truncation_value, n_on)
mu_on = np.where(mu_on <= truncation_value, truncation_value, mu_on)

term1 = np.log(n_on) - np.log(mu_on)
stat = 2 * (mu_on - n_on + n_on * term1)
stat = np.where(mu_on > 0, stat, 0)

return stat

[docs]def wstat(n_on, n_off, alpha, mu_sig, mu_bkg=None, extra_terms=True):
r"""W statistic, for Poisson data with Poisson background.

For a definition of WStat see :ref:wstat. If mu_bkg is not provided
it will be calculated according to the profile likelihood formula.

Parameters
----------
n_on : array_like
Total observed counts
n_off : array_like
Total observed background counts
alpha : array_like
Exposure ratio between on and off region
mu_sig : array_like
Signal expected counts
mu_bkg : array_like, optional
Background expected counts
extra_terms : bool, optional
Add model independent terms to convert stat into goodness-of-fit
parameter, default: True

Returns
-------
stat : ndarray
Statistic per bin

References
----------
* Habilitation M. de Naurois, p. 141
<http://inspirehep.net/record/1122589/files/these_short.pdf>_
* XSPEC page on Poisson data with Poisson background
<https://heasarc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html>_
"""
# Note: This is equivalent to what's defined on the XSPEC page under the
# following assumptions
# t_s * m_i = mu_sig
# t_b * m_b = mu_bkg
# t_s / t_b = alpha

n_on = np.asanyarray(n_on, dtype=np.float64)
n_off = np.asanyarray(n_off, dtype=np.float64)
alpha = np.asanyarray(alpha, dtype=np.float64)
mu_sig = np.asanyarray(mu_sig, dtype=np.float64)

if mu_bkg is None:
mu_bkg = get_wstat_mu_bkg(n_on, n_off, alpha, mu_sig)

term1 = mu_sig + (1 + alpha) * mu_bkg

# suppress zero division warnings, they are corrected below
with np.errstate(divide="ignore", invalid="ignore"):
# This is a false positive error from pylint
# See https://github.com/PyCQA/pylint/issues/2436
term2_ = -n_on * np.log(
mu_sig + alpha * mu_bkg
)  # pylint:disable=invalid-unary-operand-type
# Handle n_on == 0
condition = n_on == 0
term2 = np.where(condition, 0, term2_)

# suppress zero division warnings, they are corrected below
with np.errstate(divide="ignore", invalid="ignore"):
# This is a false positive error from pylint
# See https://github.com/PyCQA/pylint/issues/2436
term3_ = -n_off * np.log(mu_bkg)  # pylint:disable=invalid-unary-operand-type
# Handle n_off == 0
condition = n_off == 0
term3 = np.where(condition, 0, term3_)

stat = 2 * (term1 + term2 + term3)

if extra_terms:
stat += get_wstat_gof_terms(n_on, n_off)

return stat

[docs]def get_wstat_mu_bkg(n_on, n_off, alpha, mu_sig):
"""Background estimate mu_bkg for WSTAT.

See :ref:wstat.
"""
n_on = np.asanyarray(n_on, dtype=np.float64)
n_off = np.asanyarray(n_off, dtype=np.float64)
alpha = np.asanyarray(alpha, dtype=np.float64)
mu_sig = np.asanyarray(mu_sig, dtype=np.float64)

# NOTE: Corner cases in the docs are all handled correctly by this formula
C = alpha * (n_on + n_off) - (1 + alpha) * mu_sig
D = np.sqrt(C ** 2 + 4 * alpha * (alpha + 1) * n_off * mu_sig)
with np.errstate(invalid="ignore", divide="ignore"):
mu_bkg = (C + D) / (2 * alpha * (alpha + 1))

return mu_bkg

[docs]def get_wstat_gof_terms(n_on, n_off):
"""Goodness of fit terms for WSTAT.

See :ref:wstat.
"""
term = np.zeros(n_on.shape)

# suppress zero division warnings, they are corrected below
with np.errstate(divide="ignore", invalid="ignore"):
term1 = -n_on * (1 - np.log(n_on))
term2 = -n_off * (1 - np.log(n_off))

term += np.where(n_on == 0, 0, term1)
term += np.where(n_off == 0, 0, term2)

return 2 * term