.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/details/parameter_ul.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_details_parameter_ul.py: Constraining parameter limits ============================= Explore how to deal with upper limits on parameters. Prerequisites ------------- It is advisable to understand the general Gammapy modelling and fitting framework before proceeding with this notebook, e.g. see :doc:`/user-guide/modeling`. Context ------- Even with significant detection of a source, constraining specific model parameters may remain difficult, allowing only for the calculation of confidence intervals. Proposed approach ----------------- In this section, we will use 6 observations of the blazar PKS 2155-304, taken in 2008 by H.E.S.S, to constrain the curvature in the spectrum. .. GENERATED FROM PYTHON SOURCE LINES 28-33 Setup ----- As usual, let’s start with some general imports… .. GENERATED FROM PYTHON SOURCE LINES 33-44 .. code-block:: Python # %matplotlib inline import matplotlib.pyplot as plt import numpy as np import astropy.units as u from gammapy.datasets import SpectrumDatasetOnOff, Datasets from gammapy.modeling import Fit, select_nested_models from gammapy.modeling.models import SkyModel, LogParabolaSpectralModel from gammapy.estimators import FluxPointsEstimator .. GENERATED FROM PYTHON SOURCE LINES 45-53 Load observation ---------------- We will use a `~gammapy.datasets.SpectrumDatasetOnOff` to see how to constrain model parameters. This dataset was obtained from H.E.S.S. observation of the blazar PKS 2155-304. Detailed modeling of this dataset can be found in the :doc:`/tutorials/astrophysics/ebl` notebook. .. GENERATED FROM PYTHON SOURCE LINES 53-61 .. code-block:: Python dataset_onoff = SpectrumDatasetOnOff.read( "$GAMMAPY_DATA/PKS2155-steady/pks2155-304_steady.fits.gz" ) dataset_onoff.peek() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_parameter_ul_001.png :alt: Counts, Exposure, Energy Dispersion :srcset: /tutorials/details/images/sphx_glr_parameter_ul_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 62-68 Fit spectrum ------------ We will investigate the presence of spectral curvature by modeling the observed spectrum using a `~gammapy.modeling.models.LogParabolaSpectralModel`. .. GENERATED FROM PYTHON SOURCE LINES 68-81 .. code-block:: Python spectral_model = LogParabolaSpectralModel( amplitude="5e-12 TeV-1 s-1 cm-2", alpha=2, beta=0.5, reference=1.0 * u.TeV ) model_pks = SkyModel(spectral_model, name="model_pks") dataset_onoff.models = model_pks fit = Fit() result_pks = fit.run(dataset_onoff) print(result_pks.models) .. rst-class:: sphx-glr-script-out .. code-block:: none DatasetModels Component 0: SkyModel Name : model_pks Datasets names : None Spectral model type : LogParabolaSpectralModel Spatial model type : Temporal model type : Parameters: amplitude : 4.23e-12 +/- 9.5e-13 1 / (TeV s cm2) reference (frozen): 1.000 TeV alpha : 3.351 +/- 0.35 beta : 0.263 +/- 0.54 .. GENERATED FROM PYTHON SOURCE LINES 82-87 We see that the parameter ``beta`` (the curvature parameter) is poorly constrained as the errors are very large. Therefore, we will perform a likelihood ratio test to evaluate the significance of the curvature compared to the null hypothesis of no curvature. In the null hypothesis, ``beta=0``. .. GENERATED FROM PYTHON SOURCE LINES 87-96 .. code-block:: Python LLR = select_nested_models( datasets=Datasets(dataset_onoff), parameters=[model_pks.parameters["beta"]], null_values=[0], ) print(LLR) .. rst-class:: sphx-glr-script-out .. code-block:: none {'ts': np.float64(0.32939559496572635), 'fit_results': , 'fit_results_null': } .. GENERATED FROM PYTHON SOURCE LINES 97-103 We can see that the improvement in the test statistic after including the curvature is only ~0.3, which corresponds to a significance of only 0.5. We can safely conclude that the addition of the curvature parameter does not significantly improve the fit. As a result, the function has internally updated the best fit model to the one corresponding to the null hypothesis (i.e. ``beta=0``). .. GENERATED FROM PYTHON SOURCE LINES 103-107 .. code-block:: Python print(dataset_onoff.models) .. rst-class:: sphx-glr-script-out .. code-block:: none DatasetModels Component 0: SkyModel Name : model_pks Datasets names : None Spectral model type : LogParabolaSpectralModel Spatial model type : Temporal model type : Parameters: amplitude : 3.85e-12 +/- 5.8e-13 1 / (TeV s cm2) reference (frozen): 1.000 TeV alpha : 3.405 +/- 0.30 beta (frozen): 0.000 .. GENERATED FROM PYTHON SOURCE LINES 108-115 Compute parameter asymmetric errors and upper limits ---------------------------------------------------- In such a case, it can still be useful to be able to constrain the allowed range of the non-significant parameter (e.g.: to rule out parameter values, to compare from theoretical predications, etc.). First, we reset the alternative model on the dataset: .. GENERATED FROM PYTHON SOURCE LINES 115-119 .. code-block:: Python dataset_onoff.models = LLR["fit_results"].models parameter = dataset_onoff.models.parameters["beta"] .. GENERATED FROM PYTHON SOURCE LINES 120-123 We can then compute the asymmetric errors and upper limits on the parameter of interest. It is always useful to ensure that the fit the converged by looking at the ``success`` and ``message`` keywords. .. GENERATED FROM PYTHON SOURCE LINES 123-127 .. code-block:: Python res_1sig = fit.confidence(datasets=dataset_onoff, parameter=parameter, sigma=1) print(res_1sig) .. rst-class:: sphx-glr-script-out .. code-block:: none {'success': True, 'message': 'Minos terminated successfully.', 'errp': np.float64(0.8310464830801719), 'errn': np.float64(0.41346440407107254), 'nfev': 126} .. GENERATED FROM PYTHON SOURCE LINES 128-129 We can directly use this to compute :math:`n\sigma` upper limits on the parameter: .. GENERATED FROM PYTHON SOURCE LINES 129-137 .. code-block:: Python res_2sig = fit.confidence(datasets=dataset_onoff, parameter=parameter, sigma=2) ll_2sigma = parameter.value - res_2sig["errn"] ul_2sigma = parameter.value + res_2sig["errp"] print(f"2-sigma lower limit on beta is {ll_2sigma:.2f}") print(f"2-sigma upper limit on beta is {ul_2sigma:.2f}") .. rst-class:: sphx-glr-script-out .. code-block:: none 2-sigma lower limit on beta is -0.43 2-sigma upper limit on beta is 3.71 .. GENERATED FROM PYTHON SOURCE LINES 138-143 Likelihood profile ------------------ We can also compute the likelihood profile of the parameter. First we define the scan range such that it encompasses more than the 2-sigma parameter limits. Then we call `~gammapy.modeling.Fit.stat_profile` : .. GENERATED FROM PYTHON SOURCE LINES 143-150 .. code-block:: Python parameter.scan_n_values = 25 parameter.scan_min = parameter.value - 2.5 * res_2sig["errn"] parameter.scan_max = parameter.value + 2.5 * res_2sig["errp"] parameter.interp = "lin" profile = fit.stat_profile(datasets=dataset_onoff, parameter=parameter, reoptimize=True) .. GENERATED FROM PYTHON SOURCE LINES 151-153 The resulting `profile` is a dictionary that stores the likelihood value and the fit result for each value of beta. .. GENERATED FROM PYTHON SOURCE LINES 153-157 .. code-block:: Python print(profile) .. rst-class:: sphx-glr-script-out .. code-block:: none {'model_pks.spectral.beta_scan': array([-1.4588051 , -1.02776659, -0.59672808, -0.16568956, 0.26534895, 0.69638746, 1.12742597, 1.55846448, 1.98950299, 2.4205415 , 2.85158001, 3.28261852, 3.71365703, 4.14469554, 4.57573405, 5.00677256, 5.43781107, 5.86884958, 6.29988809, 6.7309266 , 7.16196511, 7.59300362, 8.02404213, 8.45508064, 8.88611915]), 'stat_scan': array([50.67137298, 28.40382376, 13.59596178, 7.09723872, 6.00318251, 6.40313003, 7.05413881, 7.6675433 , 8.21483666, 8.71182744, 9.17165555, 9.6013471 , 10.00441683, 10.38274403, 10.73753108, 11.06974044, 11.38028286, 11.67009235, 11.94014857, 12.1914764 , 12.42513409, 12.64219749, 12.84374351, 13.03083462, 13.20450553]), 'fit_results': [, , , , , , , , , , , , , , , , , , , , , , , , ]} .. GENERATED FROM PYTHON SOURCE LINES 158-159 Let's plot everything together .. GENERATED FROM PYTHON SOURCE LINES 159-188 .. code-block:: Python values = profile["model_pks.spectral.beta_scan"] loglike = profile["stat_scan"] ax = plt.gca() ax.plot(values, loglike - np.min(loglike)) ax.set_xlabel("Beta") ax.set_ylabel(r"$\Delta$TS") ax.set_title(r"$\beta$-parameter likelihood profile") ax.fill_betweenx( x1=parameter.value - res_2sig["errn"], x2=parameter.value + res_2sig["errp"], y=[-0.5, 25], alpha=0.3, color="pink", label="1-sigma range", ) ax.fill_betweenx( x1=parameter.value - res_1sig["errn"], x2=parameter.value + res_1sig["errp"], y=[-0.5, 25], alpha=0.3, color="salmon", label="2-sigma range", ) ax.set_ylim(-0.5, 25) plt.legend() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_parameter_ul_002.png :alt: $\beta$-parameter likelihood profile :srcset: /tutorials/details/images/sphx_glr_parameter_ul_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 190-196 Impact of the model choice on the flux upper limits --------------------------------------------------- The flux points depends on the underlying model assumption. This can have a non-negligible impact on the flux upper limits in the energy range where the model is not well constrained as illustrated in the following figure. So quote preferably upper limits from the model which is the most supported by the data. .. GENERATED FROM PYTHON SOURCE LINES 196-208 .. code-block:: Python energies = dataset_onoff.geoms["geom"].axes["energy"].edges fpe = FluxPointsEstimator(energy_edges=energies, n_jobs=4, selection_optional=["ul"]) # Null hypothesis -- no curvature dataset_onoff.models = LLR["fit_results_null"].models fp_null = fpe.run(dataset_onoff) # Alternative hypothesis -- with curvature dataset_onoff.models = LLR["fit_results"].models fp_alt = fpe.run(dataset_onoff) .. GENERATED FROM PYTHON SOURCE LINES 209-210 Plot them together .. GENERATED FROM PYTHON SOURCE LINES 210-247 .. code-block:: Python ax = fp_null.plot(sed_type="e2dnde", color="blue") LLR["fit_results_null"].models[0].spectral_model.plot( ax=ax, energy_bounds=(energies[0], energies[-1]), sed_type="e2dnde", color="blue", label="No curvature", ) LLR["fit_results_null"].models[0].spectral_model.plot_error( ax=ax, energy_bounds=(energies[0], energies[-1]), sed_type="e2dnde", facecolor="blue", alpha=0.2, ) fp_alt.plot(ax=ax, sed_type="e2dnde", color="red") LLR["fit_results"].models[0].spectral_model.plot( ax=ax, energy_bounds=(energies[0], energies[-1]), sed_type="e2dnde", color="red", label="with curvature", ) LLR["fit_results"].models[0].spectral_model.plot_error( ax=ax, energy_bounds=(energies[0], energies[-1]), sed_type="e2dnde", facecolor="red", alpha=0.2, ) plt.legend() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_parameter_ul_003.png :alt: parameter ul :srcset: /tutorials/details/images/sphx_glr_parameter_ul_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 248-252 This logic can be extended to any spectral or spatial feature. As an exercise, try to compute the 95% spatial extent on the MSH 15-52 dataset used for the ring background notebook. .. _sphx_glr_download_tutorials_details_parameter_ul.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v2.0?urlpath=lab/tree/notebooks/2.0/tutorials/details/parameter_ul.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: parameter_ul.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: parameter_ul.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: parameter_ul.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_