select_nested_models#

gammapy.modeling.select_nested_models(datasets, parameters, null_values, n_sigma=2, n_free_parameters=None, fit=None)[source]#

Compute the test statistic (TS) between two nested hypothesis.

The null hypothesis is the minimal one, for which a set of parameters are frozen to given values. The model is updated to the alternative hypothesis if there is a significant improvement (larger than the given threshold).

Parameters:
datasetsDatasets

Datasets.

parametersParameters or list of Parameter

List of parameters frozen for the null hypothesis but free for the test hypothesis.

null_valueslist of float or Parameters

Values of the parameters frozen for the null hypothesis. If a Parameters object or a list of Parameters is given the null hypothesis follows the values of these parameters, so this tests linked parameters versus unliked.

n_sigmafloat, optional

Threshold in number of sigma to switch from the null hypothesis to the alternative one. Default is 2. The TS is converted to sigma assuming that the Wilk’s theorem is verified.

n_free_parametersint, optional

Number of free parameters to consider between the two hypothesis in order to estimate the ts_threshold from the n_sigma threshold. Default is len(parameters).

fitFit, optional

Fit instance specifying the backend and fit options. Default is None, which utilises the “minuit” backend with tol=0.1 and strategy=1.

Returns:
resultdict

Dictionary with the TS of the best fit value compared to the null hypothesis and fit results for the two hypotheses. Entries are:

  • “ts” : fit statistic difference with null hypothesis

  • “fit_results” : results for the best fit

  • “fit_results_null” : fit results for the null hypothesis

Examples

from gammapy.modeling.selection import select_nested_models
from gammapy.datasets import Datasets, SpectrumDatasetOnOff
from gammapy.modeling.models import SkyModel

# Test if cutoff is significant
dataset = SpectrumDatasetOnOff.read("$GAMMAPY_DATA/joint-crab/spectra/hess/pha_obs23523.fits")
datasets = Datasets(dataset)
model = SkyModel.create(spectral_model="ecpl", spatial_model="point", name='hess')
datasets.models = model
result = select_nested_models(datasets,
                              parameters=[model.spectral_model.lambda_],
                              null_values=[0],
                              )

# Test if source is significant
filename = "$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml"
filename_models = "$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_models.yaml"
fermi_datasets = Datasets.read(filename=filename, filename_models=filename_models)
model = fermi_datasets.models["Crab Nebula"]
# Number of parameters previously fit for the source of interest
n_free_parameters = len(model.parameters.free_parameters)
# Freeze spatial parameters to ensure another weaker source does not move from its position
# to replace the source of interest during the null hypothesis test.
# (with all parameters free you test N vs. N+1 models and not the detection of a specific source.)
fermi_datasets.models.freeze(model_type='spatial')
results = select_nested_models(fermi_datasets,
                              parameters=[model.spectral_model.amplitude],
                              null_values=[0],
                              n_free_parameters=n_free_parameters,
                              n_sigma=4,
                              )