{ "cells": [ { "cell_type": "markdown", "id": "4d8b5aaf", "metadata": {}, "source": [ "# Tutorial: fitting a BL Lac broad-band SED using MCMC fit with Gammapy\n", "\n", "In this tutorial, we will perform again the fit of the MWL SED of Mrk421 [illustrated in this notebook](ssc_gammapy_fit.ipynb). This time, instead of the simple $\\chi^2$ minimisation with `iminuit`, we will perform a fit with a Monte Carlo Markov Chain (MCMC) method. As for the case of the simple $\\chi^2$ minimisation, we will rely on `Gammapy`'s functionalities for data handling and fitting.\n", "\n", "Some `Gammapy` functions wrapping the MCMC package, [emcee](https://emcee.readthedocs.io/en/stable/), are available in the [gammapy-recipes](https://gammapy.github.io/gammapy-recipes/_build/html/index.html). Please download the [gammapy-recipes repository](https://github.com/gammapy/gammapy-recipes) to run the code in this notebook.\n", "\n", "An excellent theoretical introduction to MCMC methods can be found in [these notes](https://arxiv.org/abs/1909.12313)." ] }, { "cell_type": "code", "execution_count": 1, "id": "c39e0624", "metadata": {}, "outputs": [], "source": [ "# import numpy, astropy and matplotlib for basic functionalities\n", "import numpy as np\n", "import astropy.units as u\n", "import matplotlib.pyplot as plt\n", "import pkg_resources\n", "\n", "# import agnpy classes\n", "from agnpy.spectra import BrokenPowerLaw\n", "from agnpy.fit import SynchrotronSelfComptonModel, load_gammapy_flux_points\n", "from agnpy.utils.plot import load_mpl_rc, sed_y_label\n", "\n", "load_mpl_rc()\n", "\n", "# import gammapy classes\n", "from gammapy.modeling.models import SkyModel\n", "from gammapy.modeling import Fit" ] }, { "cell_type": "code", "execution_count": 2, "id": "6df4e5ea", "metadata": {}, "outputs": [], "source": [ "# the run_mcmc method has been moved from gammapy/gammapy to gammapy/recipes\n", "# please download the gammapy-recipes\n", "# https://github.com/gammapy/gammapy-recipes\n", "# and change the following path accordingly\n", "import sys\n", "\n", "sys.path.append(\"/Users/cosimo/work/gammapy-recipes/recipes/mcmc-sampling-emcee\")\n", "from sampling import (\n", " run_mcmc,\n", " par_to_model,\n", " plot_corner,\n", " plot_trace,\n", ")" ] }, { "cell_type": "markdown", "id": "c68c3e92", "metadata": {}, "source": [ "### `gammapy` wrapper of agnpy synchrotron and SSC\n", "Please refer to the [this notebook](ssc_gammapy_fit.ipynb) for a description of the `Gammapy` wrapper of the synchrotron and SSC radiative processes." ] }, { "cell_type": "code", "execution_count": 3, "id": "410bb604", "metadata": {}, "outputs": [], "source": [ "# electron energy distribution\n", "n_e = BrokenPowerLaw(\n", " k=1.3e-8 * u.Unit(\"cm-3\"),\n", " p1=2.02,\n", " p2=3.43,\n", " gamma_b=1e5,\n", " gamma_min=500,\n", " gamma_max=1e6,\n", ")\n", "\n", "# initialise the Gammapy SpectralModel\n", "ssc_model = SynchrotronSelfComptonModel(n_e, backend=\"gammapy\")" ] }, { "cell_type": "code", "execution_count": 4, "id": "c3547014", "metadata": {}, "outputs": [], "source": [ "# set intial parameters\n", "ssc_model.parameters[\"z\"].value = 0.0308\n", "ssc_model.parameters[\"delta_D\"].value = 18\n", "ssc_model.parameters[\"t_var\"].value = (1 * u.d).to_value(\"s\")\n", "ssc_model.parameters[\"t_var\"].frozen = True\n", "ssc_model.parameters[\"log10_B\"].value = -1.3" ] }, { "cell_type": "code", "execution_count": 5, "id": "162252dc", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
| type | name | value | unit | error | min | max | frozen | is_norm | link |
|---|---|---|---|---|---|---|---|---|---|
| str8 | str15 | float64 | str1 | int64 | float64 | float64 | bool | bool | str1 |
| spectral | log10_k | -7.8861e+00 | 0.000e+00 | -1.000e+01 | 1.000e+01 | False | False | ||
| spectral | p1 | 2.0200e+00 | 0.000e+00 | 1.000e+00 | 5.000e+00 | False | False | ||
| spectral | p2 | 3.4300e+00 | 0.000e+00 | 1.000e+00 | 5.000e+00 | False | False | ||
| spectral | log10_gamma_b | 5.0000e+00 | 0.000e+00 | 2.000e+00 | 6.000e+00 | False | False | ||
| spectral | log10_gamma_min | 2.6990e+00 | 0.000e+00 | 0.000e+00 | 4.000e+00 | True | False | ||
| spectral | log10_gamma_max | 6.0000e+00 | 0.000e+00 | 4.000e+00 | 8.000e+00 | True | False | ||
| spectral | z | 3.0800e-02 | 0.000e+00 | 1.000e-03 | 1.000e+01 | True | False | ||
| spectral | delta_D | 1.8000e+01 | 0.000e+00 | 1.000e+00 | 1.000e+02 | False | False | ||
| spectral | log10_B | -1.3000e+00 | 0.000e+00 | -4.000e+00 | 2.000e+00 | False | False | ||
| spectral | t_var | 8.6400e+04 | s | 0.000e+00 | 1.000e+01 | 3.142e+07 | True | False | |
| spectral | norm | 1.0000e+00 | 0.000e+00 | 1.000e-01 | 1.000e+01 | True | True |