# Licensed under a 3-clause BSD style license - see LICENSE.rstimportabcimporthtmlimportnumpyasnpfromscipy.specialimportlambertwfromscipy.statsimportchi2fromgammapy.utils.rootsimportfind_rootsfrom.fit_statisticsimportcash,wstat__all__=["WStatCountsStatistic","CashCountsStatistic"]classCountsStatistic(abc.ABC):"""Counts statistics base class."""@property@abc.abstractmethoddefstat_null(self):pass@property@abc.abstractmethoddefstat_max(self):pass@property@abc.abstractmethoddefn_sig(self):pass@property@abc.abstractmethoddefn_bkg(self):pass@property@abc.abstractmethoddeferror(self):pass@abc.abstractmethoddef_stat_fcn(self):pass@propertydefts(self):"""Return stat difference (TS) of measured excess versus no excess."""# Remove (small) negative TS due to error in root findingts=np.clip(self.stat_null-self.stat_max,0,None)returnts@propertydefsqrt_ts(self):"""Return statistical significance of measured excess. The sign of the excess is applied to distinguish positive and negative fluctuations. """returnnp.sign(self.n_sig)*np.sqrt(self.ts)@propertydefp_value(self):"""Return p_value of measured excess. Here the value accounts only for the positive excess significance (i.e. one-sided). """return0.5*chi2.sf(self.ts,1)def__str__(self):str_="\t{:32}: {{n_on:.2f}} \n".format("Total counts")str_+="\t{:32}: {{background:.2f}}\n".format("Total background counts")str_+="\t{:32}: {{excess:.2f}}\n".format("Total excess counts")str_+="\t{:32}: {{significance:.2f}}\n".format("Total significance")str_+="\t{:32}: {{p_value:.3f}}\n".format("p - value")str_+="\t{:32}: {{n_bins:.0f}}\n".format("Total number of bins")info=self.info_dict()info["n_bins"]=np.array(self.n_on).sizestr_=str_.format(**info)returnstr_.expandtabs(tabsize=2)def_repr_html_(self):try:returnself.to_html()exceptAttributeError:returnf"<pre>{html.escape(str(self))}</pre>"definfo_dict(self):"""A dictionary of the relevant quantities. Returns ------- info_dict : dict Dictionary with summary information. """info_dict={}info_dict["n_on"]=self.n_oninfo_dict["background"]=self.n_bkginfo_dict["excess"]=self.n_siginfo_dict["significance"]=self.sqrt_tsinfo_dict["p_value"]=self.p_valuereturninfo_dictdefcompute_errn(self,n_sigma=1.0):"""Compute downward excess uncertainties. Searches the signal value for which the test statistics is n_sigma**2 away from the maximum. Parameters ---------- n_sigma : float Confidence level of the uncertainty expressed in number of sigma. Default is 1. """errn=np.zeros_like(self.n_sig,dtype="float")min_range=self.n_sig-2*n_sigma*(self.error+1)it=np.nditer(errn,flags=["multi_index"])whilenotit.finished:roots,res=find_roots(self._stat_fcn,min_range[it.multi_index],self.n_sig[it.multi_index],nbin=1,args=(self.stat_max[it.multi_index]+n_sigma**2,it.multi_index),)ifnp.isnan(roots[0]):errn[it.multi_index]=self.n_on[it.multi_index]else:errn[it.multi_index]=self.n_sig[it.multi_index]-roots[0]it.iternext()returnerrndefcompute_errp(self,n_sigma=1):"""Compute upward excess uncertainties. Searches the signal value for which the test statistics is n_sigma**2 away from the maximum. Parameters ---------- n_sigma : float Confidence level of the uncertainty expressed in number of sigma. Default is 1. """errp=np.zeros_like(self.n_on,dtype="float")max_range=self.n_sig+2*n_sigma*(self.error+1)it=np.nditer(errp,flags=["multi_index"])whilenotit.finished:roots,res=find_roots(self._stat_fcn,self.n_sig[it.multi_index],max_range[it.multi_index],nbin=1,args=(self.stat_max[it.multi_index]+n_sigma**2,it.multi_index),)errp[it.multi_index]=roots[0]it.iternext()returnerrp-self.n_sigdefcompute_upper_limit(self,n_sigma=3):"""Compute upper limit on the signal. Searches the signal value for which the test statistics is n_sigma**2 away from the maximum or from 0 if the measured excess is negative. Parameters ---------- n_sigma : float Confidence level of the upper limit expressed in number of sigma. Default is 3. """ul=np.zeros_like(self.n_on,dtype="float")min_range=self.n_sigmax_range=min_range+2*n_sigma*(self.error+1)it=np.nditer(ul,flags=["multi_index"])whilenotit.finished:ts_ref=self._stat_fcn(min_range[it.multi_index],0.0,it.multi_index)roots,res=find_roots(self._stat_fcn,min_range[it.multi_index],max_range[it.multi_index],nbin=1,args=(ts_ref+n_sigma**2,it.multi_index),)ul[it.multi_index]=roots[0]it.iternext()returnul@abc.abstractmethoddef_n_sig_matching_significance_fcn(self):passdefn_sig_matching_significance(self,significance):"""Compute excess matching a given significance. This function is the inverse of `significance`. Parameters ---------- significance : float Significance. Returns ------- n_sig : `numpy.ndarray` Excess. """n_sig=np.zeros_like(self.n_bkg,dtype="float")it=np.nditer(n_sig,flags=["multi_index"])whilenotit.finished:lower_bound=np.sqrt(self.n_bkg[it.multi_index])*significance# find upper bounds for secant method as in scipyeps=1e-4upper_bound=lower_bound*(1+eps)upper_bound+=epsifupper_bound>=0else-epsroots,res=find_roots(self._n_sig_matching_significance_fcn,lower_bound=lower_bound,upper_bound=upper_bound,args=(significance,it.multi_index),nbin=1,method="secant",)n_sig[it.multi_index]=roots[0]# return NaN if failit.iternext()returnn_sig@abc.abstractmethoddefsum(self,axis=None):"""Return summed CountsStatistics. Parameters ---------- axis : None or int or tuple of ints, optional Axis or axes on which to perform the summation. Default, axis=None, will perform the sum over the whole array. Returns ------- stat : `~gammapy.stats.CountsStatistics` The summed stat object. """pass
[docs]classCashCountsStatistic(CountsStatistic):"""Class to compute statistics for Poisson distributed variable with known background. Parameters ---------- n_on : int Measured counts. mu_bkg : float Known level of background. """def__init__(self,n_on,mu_bkg):self.n_on=np.asanyarray(n_on)self.mu_bkg=np.asanyarray(mu_bkg)@propertydefn_bkg(self):"""Expected background counts."""returnself.mu_bkg@propertydefn_sig(self):"""Excess."""returnself.n_on-self.n_bkg@propertydeferror(self):"""Approximate error from the covariance matrix."""returnnp.sqrt(self.n_on)@propertydefstat_null(self):"""Stat value for null hypothesis, i.e. 0 expected signal counts."""returncash(self.n_on,self.mu_bkg+0)@propertydefstat_max(self):"""Stat value for best fit hypothesis, i.e. expected signal mu = n_on - mu_bkg."""returncash(self.n_on,self.n_on)
[docs]definfo_dict(self):"""A dictionary of the relevant quantities. Returns ------- info_dict : dict Dictionary with summary info. """info_dict=super().info_dict()info_dict["mu_bkg"]=self.mu_bkgreturninfo_dict
[docs]classWStatCountsStatistic(CountsStatistic):"""Class to compute statistics for Poisson distributed variable with unknown background. Parameters ---------- n_on : int Measured counts in on region. n_off : int Measured counts in off region. alpha : float Acceptance ratio of on and off measurements. mu_sig : float Expected signal counts in on region. """def__init__(self,n_on,n_off,alpha,mu_sig=None):self.n_on=np.asanyarray(n_on)self.n_off=np.asanyarray(n_off)self.alpha=np.asanyarray(alpha)ifmu_sigisNone:self.mu_sig=np.zeros_like(self.n_on)else:self.mu_sig=np.asanyarray(mu_sig)@propertydefn_bkg(self):"""Known background computed alpha * n_off."""returnself.alpha*self.n_off@propertydefn_sig(self):"""Excess"""returnself.n_on-self.n_bkg-self.mu_sig@propertydeferror(self):"""Approximate error from the covariance matrix."""returnnp.sqrt(self.n_on+self.alpha**2*self.n_off)@propertydefstat_null(self):"""Stat value for null hypothesis, i.e. mu_sig expected signal counts."""returnwstat(self.n_on,self.n_off,self.alpha,self.mu_sig)@propertydefstat_max(self):"""Stat value for best fit hypothesis. i.e. expected signal mu = n_on - alpha * n_off - mu_sig. """returnwstat(self.n_on,self.n_off,self.alpha,self.n_sig+self.mu_sig)
[docs]definfo_dict(self):"""A dictionary of the relevant quantities. Returns ------- info_dict : dict Dictionary with summary info. """info_dict=super().info_dict()info_dict["n_off"]=self.n_offinfo_dict["alpha"]=self.alphainfo_dict["mu_sig"]=self.mu_sigreturninfo_dict
def__str__(self):str_=f"{self.__class__.__name__}\n"str_+=super().__str__()info_dict=self.info_dict()str_+="\t{:32}: {:.2f}\n".format("Off counts",info_dict["n_off"])str_+="\t{:32}: {:.2f}\n".format("alpha ",info_dict["alpha"])str_+="\t{:32}: {:.2f}\n".format("Predicted signal counts",info_dict["mu_sig"])returnstr_.expandtabs(tabsize=2)def_stat_fcn(self,mu,delta=0,index=None):return(wstat(self.n_on[index],self.n_off[index],self.alpha[index],(mu+self.mu_sig[index]),)-delta)def_n_sig_matching_significance_fcn(self,n_sig,significance,index):stat0=wstat(n_sig+self.n_bkg[index],self.n_off[index],self.alpha[index],0)stat1=wstat(n_sig+self.n_bkg[index],self.n_off[index],self.alpha[index],n_sig,)returnnp.sign(n_sig)*np.sqrt(np.clip(stat0-stat1,0,None))-significance