Source code for gammapy.estimators.energydependentmorphology
# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Implementation of energy-dependent morphology estimator tool."""importnumpyasnpfromgammapy.datasetsimportDatasetsfromgammapy.modelingimportFitfromgammapy.modeling.modelsimportFoVBackgroundModel,Modelsfromgammapy.modeling.selectionimportTestStatisticNestedfromgammapy.stats.utilsimportts_to_sigmafrom.coreimportEstimator__all__=["weighted_chi2_parameter","EnergyDependentMorphologyEstimator"]defweighted_chi2_parameter(results_edep,parameters=["sigma"]):r"""Calculate the weighted chi-squared value for the parameters of interest. The chi-squared parameter is defined as follows: .. math:: \chi^2 = \sum_i \frac{(x_i - \bar{\mu})^2}{\sigma_i ^ 2} where the :math:`x_i` and :math:`\sigma_i` are the value and error of the parameter of interest, and the weighted mean is .. math:: \bar{\mu} = \sum_i \frac{(x_i/ \sigma_i ^ 2)}{(1/\sigma_i ^ 2)} Parameters ---------- result_edep : `dict` Dictionary of results for the energy-dependent estimator. parameters : list of str, optional The model parameters to calculate the chi-squared value for. Default is ["sigma"]. Returns ------- chi2_result : `dict` Dictionary with the chi-squared values for the parameters of interest. Notes ----- This chi-square should be utilised with caution as it does not take into account any correlation between the parameters. To properly utilise the chi-squared parameter one must ensure each of the parameters are independent, which cannot be guaranteed in this use case. """chi2_value=[]df=[]forparameterinparameters:values=results_edep[parameter][1:]errors=results_edep[f"{parameter}_err"][1:]weights=1/errors**2avg=np.average(values,weights=weights)chi2_value+=[np.sum((values-avg)**2/errors**2).to_value()]df+=[len(values)-1]significance=[ts_to_sigma(chi2_value[i],df[i])foriinrange(len(chi2_value))]chi2_result={}chi2_result["parameter"]=parameterschi2_result["chi2"]=chi2_valuechi2_result["df"]=dfchi2_result["significance"]=significancereturnchi2_result
[docs]classEnergyDependentMorphologyEstimator(Estimator):"""Test if there is any energy-dependent morphology in a map dataset for a given set of energy bins. Parameters ---------- energy_edges : `~astropy.units.Quantity` Energy edges for the energy-dependence test. source : str or int For which source in the model to compute the estimator. fit : `~gammapy.modeling.Fit`, optional Fit instance specifying the backend and fit options. If None, the fit backend default is minuit. Default is None. """tag="EnergyDependentMorphologyEstimator"def__init__(self,energy_edges,source=0,fit=None):self.energy_edges=energy_edgesself.source=sourceself.num_energy_bands=len(self.energy_edges)-1iffitisNone:fit=Fit()self.fit=fitdef_slice_datasets(self,datasets):"""Calculate a dataset for each energy slice. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Input datasets to use. Returns ------- slices_src : `~gammapy.datasets.Datasets` Sliced datasets. """model=datasets.models[self.source]filtered_names=[namefornameindatasets.models.namesifname!=self.source]other_models=Models()fornameinfiltered_names:other_models.append(datasets.models[name])slices_src=Datasets()fori,(emin,emax)inenumerate(zip(self.energy_edges[:-1],self.energy_edges[1:])):fordatasetindatasets:sliced_src=dataset.slice_by_energy(emin,emax,name=f"{self.source}_{i}")bkg_sliced_model=FoVBackgroundModel(dataset_name=sliced_src.name)sliced_src.models=[model.copy(name=f"{sliced_src.name}-model"),*other_models,bkg_sliced_model,]slices_src.append(sliced_src)returnslices_srcdef_estimate_source_significance(self,datasets):"""Estimate the significance of the source above the background. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Input datasets to use. Returns ------- result_bkg_src : `dict` Dictionary with the results of the null hypothesis with no source, and alternative hypothesis with the source added in. Entries are: * "Emin" : the minimum energy of the energy band * "Emax" : the maximum energy of the energy band * "delta_ts" : difference in ts * "df" : the degrees of freedom between null and alternative hypothesis * "significance" : significance of the result """slices_src=self._slice_datasets(datasets)# Norm is free and fittest_results=[]forslicedinslices_src:parameters=[paramforparaminsliced.models[f"{sliced.name}-model"].parameters.free_parameters]null_values=[0]+[param.valueforparaminsliced.models[f"{sliced.name}-model"].spatial_model.parameters.free_parameters]test=TestStatisticNested(parameters=parameters,null_values=null_values,n_sigma=-np.inf,fit=self.fit,)test_results.append(test.run(sliced))delta_ts_bkg_src=[_["ts"]for_intest_results]df_src=[len(_["fit_results"].parameters.free_parameters.names)for_intest_results]df_bkg=1df_bkg_src=df_src[0]-df_bkgsigma_ts_bkg_src=ts_to_sigma(delta_ts_bkg_src,df=df_bkg_src)# Prepare results dictionary for signal above backgroundresult_bkg_src={}result_bkg_src["Emin"]=self.energy_edges[:-1]result_bkg_src["Emax"]=self.energy_edges[1:]result_bkg_src["delta_ts"]=delta_ts_bkg_srcresult_bkg_src["df"]=[df_bkg_src]*self.num_energy_bandsresult_bkg_src["significance"]=[elemforeleminsigma_ts_bkg_src]returnresult_bkg_src
[docs]defestimate_energy_dependence(self,datasets):"""Estimate the potential of energy-dependent morphology. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Input datasets to use. Returns ------- results : `dict` Dictionary with results of the energy-dependence test. Entries are: * "delta_ts" : difference in ts between fitting each energy band individually (sliced fit) and the joint fit * "df" : the degrees of freedom between fitting each energy band individually (sliced fit) and the joint fit * "result" : the results for the fitting each energy band individually (sliced fit) and the joint fit """model=datasets.models[self.source]# Calculate the individually sliced componentsslices_src=self._slice_datasets(datasets)results_src=[]forslicedinslices_src:results_src.append(self.fit.run(sliced))results_src_total_stat=[result.total_statforresultinresults_src]free_x,free_y=np.shape([result.parameters.free_parameters.namesforresultinresults_src])df_src=free_x*free_y# Calculate the joint fitparameters=model.spatial_model.parameters.free_parameters.namesslice0=slices_src[0]fori,slice_jinenumerate(slices_src[1:]):forparaminparameters:setattr(slice_j.models[f"{self.source}_{i+1}-model"].spatial_model,param,slice0.models[f"{self.source}_0-model"].spatial_model.parameters[param],)result_joint=self.fit.run(slices_src)# Compare fit of individual energy slices to the results with joint fitdelta_ts_joint=result_joint.total_stat-np.sum(results_src_total_stat)df_joint=len(slices_src.parameters.free_parameters.names)df=df_src-df_joint# Prepare results dictionaryjoint_values=[result_joint.parameters[param].valueforparaminparameters]joint_errors=[result_joint.parameters[param].errorforparaminparameters]parameter_values=np.empty((len(parameters),self.num_energy_bands))parameter_errors=np.empty((len(parameters),self.num_energy_bands))foriinrange(self.num_energy_bands):parameter_values[:,i]=[results_src[i].parameters[param].valueforparaminparameters]parameter_errors[:,i]=[results_src[i].parameters[param].errorforparaminparameters]result={}result["Hypothesis"]=["H0"]+["H1"]*self.num_energy_bandsresult["Emin"]=np.append(self.energy_edges[0],self.energy_edges[:-1])result["Emax"]=np.append(self.energy_edges[-1],self.energy_edges[1:])units=[result_joint.parameters[param].unitforparaminparameters]# Results for H0 in the first row and then H1 -- i.e. individual bands in other rowsforiinrange(len(parameters)):result[f"{parameters[i]}"]=np.append(joint_values[i]*units[i],parameter_values[i]*units[i])result[f"{parameters[i]}_err"]=np.append(joint_errors[i]*units[i],parameter_errors[i]*units[i])returndict(delta_ts=delta_ts_joint,df=df,result=result)
[docs]defrun(self,datasets):"""Run the energy-dependence estimator. Parameters ---------- datasets : `~gammapy.datasets.Datasets` Input datasets to use. Returns ------- results : `dict` Dictionary with the various energy-dependence estimation values. """ifnotisinstance(datasets,Datasets)ordatasets.is_all_same_typeisFalse:raiseValueError("Unsupported datasets type.")results={}results["energy_dependence"]=self.estimate_energy_dependence(datasets)results["src_above_bkg"]=self._estimate_source_significance(datasets)returnresults