Add truncation to cstat as well
@@ -8,8 +8,6 @@ from gammapy.stats.fit_statistics_cython import TRUNCATION_VALUE
__all__ = ["cash", "cstat", "wstat", "get_wstat_mu_bkg", "get_wstat_gof_terms"]
-N_ON_MIN = 1e-25
-
def cash(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
r"""Cash statistic, for Poisson data.
@@ -30,6 +28,7 @@ def cash(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
truncation_value : array_like
Minimum value use for ``mu_on``
``mu_on`` = ``truncation_value`` where ``n_on`` <= ``truncation_value``
+ Default is 1e-25.
Returns
-------
@@ -59,7 +58,7 @@ def cash(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
return stat
-def cstat(n_on, mu_on, n_on_min=N_ON_MIN):
+def cstat(n_on, mu_on, truncation_value=TRUNCATION_VALUE):
r"""C statistic, for Poisson data.
The C statistic is defined as
@@ -70,7 +69,7 @@ def cstat(n_on, mu_on, n_on_min=N_ON_MIN):
and :math:`C = 0` where :math:`\mu_{on} <= 0`.
- ``n_on_min`` handles the case where ``n_on`` is 0 or less and
+ ``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`
@@ -80,8 +79,10 @@ def cstat(n_on, mu_on, n_on_min=N_ON_MIN):
Observed counts
mu_on : array_like
Expected counts
- n_on_min : array_like
- ``n_on`` = ``n_on_min`` where ``n_on`` <= ``n_on_min.``
+ 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
-------
@@ -99,9 +100,10 @@ def cstat(n_on, mu_on, n_on_min=N_ON_MIN):
"""
n_on = np.asanyarray(n_on, dtype=np.float64)
mu_on = np.asanyarray(mu_on, dtype=np.float64)
- n_on_min = np.asanyarray(n_on_min, dtype=np.float64)
+ truncation_value = np.asanyarray(truncation_value, dtype=np.float64)
- n_on = np.where(n_on <= n_on_min, n_on_min, n_on)
+ 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)
