.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/api/estimators.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_api_estimators.py: Estimators ========== This tutorial provides an overview of the `Estimator` API. All estimators live in the `gammapy.estimators` sub-module, offering a range of algorithms and classes for high-level flux and significance estimation. This is accomplished through a common functionality allowing the estimation of flux points, light curves, flux maps and profiles via a common API. Key Features ------------ - **Hypothesis Testing**: Estimations are based on testing a reference model against a null hypothesis, deriving flux and significance values. - **Estimation via Two Methods**: - **Model Fitting (Forward Folding)**: Refit the flux of a model component within specified energy, time, or spatial regions. - **Excess Calculation (Backward Folding)**: Use the analytical solution by Li and Ma for significance based on excess counts, currently available in `~gammapy.estimators.ExcessMapEstimator`. For further information on these details please refer to :doc:`/user-guide/estimators`. The setup --------- .. GENERATED FROM PYTHON SOURCE LINES 31-47 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from astropy import units as u from IPython.display import display from gammapy.datasets import SpectrumDatasetOnOff, Datasets, MapDataset from gammapy.estimators import ( FluxPointsEstimator, ExcessMapEstimator, FluxPoints, ) from gammapy.modeling import Fit, Parameter from gammapy.modeling.models import SkyModel, PowerLawSpectralModel from gammapy.utils.scripts import make_path .. GENERATED FROM PYTHON SOURCE LINES 48-55 Flux Points Estimation ---------------------- We start with a simple example for flux points estimation taking multiple datasets into account. In this section we show the steps to estimate the flux points. First we read the pre-computed datasets from `$GAMMAPY_DATA`. .. GENERATED FROM PYTHON SOURCE LINES 55-63 .. code-block:: Python datasets = Datasets() path = make_path("$GAMMAPY_DATA/joint-crab/spectra/hess/") for filename in path.glob("pha_obs*.fits"): dataset = SpectrumDatasetOnOff.read(filename) datasets.append(dataset) .. GENERATED FROM PYTHON SOURCE LINES 64-66 Next we define a spectral model and set it on the datasets: .. GENERATED FROM PYTHON SOURCE LINES 66-70 .. code-block:: Python pwl = PowerLawSpectralModel(index=2.7, amplitude="5e-11 cm-2 s-1 TeV-1") datasets.models = SkyModel(spectral_model=pwl, name="crab") .. GENERATED FROM PYTHON SOURCE LINES 71-75 Before using the estimators, it is necessary to first ensure that the model is properly fitted. This applies to all scenarios, including light curve estimation. To optimize the model parameters to best fit the data we utilise the following: .. GENERATED FROM PYTHON SOURCE LINES 75-80 .. code-block:: Python fit = Fit() fit_result = fit.optimize(datasets=datasets) print(fit_result) .. rst-class:: sphx-glr-script-out .. code-block:: none OptimizeResult backend : minuit method : migrad success : True message : Optimization terminated successfully. nfev : 37 total stat : 156.99 .. GENERATED FROM PYTHON SOURCE LINES 81-90 A fully configured Flux Points Estimation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The `~gammapy.estimators.FluxPointsEstimator` estimates flux points for a given list of datasets, energies and spectral model. The most simple way to call the estimator is by defining both the name of the ``source`` and its ``energy_edges``. Here we prepare a full configuration of the flux point estimation. Firstly we define the ``backend`` for the fit: .. GENERATED FROM PYTHON SOURCE LINES 90-96 .. code-block:: Python fit = Fit( optimize_opts={"backend": "minuit"}, confidence_opts={"backend": "scipy"}, ) .. GENERATED FROM PYTHON SOURCE LINES 97-99 Define the fully configured flux points estimator: .. GENERATED FROM PYTHON SOURCE LINES 99-113 .. code-block:: Python energy_edges = np.geomspace(0.7, 100, 9) * u.TeV norm = Parameter(name="norm", value=1.0) fp_estimator = FluxPointsEstimator( source="crab", energy_edges=energy_edges, n_sigma=1, n_sigma_ul=2, selection_optional="all", fit=fit, norm=norm, ) .. GENERATED FROM PYTHON SOURCE LINES 114-117 The ``norm`` parameter can be adjusted in a few different ways. For example, we can change its minimum and maximum values that it scans over, as follows. .. GENERATED FROM PYTHON SOURCE LINES 117-121 .. code-block:: Python fp_estimator.norm.scan_min = 0.1 fp_estimator.norm.scan_max = 3 .. GENERATED FROM PYTHON SOURCE LINES 122-142 Note: The default scan range of the norm parameter is between 0.1 to 10. In case the upper limit values lie outside this range, nan values will be returned. It may thus be useful to increase this range, specially for the computation of upper limits from weak sources. The various quantities utilised in this tutorial are described here: - ``source``: which source from the model to compute the flux points for - ``energy_edges``: edges of the flux points energy bins - ``n_sigma``: number of sigma for the flux error - ``n_sigma_ul``: the number of sigma for the flux upper limits - ``selection_optional``: what additional maps to compute - ``fit``: the fit instance (as defined above) - ``reoptimize``: whether to reoptimize the flux points with other model parameters, aside from the ``norm`` - ``norm``: normalisation parameter for the fit **Important note**: the output ``energy_edges`` are taken from the parent dataset energy bins, selecting the bins closest to the requested ``energy_edges``. To match the input bins directly, specific binning must be defined based on the parent dataset geometry. This could be done in the following way: `energy_edges = datasets[0].geoms["geom"].axes["energy"].downsample(factor=5).edges` .. GENERATED FROM PYTHON SOURCE LINES 145-147 .. code-block:: Python fp_result = fp_estimator.run(datasets=datasets) .. GENERATED FROM PYTHON SOURCE LINES 148-150 Accessing and visualising the results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 150-153 .. code-block:: Python print(fp_result) .. rst-class:: sphx-glr-script-out .. code-block:: none FluxPoints ---------- geom : RegionGeom axes : ['lon', 'lat', 'energy'] shape : (1, 1, 8) quantities : ['norm', 'norm_err', 'norm_errn', 'norm_errp', 'norm_ul', 'ts', 'npred', 'npred_excess', 'stat', 'stat_null', 'stat_scan', 'counts', 'success'] ref. model : pl n_sigma : 1 n_sigma_ul : 2 sqrt_ts_threshold_ul : 2 sed type init : likelihood .. GENERATED FROM PYTHON SOURCE LINES 154-156 We can specify the SED type to plot: .. GENERATED FROM PYTHON SOURCE LINES 156-159 .. code-block:: Python fp_result.plot(sed_type="dnde") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_001.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 160-164 We can also access the quantities names through ``fp_result.available_quantities``. Here we show how you can plot a different plot type and define the axes units, we also overlay the TS profile. .. GENERATED FROM PYTHON SOURCE LINES 164-172 .. code-block:: Python ax = plt.subplot() ax.xaxis.set_units(u.eV) ax.yaxis.set_units(u.Unit("TeV cm-2 s-1")) fp_result.plot(ax=ax, sed_type="e2dnde", color="tab:orange") fp_result.plot_ts_profiles(sed_type="e2dnde") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_002.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 173-175 The actual data members are N-dimensional `~gammapy.maps.region.ndmap.RegionNDMap` objects. So you can also plot them: .. GENERATED FROM PYTHON SOURCE LINES 175-178 .. code-block:: Python print(type(fp_result.dnde)) .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 180-183 .. code-block:: Python fp_result.dnde.plot() plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_003.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 184-185 From the above, we can see that we access to many quantities. .. GENERATED FROM PYTHON SOURCE LINES 188-189 Access the data: .. GENERATED FROM PYTHON SOURCE LINES 189-192 .. code-block:: Python print(fp_result.e2dnde.quantity.to("TeV cm-2 s-1")) .. rst-class:: sphx-glr-script-out .. code-block:: none [[[4.45517415e-11]] [[3.25823909e-11]] [[2.13944023e-11]] [[1.51958055e-11]] [[6.92052435e-12]] [[2.03447394e-12]] [[1.73951894e-18]] [[4.71387207e-12]]] TeV / (cm2 s) .. GENERATED FROM PYTHON SOURCE LINES 194-196 .. code-block:: Python print(fp_result.dnde.quantity.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (8, 1, 1) .. GENERATED FROM PYTHON SOURCE LINES 198-200 .. code-block:: Python print(fp_result.dnde.quantity[:, 0, 0]) .. rst-class:: sphx-glr-script-out .. code-block:: none [4.99878762e-11 1.03034567e-11 1.90677811e-12 4.28275990e-13 5.49716789e-14 4.55461982e-15 1.09756225e-21 8.38258165e-16] 1 / (cm2 s TeV) .. GENERATED FROM PYTHON SOURCE LINES 201-202 Or even extract an energy range: .. GENERATED FROM PYTHON SOURCE LINES 202-206 .. code-block:: Python fp_result.dnde.slice_by_idx({"energy": slice(3, 10)}) .. raw:: html 
RegionNDMap

    	geom  : RegionGeom 
     	axes  : ['lon', 'lat', 'energy']
    	shape : (1, 1, 5)
    	ndim  : 3
    	unit  : 1 / (cm2 s TeV)
    	dtype : float64


.. GENERATED FROM PYTHON SOURCE LINES 207-212 A note on the internal representation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The result contains a reference spectral model, which defines the spectral shape. Typically, it is the best fit model: .. GENERATED FROM PYTHON SOURCE LINES 212-215 .. code-block:: Python print(fp_result.reference_model) .. rst-class:: sphx-glr-script-out .. code-block:: none SkyModel Name : Z1cubNWV Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : Temporal model type : Parameters: index : 2.700 +/- 0.00 amplitude : 4.58e-11 +/- 0.0e+00 1 / (cm2 s TeV) reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 216-218 `~gammapy.estimators.FluxPoints` are the represented by the "norm" scaling factor with respect to the reference model: .. GENERATED FROM PYTHON SOURCE LINES 218-222 .. code-block:: Python fp_result.norm.plot() plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_004.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 223-237 Dataset specific quantities ("counts like") ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ While the flux estimate and associated errors are common to all datasets, the result also stores some dataset specific quantities, which can be useful for debugging. Here we remind the user of the meaning of the forthcoming quantities: - ``counts``: predicted counts from the null hypothesis, - ``npred``: predicted number of counts from best fit hypothesis, - ``npred_excess``: predicted number of excess counts from best fit hypothesis. The `~gammapy.maps.region.ndmap.RegionNDMap` allows for plotting of multidimensional data as well, by specifying the primary ``axis_name``: .. GENERATED FROM PYTHON SOURCE LINES 237-242 .. code-block:: Python fp_result.counts.plot(axis_name="energy") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_005.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 244-247 .. code-block:: Python fp_result.npred.plot(axis_name="energy") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_006.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 249-252 .. code-block:: Python fp_result.npred_excess.plot(axis_name="energy") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_estimators_007.png :alt: estimators :srcset: /tutorials/api/images/sphx_glr_estimators_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 253-258 Table conversion ~~~~~~~~~~~~~~~~ Flux points can be converted to tables: .. GENERATED FROM PYTHON SOURCE LINES 258-262 .. code-block:: Python table = fp_result.to_table(sed_type="flux", format="gadf-sed") display(table) .. rst-class:: sphx-glr-script-out .. code-block:: none e_ref e_min e_max flux flux_err ... is_ul counts success norm_scan keV keV keV keV / (cm2 s TeV) keV / (cm2 s TeV) ... ------------------ ------------------ ------------------ ---------------------- ---------------------- ... ----- ------------ ------- ---------- 944060876.2859229 707945784.3841385 1258925411.7941697 0.028262201290920164 0.0022082580047718126 ... False 64.0 .. 48.0 True 0.1 .. 3.0 1778279410.03893 1258925411.7941697 2511886431.509587 0.013396234895504788 0.0011199181643163176 ... False 39.0 .. 43.0 True 0.1 .. 3.0 3349654391.5782814 2511886431.509587 4466835921.509632 0.0038250882377506438 0.0005036517102677246 ... False 8.0 .. 22.0 True 0.1 .. 3.0 5956621435.290098 4466835921.509632 7943282347.2428255 0.0015277947990147364 0.0002825711636465431 ... False 2.0 .. 10.0 True 0.1 .. 3.0 11220184543.019672 7943282347.2428255 15848931924.611168 0.0004509608142981447 0.00015061726900511286 ... False 5.0 .. 2.0 True 0.1 .. 3.0 21134890398.366486 15848931924.611168 28183829312.644527 5.764921943853751e-05 5.292694386684509e-05 ... True 0.0 .. 1.0 True 0.1 .. 3.0 39810717055.349655 28183829312.644527 56234132519.03493 3.1946917901574355e-11 1.8815491773267503e-08 ... True 0.0 .. 0.0 True 0.1 .. 3.0 74989420933.24579 56234132519.03493 100000000000.00015 3.764602201409095e-05 3.764600555874051e-05 ... True 1.0 .. 0.0 True 0.1 .. 3.0 .. GENERATED FROM PYTHON SOURCE LINES 264-268 .. code-block:: Python table = fp_result.to_table(sed_type="likelihood", format="gadf-sed", formatted=True) display(table) .. rst-class:: sphx-glr-script-out .. code-block:: none e_ref e_min e_max ref_dnde ref_flux ... stat_scan is_ul counts success norm_scan keV keV keV 1 / (cm2 s TeV) keV / (cm2 s TeV) ... --------------- --------------- ---------------- --------------- ----------------- ... ------------------ ----- ------------ ------- -------------- 944060876.286 707945784.384 1258925411.794 5.349e-11 3.024e-02 ... 380.466 .. 381.242 False 64.0 .. 48.0 True 0.100 .. 3.000 1778279410.039 1258925411.794 2511886431.510 9.681e-12 1.259e-02 ... 352.415 .. 248.419 False 39.0 .. 43.0 True 0.100 .. 3.000 3349654391.578 2511886431.510 4466835921.510 1.752e-12 3.514e-03 ... 162.364 .. 126.794 False 8.0 .. 22.0 True 0.100 .. 3.000 5956621435.290 4466835921.510 7943282347.243 3.703e-13 1.321e-03 ... 99.705 .. 66.357 False 2.0 .. 10.0 True 0.100 .. 3.000 11220184543.020 7943282347.243 15848931924.611 6.702e-14 5.498e-04 ... 38.534 .. 49.779 False 5.0 .. 2.0 True 0.100 .. 3.000 21134890398.366 15848931924.611 28183829312.645 1.213e-14 1.535e-04 ... 8.144 .. 21.204 True 0.0 .. 1.0 True 0.100 .. 3.000 39810717055.350 28183829312.645 56234132519.035 2.195e-15 6.389e-05 ... 0.601 .. 10.884 True 0.0 .. 0.0 True 0.100 .. 3.000 74989420933.246 56234132519.035 100000000000.000 3.972e-16 1.784e-05 ... 7.361 .. 6.504 True 1.0 .. 0.0 True 0.100 .. 3.000 .. GENERATED FROM PYTHON SOURCE LINES 269-276 Common API ---------- In `GAMMAPY_DATA` we have access to other `~gammapy.estimators.FluxPoints` objects which have been created utilising the above method. Here we read the PKS 2155-304 light curve and create a `~gammapy.estimators.FluxMaps` object and show the data structure of such objects. We emphasize that these follow a very similar structure. .. GENERATED FROM PYTHON SOURCE LINES 278-281 Load the light curve for the PKS 2155-304 as a `~gammapy.estimators.FluxPoints` object. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 281-289 .. code-block:: Python lightcurve = FluxPoints.read( "$GAMMAPY_DATA/estimators/pks2155_hess_lc/pks2155_hess_lc.fits", format="lightcurve" ) display(lightcurve.available_quantities) .. rst-class:: sphx-glr-script-out .. code-block:: none ['norm', 'norm_err', 'norm_errn', 'norm_errp', 'norm_ul', 'ts', 'sqrt_ts', 'npred', 'npred_excess', 'stat', 'stat_null', 'stat_scan', 'is_ul', 'counts', 'success'] .. GENERATED FROM PYTHON SOURCE LINES 290-293 Create a `~gammapy.estimators.FluxMaps` object through one of the estimators. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 293-299 .. code-block:: Python dataset = MapDataset.read("$GAMMAPY_DATA/cta-1dc-gc/cta-1dc-gc.fits.gz") estimator = ExcessMapEstimator(correlation_radius="0.1 deg") result = estimator.run(dataset) display(result) .. rst-class:: sphx-glr-script-out .. code-block:: none FluxMaps -------- geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (320, 240, 1) quantities : ['npred', 'npred_excess', 'counts', 'ts', 'sqrt_ts', 'norm', 'norm_err'] ref. model : pl n_sigma : 1 n_sigma_ul : 2 sqrt_ts_threshold_ul : 2 sed type init : likelihood .. GENERATED FROM PYTHON SOURCE LINES 301-302 .. code-block:: Python display(result.available_quantities) .. rst-class:: sphx-glr-script-out .. code-block:: none ['npred', 'npred_excess', 'counts', 'ts', 'sqrt_ts', 'norm', 'norm_err'] .. _sphx_glr_download_tutorials_api_estimators.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/main?urlpath=lab/tree/notebooks/dev/tutorials/api/estimators.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: estimators.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: estimators.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: estimators.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_