.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-time/Variability_estimation.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_analysis-time_Variability_estimation.py: Estimation of time variability in a lightcurve ============================================== Compute a series of time variability significance estimators for a lightcurve. Prerequisites ------------- Understanding the light curve estimator, please refer to the :doc:`/tutorials/analysis-time/light_curve` tutorial. For more in-depth explanation on the creation of smaller observations for exploring time variability, refer to the :doc:`/tutorials/analysis-time/light_curve_flare` tutorial. Context ------- Frequently, after computing a lightcurve, we need to quantify its variability in the time domain, for example in the case of a flare, burst, decaying light curve in GRBs or heightened activity in general. There are many ways to define the significance of the variability. **Objective: Estimate the level time variability in a lightcurve through different methods.** Proposed approach ----------------- We will start by reading the pre-computed light curve for PKS 2155-304 that is stored in `$GAMMAPY_DATA` To learn how to compute such an object, see the :doc:`/tutorials/analysis-time/light_curve_flare` tutorial. This tutorial will demonstrate how to compute different estimates which measure the significance of variability. These estimators range from basic ones that calculate the peak-to-trough variation, to more complex ones like fractional excess and point-to-point fractional variance, which consider the entire light curve. We also show an approach which utilises the change points in Bayesian blocks as indicators of variability. .. GENERATED FROM PYTHON SOURCE LINES 37-40 Setup ----- As usual, we’ll start with some general imports… .. GENERATED FROM PYTHON SOURCE LINES 40-54 .. code-block:: Python import numpy as np from astropy.stats import bayesian_blocks from astropy.time import Time import matplotlib.pyplot as plt from gammapy.estimators import FluxPoints from gammapy.estimators.utils import ( compute_lightcurve_doublingtime, compute_lightcurve_fpp, compute_lightcurve_fvar, ) from gammapy.maps import TimeMapAxis .. GENERATED FROM PYTHON SOURCE LINES 55-56 Load the light curve for the PKS 2155-304 flare directly from `$GAMMAPY_DATA/estimators`. .. GENERATED FROM PYTHON SOURCE LINES 56-65 .. code-block:: Python lc_1d = FluxPoints.read( "$GAMMAPY_DATA/estimators/pks2155_hess_lc/pks2155_hess_lc.fits", format="lightcurve" ) plt.figure(figsize=(8, 6)) lc_1d.plot(marker="o") plt.show() .. image-sg:: /tutorials/analysis-time/images/sphx_glr_Variability_estimation_001.png :alt: Variability estimation :srcset: /tutorials/analysis-time/images/sphx_glr_Variability_estimation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 66-80 Methods to characterize variability ----------------------------------- The three methods shown here are: - amplitude maximum variation - relative variability amplitude - variability amplitude. The amplitude maximum variation is the simplest method to define variability (as described in `Boller et al., 2016 `__) as it just computes the level of tension between the lowest and highest measured fluxes in the lightcurve. This estimator requires fully Gaussian errors. .. GENERATED FROM PYTHON SOURCE LINES 80-100 .. code-block:: Python flux = lc_1d.flux.quantity flux_err = lc_1d.flux_err.quantity f_mean = np.mean(flux) f_mean_err = np.mean(flux_err) f_max = flux.max() f_max_err = flux_err[flux.argmax()] f_min = flux.min() f_min_err = flux_err[flux.argmin()] amplitude_maximum_variation = (f_max - f_max_err) - (f_min + f_min_err) amplitude_maximum_significance = amplitude_maximum_variation / np.sqrt( f_max_err**2 + f_min_err**2 ) print(amplitude_maximum_significance) .. rst-class:: sphx-glr-script-out .. code-block:: none [[[12.41584196]]] .. GENERATED FROM PYTHON SOURCE LINES 101-104 There are other methods based on the peak-to-trough difference to assess the variability in a lightcurve. Here we present as example the relative variability amplitude as presented in `Kovalev et al., 2004 `__: .. GENERATED FROM PYTHON SOURCE LINES 104-119 .. code-block:: Python relative_variability_amplitude = (f_max - f_min) / (f_max + f_min) relative_variability_error = ( 2 * np.sqrt((f_max * f_min_err) ** 2 + (f_min * f_max_err) ** 2) / (f_max + f_min) ** 2 ) relative_variability_significance = ( relative_variability_amplitude / relative_variability_error ) print(relative_variability_significance) .. rst-class:: sphx-glr-script-out .. code-block:: none [[[19.613114]]] .. GENERATED FROM PYTHON SOURCE LINES 120-121 The variability amplitude as presented in `Heidt & Wagner, 1996 `__ is: .. GENERATED FROM PYTHON SOURCE LINES 121-143 .. code-block:: Python variability_amplitude = np.sqrt((f_max - f_min) ** 2 - 2 * f_mean_err**2) variability_amplitude_100 = 100 * variability_amplitude / f_mean variability_amplitude_error = ( 100 * ((f_max - f_min) / (f_mean * variability_amplitude_100 / 100)) * np.sqrt( (f_max_err / f_mean) ** 2 + (f_min_err / f_mean) ** 2 + ((np.std(flux, ddof=1) / np.sqrt(len(flux))) / (f_max - f_mean)) ** 2 * (variability_amplitude_100 / 100) ** 4 ) ) variability_amplitude_significance = ( variability_amplitude_100 / variability_amplitude_error ) print(variability_amplitude_significance) .. rst-class:: sphx-glr-script-out .. code-block:: none [[[6.12525306]]] .. GENERATED FROM PYTHON SOURCE LINES 144-151 Fractional excess variance, point-to-point fractional variance and doubling/halving time ----------------------------------------------------------------------------------------- The fractional excess variance, as presented by `Vaughan et al., 2003 `__, is a simple but effective method to assess the significance of a time variability feature in an object, for example an AGN flare. It is important to note that it requires Gaussian errors to be applicable. The excess variance computation is implemented in `~gammapy.estimators.utils`. .. GENERATED FROM PYTHON SOURCE LINES 151-155 .. code-block:: Python fvar_table = compute_lightcurve_fvar(lc_1d) print(fvar_table) .. rst-class:: sphx-glr-script-out .. code-block:: none min_energy max_energy ... fvar_err significance TeV TeV ... ------------------ ----------------- ... ------------------- ------------------ 0.5915030546513255 6.184989894219835 ... 0.01700709977114979 33.082691868487906 .. GENERATED FROM PYTHON SOURCE LINES 156-161 A similar estimator is the point-to-point fractional variance, as defined by `Edelson et al., 2002 `__, which samples the lightcurve on smaller time granularity. In general, the point-to-point fractional variance being higher than the fractional excess variance is indicative of the presence of very short timescale variability. .. GENERATED FROM PYTHON SOURCE LINES 161-165 .. code-block:: Python fpp_table = compute_lightcurve_fpp(lc_1d) print(fpp_table) .. rst-class:: sphx-glr-script-out .. code-block:: none min_energy max_energy ... fpp_err significance TeV TeV ... ------------------ ----------------- ... -------------------- ---------------- 0.5915030546513255 6.184989894219835 ... 0.017442925431194484 15.2995484265169 .. GENERATED FROM PYTHON SOURCE LINES 166-169 The characteristic doubling and halving time of the light curve, as introduced by `Brown, 2013 `__, can also be computed. This provides information on the shape of the variability feature, in particular how quickly it rises and falls. .. GENERATED FROM PYTHON SOURCE LINES 169-173 .. code-block:: Python dtime_table = compute_lightcurve_doublingtime(lc_1d, flux_quantity="flux") print(dtime_table) .. rst-class:: sphx-glr-script-out .. code-block:: none min_energy max_energy ... halving_err halving_coord TeV TeV ... s ------------------ ----------------- ... ------------------ ----------------- 0.5915030546513255 6.184989894219835 ... 20.935826709880043 53946.00422666667 .. GENERATED FROM PYTHON SOURCE LINES 174-185 Bayesian blocks --------------- The presence of temporal variability in a lightcurve can be assessed by using bayesian blocks (`Scargle et al., 2013 `__). A good and simple-to-use implementation of the algorithm is found in `astropy.stats.bayesian_blocks`. This implementation uses Gaussian statistics, as opposed to the `first introductory paper `__ which is based on Poissonian statistics. By passing the flux and error on the flux as ``measures`` to the method we can obtain the list of optimal bin edges defined by the bayesian blocks algorithm. .. GENERATED FROM PYTHON SOURCE LINES 185-192 .. code-block:: Python time = lc_1d.geom.axes["time"].time_mid.mjd bayesian_edges = bayesian_blocks( t=time, x=flux.flatten(), sigma=flux_err.flatten(), fitness="measures" ) .. GENERATED FROM PYTHON SOURCE LINES 193-195 The result giving a significance estimation for variability in the lightcurve is the number of *change points*, i.e. the number of internal bin edges: if at least one change point is identified by the algorithm, there is significant variability. .. GENERATED FROM PYTHON SOURCE LINES 195-199 .. code-block:: Python ncp = len(bayesian_edges) - 2 print(ncp) .. rst-class:: sphx-glr-script-out .. code-block:: none 7 .. GENERATED FROM PYTHON SOURCE LINES 200-203 We can rebin the lightcurve to compute the one expected with bayesian edges First, we adjust the first and last bins of the bayesian_edges to coincide with the original light curve start and end points. .. GENERATED FROM PYTHON SOURCE LINES 203-218 .. code-block:: Python # create a new axis axis_original = lc_1d.geom.axes["time"] bayesian_edges[0] = axis_original.time_edges[0].value bayesian_edges[-1] = axis_original.time_edges[-1].value edges = Time(bayesian_edges, format="mjd", scale=axis_original.reference_time.scale) axis_new = TimeMapAxis.from_time_edges(edges[:-1], edges[1:]) # rebin the lightcurve resample = lc_1d.resample_axis(axis_new) # plot the new lightcurve on top of the old one ax = lc_1d.plot(label="original") resample.plot(ax=ax, marker="s", label="rebinned") plt.legend() .. image-sg:: /tutorials/analysis-time/images/sphx_glr_Variability_estimation_002.png :alt: Variability estimation :srcset: /tutorials/analysis-time/images/sphx_glr_Variability_estimation_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. _sphx_glr_download_tutorials_analysis-time_Variability_estimation.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/analysis-time/Variability_estimation.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Variability_estimation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Variability_estimation.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_