# Licensed under a 3-clause BSD style license - see LICENSE.rstimportloggingimportnumpyasnpfromgammapy.datasetsimportDatasetsfromgammapy.datasets.actorsimportDatasetsActorfromgammapy.modelingimportFitfrom.coreimportEstimatorlog=logging.getLogger(__name__)
[docs]classParameterEstimator(Estimator):"""Model parameter estimator. Estimates a model parameter for a group of datasets. Compute best fit value, symmetric and delta(TS) for a given null value. Additionally asymmetric errors as well as parameter upper limit and fit statistic profile can be estimated. Parameters ---------- n_sigma : int Sigma to use for asymmetric error computation. Default is 1. n_sigma_ul : int Sigma to use for upper limit computation. Default is 2. null_value : float Which null value to use for the parameter. selection_optional : list of str, optional Which additional quantities to estimate. Available options are: * "all": all the optional steps are executed. * "errn-errp": estimate asymmetric errors on parameter best fit value. * "ul": estimate upper limits. * "scan": estimate fit statistic profiles. Default is None so the optional steps are not executed. fit : `Fit` Fit instance specifying the backend and fit options. reoptimize : bool Re-optimize other free model parameters. Default is True. """tag="ParameterEstimator"_available_selection_optional=["errn-errp","ul","scan"]def__init__(self,n_sigma=1,n_sigma_ul=2,null_value=1e-150,selection_optional=None,fit=None,reoptimize=True,):self.n_sigma=n_sigmaself.n_sigma_ul=n_sigma_ulself.null_value=null_valueself.selection_optional=selection_optionaliffitisNone:fit=Fit()self.fit=fitself.reoptimize=reoptimize
[docs]defestimate_best_fit(self,datasets,parameter):"""Estimate parameter asymmetric errors. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Datasets. parameter : `Parameter` For which parameter to get the value. Returns ------- result : dict Dictionary with the various parameter estimation values. Entries are: * parameter.name: best fit parameter value. * "stat": best fit total stat. * "success": boolean flag for fit success. * parameter.name_err: covariance-based error estimate on parameter value. """value,total_stat,success,error=np.nan,0.0,False,np.nanifnp.any(datasets.contributes_to_stat):result=self.fit.run(datasets=datasets)value,error=parameter.value,parameter.errortotal_stat=result.optimize_result.total_statsuccess=result.successreturn{f"{parameter.name}":value,"stat":total_stat,"success":success,f"{parameter.name}_err":error*self.n_sigma,}
[docs]defestimate_ts(self,datasets,parameter):"""Estimate parameter ts. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Datasets. parameter : `Parameter` For which parameter to get the value. Returns ------- result : dict Dictionary with the test statistic of the best fit value compared to the null hypothesis. Entries are: * "ts" : fit statistic difference with null hypothesis. * "npred" : predicted number of counts per dataset. * "stat_null" : total stat corresponding to the null hypothesis """npred=self.estimate_npred(datasets=datasets)ifnotnp.any(datasets.contributes_to_stat):stat=np.nannpred["npred"][...]=np.nanelse:stat=datasets.stat_sum()withdatasets.parameters.restore_status():# compute ts valueparameter.value=self.null_valueifself.reoptimize:parameter.frozen=True_=self.fit.optimize(datasets=datasets)ts=datasets.stat_sum()-statstat_null=datasets.stat_sum()return{"ts":ts,"npred":npred["npred"],"stat_null":stat_null}
[docs]defestimate_errn_errp(self,datasets,parameter):"""Estimate parameter asymmetric errors. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Datasets. parameter : `Parameter` For which parameter to get the value. Returns ------- result : dict Dictionary with the parameter asymmetric errors. Entries are: * {parameter.name}_errp : positive error on parameter value. * {parameter.name}_errn : negative error on parameter value. """ifnotnp.any(datasets.contributes_to_stat):return{f"{parameter.name}_errp":np.nan,f"{parameter.name}_errn":np.nan,}self.fit.optimize(datasets=datasets)res=self.fit.confidence(datasets=datasets,parameter=parameter,sigma=self.n_sigma,reoptimize=self.reoptimize,)return{f"{parameter.name}_errp":res["errp"],f"{parameter.name}_errn":res["errn"],}
[docs]defestimate_scan(self,datasets,parameter):"""Estimate parameter statistic scan. Parameters ---------- datasets : `~gammapy.datasets.Datasets` The datasets used to estimate the model parameter. parameter : `~gammapy.modeling.Parameter` For which parameter to get the value. Returns ------- result : dict Dictionary with the parameter fit scan values. Entries are: * parameter.name_scan : parameter values scan. * "stat_scan" : fit statistic values scan. """scan_values=parameter.scan_valuesifnotnp.any(datasets.contributes_to_stat):return{f"{parameter.name}_scan":scan_values,"stat_scan":scan_values*np.nan,}self.fit.optimize(datasets=datasets)profile=self.fit.stat_profile(datasets=datasets,parameter=parameter,reoptimize=self.reoptimize)return{f"{parameter.name}_scan":scan_values,"stat_scan":profile["stat_scan"],}
[docs]defestimate_ul(self,datasets,parameter):"""Estimate parameter ul. Parameters ---------- datasets : `~gammapy.datasets.Datasets` The datasets used to estimate the model parameter. parameter : `~gammapy.modeling.Parameter` For which parameter to get the value. Returns ------- result : dict Dictionary with the parameter upper limits. Entries are: * parameter.name_ul : upper limit on parameter value. """ifnotnp.any(datasets.contributes_to_stat):return{f"{parameter.name}_ul":np.nan}self.fit.optimize(datasets=datasets)res=self.fit.confidence(datasets=datasets,parameter=parameter,sigma=self.n_sigma_ul,reoptimize=self.reoptimize,)return{f"{parameter.name}_ul":res["errp"]+parameter.value}
[docs]@staticmethoddefestimate_counts(datasets):"""Estimate counts for the flux point. Parameters ---------- datasets : Datasets Datasets. Returns ------- result : dict Dictionary with an array with one entry per dataset with the sum of the masked counts. """counts=[]fordatasetindatasets:mask=dataset.maskcounts.append(dataset.counts.data[mask].sum())return{"counts":np.array(counts,dtype=int),"datasets":datasets.names}
[docs]@staticmethoddefestimate_npred(datasets):"""Estimate npred for the flux point. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Datasets. Returns ------- result : dict Dictionary with an array with one entry per dataset with the sum of the masked npred. """npred=[]fordatasetindatasets:mask=dataset.masknpred.append(dataset.npred().data[mask].sum())return{"npred":np.array(npred),"datasets":datasets.names}
[docs]defrun(self,datasets,parameter):"""Run the parameter estimator. Parameters ---------- datasets : `~gammapy.datasets.Datasets` The datasets used to estimate the model parameter. parameter : `str` or `~gammapy.modeling.Parameter` For which parameter to run the estimator. Returns ------- result : dict Dictionary with the various parameter estimation values. """ifnotisinstance(datasets,DatasetsActor):datasets=Datasets(datasets)parameter=datasets.parameters[parameter]withdatasets.parameters.restore_status():ifnotself.reoptimize:datasets.parameters.freeze_all()parameter.frozen=Falseresult=self.estimate_best_fit(datasets,parameter)result.update(self.estimate_ts(datasets,parameter))if"errn-errp"inself.selection_optional:result.update(self.estimate_errn_errp(datasets,parameter))if"ul"inself.selection_optional:result.update(self.estimate_ul(datasets,parameter))if"scan"inself.selection_optional:result.update(self.estimate_scan(datasets,parameter))result.update(self.estimate_counts(datasets))returnresult