Where to stick your Spectral Points?¶
The gammapy.spectrum module offers a number of options for positioning data points within an energy band. This example offers a comparison between
the log center and Lafferty & Wyatt (described in [Lafferty1994]) methods. See gammapy.spectrum.compute_differential_flux_points for documentation on usage.
Lafferty & Wyatt vs. Log Center¶
In measurements of spectra which do not vary rapidly, it is usually sufficient to represent the data within a bin using the log bin center. However, it is sometimes the case where the region of interest lies within a low frequency regime such that large bin widths must be chosen to ensure statistical errors remain small. In cases where the spectrum varies rapidly (for instance, exponentially) over wide bins, simply choosing the log bin center does not offer a good representation of the true underlying spectrum.
Instead, Lafferty & Wyatt [Lafferty1994] offer an alternative approach for such situations. In a bin of width \(\Delta E\) between bounds \(E_1\) and \(E_2\) for energy \(E\), the expectation \(<g_{meas}>\) of the true underlying spectrum, \(g(E)\) is defined as
As the bin size tends to zero, the expectation of the spectrum tends to it’s true value. The value of \(E\) within a bin for which the expectation should be regarded as a measurement of the true spectrum is determined by Lafferty & Wyatt as the energy at which the expectation value is equal to the mean value of the underlying true spectrum within that bin, noted as \(E_{lw}\). Thus knowledge of the true spectrum \(g(E)\) or an estimate for this (determined by fitting) is required.
So it follows that, in setting expectation equal to \(g(E)\) at \(E_{lw}\), the position of \(E_{lw}\) is given by the following equation:
For instances where a power law of spectral index 2 is taken, it can be analytically shown that the Lafferty & Wyatt method and log center method are coincident. In the case of steeper power laws (e.g. spectral index 3), the Lafferty & Wyatt method returns a lower coordinate on the energy axis than the log bin center, and the reverse effect is seen for less steep power laws (e.g. spectral index 1).
Power Law Assumption¶
In many “real world” examples, the nature of the true underlying spectrum is unknown and can only be estimated, either by using a fitting algorithm with data points or by assuming a certain spectral form. In this example, the true spectrum (being a piecewise power law function of spectral indices 1 through to 5 for bins of increasing energy) is shown in blue for illustration. This would be unknown to the user who, in this example, assumes the spectrum follows a simple power law with spectral index 4 across all energies.
The plot demonstrates that in these cases where the true spectrum is not known, the Lafferty & Wyatt method of positioning the data points offers a closer representation than the log center method. Residuals showing the percentage difference of each data point from the true spectrum are also shown.
"""Spectral plotting with gammapy.spectrum.flux_point
"""
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table
#from gammapy.spectrum.flux_point import (compute_differential_flux_points,
# _x_lafferty, _integrate)
#from gammapy.spectrum.powerlaw import power_law_integral_flux
SPECTRAL_INDEX = 4
def my_spectrum(x):
E_1 = 1
E_2 = 10
E_3 = 100
E_4 = 1000
E_5 = 10000
g1 = -1
g2 = -2
g3 = -3
g4 = -4
g5 = -5
im = np.select([x <= E_1, x <= E_2, x <= E_3, x <= E_4, x <= E_5, x > E_5],
[(x / E_1) ** g1, 1e-2 * (x / E_2) ** g2,
1e-5 * (x / E_3) ** g3, 1e-9 * (x / E_4) ** g4,
1e-14 * (x / E_5) ** g5, 0])
return im
def get_flux_tables(table, y_method, function, spectral_index):
table1 = table.copy()
lafferty_flux = compute_differential_flux_points('lafferty', y_method, table1,
function, spectral_index)
table2 = table1.copy()
log_flux = compute_differential_flux_points('log_center', y_method, table2,
function, spectral_index)
return lafferty_flux, log_flux
def plot_flux_points(table, x, y, function, energy_min, y_method):
f, ax = plt.subplots(2, sharex=True)
lafferty_flux, log_flux = get_flux_tables(table, y_method, function,
SPECTRAL_INDEX)
ax[0].loglog(lafferty_flux['ENERGY'],
((lafferty_flux['ENERGY'] ** 2) * lafferty_flux['DIFF_FLUX']),
marker='D', linewidth=0, ms=5,
label='Lafferty Method')
residuals_lafferty = (lafferty_flux['DIFF_FLUX']
- function(lafferty_flux['ENERGY'])) / function(lafferty_flux['ENERGY']) * 100
ax[0].loglog(log_flux['ENERGY'],
(log_flux['ENERGY'] ** 2) * log_flux['DIFF_FLUX'],
marker='D', linewidth=0, ms=5,
label='Log Center Method')
ax[0].legend(loc='lower left', fontsize=10)
residuals_log = (log_flux['DIFF_FLUX'] -
function(log_flux['ENERGY'])) / function(log_flux['ENERGY']) * 100
ax[1].semilogx(lafferty_flux['ENERGY'], residuals_lafferty, marker='D',
linewidth=0, ms=5)
ax[1].semilogx(log_flux['ENERGY'], residuals_log, marker='D',
linewidth=0, ms=5)
ax[0].loglog(x, (x ** 2) * y, color='k')
indices = np.arange(len(energy_min))
for index in indices:
ax[0].axvline(energy_min[index], 0, 1e6, color='k',
linestyle=':')
ax[1].axvline(energy_min[index], 0, 1e6, color='k',
linestyle=':')
ax[1].axhline(0, 0, 10, color='k')
ax[0].set_ylabel('E^2 * Differential Flux')
ax[1].set_ylabel('Residuals/%')
ax[1].set_xlabel('Energy')
ax[0].set_xlim([0.1, 10000])
ax[0].set_ylim([1e-6, 1e1])
return f
def plot_power_law():
# Define the function
x = np.arange(0.1, 100000, 0.1)
spectral_model = my_spectrum(x)
spectral_model_function = lambda x: my_spectrum(x)
y = spectral_model
# Set the x-bins
energy_min = [0.1, 1, 10, 100, 1000]
energy_max = [1, 10, 100, 1000, 10000]
energies = np.array(_x_lafferty(energy_min, energy_max,
spectral_model_function))
diff_fluxes = spectral_model_function(energies)
indices = np.array([0, 1, 2, 3, 4])
int_flux = power_law_integral_flux(diff_fluxes, (indices + 1),
energies, energy_min, energy_max)
special = lambda x: np.log(x)
int_flux[0] = np.abs(_integrate([energy_min[0]],
[energy_max[0]], special)[0])
# Put data into table
table = Table()
table['ENERGY_MIN'] = energy_min
table['ENERGY_MAX'] = energy_max
table['INT_FLUX'] = int_flux
plot_flux_points(table, x, y, spectral_model_function,
energy_min, 'power_law')
plt.tight_layout()
plt.legend()
Method Evaluation¶
The dependence of the positioning of the x coordinate of the flux point with assumed spectral index is demonstrated below. As before, the underlying spectrum is a simple power law of index 4. The position is given as the fraction of the size of the bin in Log Energy, for instance the Log Center method always positions the flux point at the center of the bin for a Log Energy scale, and so the position is 0.5. The coincidence of the two methods for the power law case of index 2 is evident, with the Lafferty method positioning the flux point to the right of the Log Center position for power law spectra of index less than 2, and to the left for index greater than 2.
"""Plot of x position depending on spectral index estimation
"""
import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table
#from gammapy.spectrum.flux_point import _energy_lafferty_power_law
from gammapy.spectrum.powerlaw import power_law_evaluate, power_law_integral_flux
from flux_point_demo import get_flux_tables
SPECTRAL_INDEX = 4
def make_x_plot():
energy_min = np.array([300])
energy_max = np.array([1000])
energies = np.array(_energy_lafferty_power_law(energy_min, energy_max,
SPECTRAL_INDEX))
diff_flux = power_law_evaluate(energies, 1, SPECTRAL_INDEX, 1)
# `True' differential & integral fluxes
int_flux = power_law_integral_flux(diff_flux, SPECTRAL_INDEX,
energies, energy_min, energy_max)
# Put data into table
table = Table()
table['ENERGY_MIN'] = energy_min
table['ENERGY_MAX'] = energy_max
table['INT_FLUX'] = int_flux
lafferty_array = []
log_array = []
spectral_indices = np.arange(1.1, 6, 0.01)
for spectral_index in spectral_indices:
lafferty_flux, log_flux = get_flux_tables(table, 'power_law', None,
spectral_index)
dlog_energy = np.log(energy_max) - np.log(energy_min)
residuals_lafferty = np.log(lafferty_flux['ENERGY'] - np.log(energy_min)) / dlog_energy
residuals_log = np.log(log_flux['ENERGY'] - np.log(energy_min)) / dlog_energy
lafferty_array.append(residuals_lafferty[0])
log_array.append(residuals_log[0])
plt.plot(spectral_indices, lafferty_array,
linewidth=1, ms=0, label='Lafferty Method')
plt.plot(spectral_indices, log_array,
linewidth=1, ms=0, label='Log Center Method')
plt.legend()
plt.ylabel('X position in bin')
plt.xlabel('Guessed spectral Index')
The ability of the two methods to represent an underlying spectrum is demonstrated in the residual images below. Here, flux points are positioned using the Lafferty method (left image) and Log Center method (central image) for different ‘true’ and ‘assumed’ power law spectral indices. Residuals are determined as the difference between the reconstructed flux compared to the true flux of the underlying spectrum. The closer the residual value is to zero (white), the better the representation offered.
The residuals log ratio image (right) shows the log ratio of residuals of the Lafferty method to those of the Log Center method, indicating the regions of assumed/true parameter space for which the Lafferty method offers an improved representation of the underlying frequency over the Log Center method (blue regions), and where it does not (red regions).
"""Create residuals image based on the two flux point methods
"""
import numpy as np
import matplotlib.pyplot as plt
#from gammapy.spectrum.flux_point import compute_differential_flux_points
from gammapy.spectrum.powerlaw import power_law_evaluate, power_law_integral_flux
def compute_flux_error(gamma_true, gamma_reco, method):
# Let's assume a concrete true spectrum and energy bin.
# Note that the residuals computed below do *not* depend on
# these parameters.
energy_min, energy_max = 1, 10
energy_ref, diff_flux_ref = 1, 1
# Compute integral flux in the energy band assuming `gamma_true`
int_flux = power_law_integral_flux(diff_flux_ref, gamma_true,
energy_ref, energy_min, energy_max)
# Compute flux point
table = compute_differential_flux_points(method, 'power_law',
spectral_index=gamma_reco,
energy_min=energy_min, energy_max=energy_max,
int_flux=int_flux)
# Compute relative error of the flux point
energy = table['ENERGY'].data
flux_reco = table['DIFF_FLUX'].data
flux_true = power_law_evaluate(energy, diff_flux_ref * np.ones_like(energy),
np.array(gamma_true).reshape(energy.shape),
energy_ref * np.ones_like(energy))
flux_true = flux_true.reshape(gamma_true.shape)
flux_reco = flux_reco.reshape(gamma_true.shape)
flux_error = (flux_reco - flux_true) / flux_true
return flux_error
def residuals_image():
fig = plt.figure(figsize=(15, 5))
gamma_true = np.arange(1.01, 7, 1)
gamma_reco = np.arange(1.01, 7, 1)
gamma_true, gamma_reco = np.meshgrid(gamma_true, gamma_reco)
flux_error_lafferty = compute_flux_error(gamma_true, gamma_reco,
method='lafferty')
flux_error_log_center = compute_flux_error(gamma_true, gamma_reco,
method='log_center')
flux_error_ratio = np.log10(flux_error_lafferty / flux_error_log_center)
extent = [0.5, 6.5, 0.5, 6.5]
vmin, vmax = -3, 3
axes_1 = fig.add_subplot(131)
axes_1.imshow(np.array(flux_error_lafferty),
interpolation='nearest', extent=extent,
origin="lower", vmin=vmin, vmax=vmax,
cmap='RdBu')
axes_1.set_ylabel('Assumed Spectral Index', fontsize=14)
axes_1.set_title('Lafferty Method', fontsize=12)
axes_2 = fig.add_subplot(132)
axes_2.imshow(np.array(flux_error_log_center),
interpolation='nearest', extent=extent,
origin="lower", vmin=vmin, vmax=vmax,
cmap='RdBu')
axes_2.set_xlabel('True Spectral Index', fontsize=14)
axes_2.set_title('Log-Center Method', fontsize=12)
axes_3 = fig.add_subplot(133)
im = axes_3.imshow(np.array(flux_error_ratio),
interpolation='nearest', extent=extent,
origin="lower", vmin=vmin, vmax=vmax,
cmap='RdBu')
axes_3.set_title('Residual Log Ratio: \n Log(Lafferty/Log Center)',
fontsize=12)
plt.tight_layout()
fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.92, 0.11, 0.025, 0.78])
fig.colorbar(im, cax=cbar_ax)
return fig
Caveats¶
It should be noted that in some cases (for instance, a power law with spectral index lower than 2), the Lafferty & Wyatt approach does not offer a better representation of the underlying spectrum compared to simply taking the log bin center. Indeed, the improvements that are seen in the Lafferty & Wyatt approach over using the bin Log Center occurs due to Lafferty & Wyatt simply positioning the flux point further to the left of the bin. Thus there is reason to consider further methods, namely:
- Positioning the flux point such that the integral flux on either side of it within the energy bin are equal.
- Simply choosing a position towards the left of the bin, for instance at a position 0.2 bin lengths from the left of the bin in log energy.