# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Common fit statistics used in gamma-ray astronomy.see :ref:`fit-statistics`"""importnumpyasnpfromgammapy.stats.fit_statistics_cythonimportTRUNCATION_VALUE__all__=["cash","cstat","wstat","get_wstat_mu_bkg","get_wstat_gof_terms"]
[docs]defcash(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 ``mu_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)ifnp.any(truncation_value)<=0:raiseValueError("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 belowwithnp.errstate(divide="ignore",invalid="ignore"):stat=2*(mu_on-n_on*np.log(mu_on))returnstat
[docs]defcstat(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)ifnp.any(truncation_value)<=0:raiseValueError("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)returnstat
[docs]defwstat(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 is 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.gsfc.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 = alphan_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)ifmu_bkgisNone: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 belowwithnp.errstate(divide="ignore",invalid="ignore"):# This is a false positive error from pylint# See https://github.com/PyCQA/pylint/issues/2436term2_=-n_on*np.log(mu_sig+alpha*mu_bkg)# pylint:disable=invalid-unary-operand-type# Handle n_on == 0condition=n_on==0term2=np.where(condition,0,term2_)# suppress zero division warnings, they are corrected belowwithnp.errstate(divide="ignore",invalid="ignore"):# This is a false positive error from pylint# See https://github.com/PyCQA/pylint/issues/2436term3_=-n_off*np.log(mu_bkg)# pylint:disable=invalid-unary-operand-type# Handle n_off == 0condition=n_off==0term3=np.where(condition,0,term3_)stat=2*(term1+term2+term3)ifextra_terms:stat+=get_wstat_gof_terms(n_on,n_off)returnstat
[docs]defget_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 formulaC=alpha*(n_on+n_off)-(1+alpha)*mu_sigD=np.sqrt(C**2+4*alpha*(alpha+1)*n_off*mu_sig)withnp.errstate(invalid="ignore",divide="ignore"):mu_bkg=(C+D)/(2*alpha*(alpha+1))returnmu_bkg
[docs]defget_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 belowwithnp.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)return2*term