.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "tutorials/details/estimators.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_tutorials_details_estimators.py>`
        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_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 = 10








.. GENERATED FROM PYTHON SOURCE LINES 122-142

Note: The default scan range of the norm parameter is between 0.2 to 5. 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 142-146

.. code-block:: Python



    fp_result = fp_estimator.run(datasets=datasets)








.. GENERATED FROM PYTHON SOURCE LINES 147-149

Accessing and visualising the results
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. GENERATED FROM PYTHON SOURCE LINES 149-152

.. 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', 'norm_sensitivity', '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 153-155

We can specify the SED type to plot:


.. GENERATED FROM PYTHON SOURCE LINES 155-158

.. code-block:: Python

    fp_result.plot(sed_type="dnde")
    plt.show()




.. image-sg:: /tutorials/details/images/sphx_glr_estimators_001.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 159-163

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 163-171

.. 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/details/images/sphx_glr_estimators_002.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 172-174

The actual data members are N-dimensional `~gammapy.maps.region.ndmap.RegionNDMap` objects. So you can
also plot them:

.. GENERATED FROM PYTHON SOURCE LINES 174-177

.. code-block:: Python


    print(type(fp_result.dnde))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    <class 'gammapy.maps.region.ndmap.RegionNDMap'>




.. GENERATED FROM PYTHON SOURCE LINES 179-182

.. code-block:: Python

    fp_result.dnde.plot()
    plt.show()




.. image-sg:: /tutorials/details/images/sphx_glr_estimators_003.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 183-184

From the above, we can see that we access to many quantities.

.. GENERATED FROM PYTHON SOURCE LINES 187-188

Access the data:

.. GENERATED FROM PYTHON SOURCE LINES 188-191

.. 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]]

     [[-7.72411791e-28]]

     [[ 4.71387207e-12]]] TeV / (s cm2)




.. GENERATED FROM PYTHON SOURCE LINES 193-195

.. 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 197-199

.. 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 -4.87358893e-31  8.38258165e-16] 1 / (TeV s cm2)




.. GENERATED FROM PYTHON SOURCE LINES 200-201

Or even extract an energy range:

.. GENERATED FROM PYTHON SOURCE LINES 201-205

.. code-block:: Python


    fp_result.dnde.slice_by_idx({"energy": slice(3, 10)})







.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <pre>RegionNDMap

    	geom  : RegionGeom 
     	axes  : [&#x27;lon&#x27;, &#x27;lat&#x27;, &#x27;energy&#x27;]
    	shape : (1, 1, 5)
    	ndim  : 3
    	unit  : 1 / (TeV s cm2)
    	dtype : float64
    </pre>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 206-211

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 211-214

.. code-block:: Python


    print(fp_result.reference_model)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    SkyModel

      Name                      : 3b3UgCvW
      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 / (TeV s cm2)
        reference             (frozen):      1.000       TeV         






.. GENERATED FROM PYTHON SOURCE LINES 215-217

`~gammapy.estimators.FluxPoints` are the represented by the "norm" scaling factor with
respect to the reference model:

.. GENERATED FROM PYTHON SOURCE LINES 217-221

.. code-block:: Python


    fp_result.norm.plot()
    plt.show()




.. image-sg:: /tutorials/details/images/sphx_glr_estimators_004.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 222-236

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 236-241

.. code-block:: Python



    fp_result.counts.plot(axis_name="energy")
    plt.show()




.. image-sg:: /tutorials/details/images/sphx_glr_estimators_005.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 243-246

.. code-block:: Python

    fp_result.npred.plot(axis_name="energy")
    plt.show()




.. image-sg:: /tutorials/details/images/sphx_glr_estimators_006.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 248-251

.. code-block:: Python

    fp_result.npred_excess.plot(axis_name="energy")
    plt.show()




.. image-sg:: /tutorials/details/images/sphx_glr_estimators_007.png
   :alt: estimators
   :srcset: /tutorials/details/images/sphx_glr_estimators_007.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 252-257

Table conversion
~~~~~~~~~~~~~~~~

Flux points can be converted to tables:


.. GENERATED FROM PYTHON SOURCE LINES 257-261

.. 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        ... success  norm_scan 
           keV                keV                keV         ...                    
    ------------------ ------------------ ------------------ ... ------- -----------
     944060876.2859229  707945784.3841385 1258925411.7941697 ...    True 0.1 .. 10.0
      1778279410.03893 1258925411.7941697  2511886431.509587 ...    True 0.1 .. 10.0
    3349654391.5782814  2511886431.509587  4466835921.509632 ...    True 0.1 .. 10.0
     5956621435.290098  4466835921.509632 7943282347.2428255 ...    True 0.1 .. 10.0
    11220184543.019672 7943282347.2428255 15848931924.611168 ...    True 0.1 .. 10.0
    21134890398.366486 15848931924.611168 28183829312.644527 ...    True 0.1 .. 10.0
    39810717055.349655 28183829312.644527  56234132519.03493 ...    True 0.1 .. 10.0
     74989420933.24579  56234132519.03493 100000000000.00015 ...    True 0.1 .. 10.0




.. GENERATED FROM PYTHON SOURCE LINES 263-267

.. 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       ... success    norm_scan   
          keV             keV             keV        ...                        
    --------------- --------------- ---------------- ... ------- ---------------
      944060876.286   707945784.384   1258925411.794 ...    True 0.100 .. 10.000
     1778279410.039  1258925411.794   2511886431.510 ...    True 0.100 .. 10.000
     3349654391.578  2511886431.510   4466835921.510 ...    True 0.100 .. 10.000
     5956621435.290  4466835921.510   7943282347.243 ...    True 0.100 .. 10.000
    11220184543.020  7943282347.243  15848931924.611 ...    True 0.100 .. 10.000
    21134890398.366 15848931924.611  28183829312.645 ...    True 0.100 .. 10.000
    39810717055.350 28183829312.645  56234132519.035 ...    True 0.100 .. 10.000
    74989420933.246 56234132519.035 100000000000.000 ...    True 0.100 .. 10.000




.. GENERATED FROM PYTHON SOURCE LINES 268-275

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 277-280

Load the light curve for the PKS 2155-304 as a `~gammapy.estimators.FluxPoints` object.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 280-288

.. 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 289-292

Create a `~gammapy.estimators.FluxMaps` object through one of the estimators.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


.. GENERATED FROM PYTHON SOURCE LINES 292-298

.. 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                  : (np.int64(320), np.int64(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 300-301

.. 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']





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 22.044 seconds)


.. _sphx_glr_download_tutorials_details_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/v2.0?urlpath=lab/tree/notebooks/2.0/tutorials/details/estimators.ipynb
        :alt: Launch binder
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: estimators.ipynb <estimators.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: estimators.py <estimators.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: estimators.zip <estimators.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_