Source code for gammapy.stats.data

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""On-off bin stats computations."""
import numpy as np
from ..utils.random import get_random_state

__all__ = ["Stats", "make_stats", "combine_stats", "compute_total_stats"]


[docs]class Stats: """Container for an on-off observation. Parameters ---------- n_on : array_like Observed number of counts in the on region n_off : array_like Observed number of counts in the off region a_on : array_like Relative background exposure of the on region a_off : array_like Relative background exposure in the off region """ # TODO: use numpy arrays and properties # TODO: add gamma exposure def __init__(self, n_on, n_off, a_on, a_off): self.n_on = n_on self.n_off = n_off self.a_on = a_on self.a_off = a_off @property def alpha(self): r"""Background exposure ratio (float) .. math:: \alpha = a_\mathrm{on} / a_\mathrm{off} """ return self.a_on / self.a_off @property def background(self): r"""Background estimate (float) .. math:: \mu_\mathrm{bg} = \alpha\ n_\mathrm{off} """ return self.alpha * self.n_off @property def excess(self): r"""Excess (float) .. math:: n_\mathrm{ex} = n_\mathrm{on} - \mu_\mathrm{bg} """ return self.n_on - self.background def __str__(self): keys = ["n_on", "n_off", "a_on", "a_off", "alpha", "background", "excess"] values = [ self.n_on, self.n_off, self.a_on, self.a_off, self.alpha, self.background, self.excess, ] return "\n".join(["{} = {}".format(k, v) for (k, v) in zip(keys, values)])
[docs]def make_stats( signal, background, area_factor, weight_method="background", poisson_fluctuate=False, random_state="random-seed", ): """Fill using some weight method for the exposure. Parameters ---------- random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Returns ------- stats : `Stats` Statistics """ random_state = get_random_state(random_state) # Compute counts n_on = signal + background n_off = area_factor * background if poisson_fluctuate: n_on = random_state.poisson(n_on) n_off = random_state.poisson(n_off) # Compute weight if weight_method == "none": weight = 1 elif weight_method == "background": weight = background elif weight_method == "n_off": weight = n_off else: raise ValueError("Invalid weight_method: {}".format(weight_method)) # Compute exposure a_on = weight a_off = weight * area_factor return Stats(n_on, n_off, a_on, a_off)
[docs]def combine_stats(stats_1, stats_2, weight_method="none"): """Combine using some weight method for the exposure. Parameters ---------- stats_1 : `Stats` Observation 1 stats_2 : `Stats` Observation 2 weight_method : {'none', 'background', 'n_off'} Observation weighting method. Returns ------- stats : `Stats` Combined Observation 1 and 2 """ # Compute counts n_on = stats_1.n_on + stats_2.n_on n_off = stats_1.n_off + stats_2.n_off # Compute weights if weight_method == "none": weight_1 = 1 weight_2 = 1 elif weight_method == "background": weight_1 = stats_1.background() weight_2 = stats_2.background() elif weight_method == "n_off": weight_1 = stats_1.n_off weight_2 = stats_2.n_off else: raise ValueError("Invalid weight_method: {}".format(weight_method)) # Compute exposure a_on = weight_1 * stats_1.a_on + weight_2 * stats_2.a_on a_off = weight_1 * stats_1.a_off + weight_2 * stats_2.a_off return Stats(n_on, n_off, a_on, a_off)
[docs]def compute_total_stats(counts, exposure, background=None, solid_angle=None, mask=None): r"""Compute total stats for arrays of per-bin stats. The ``result`` dictionary contains a ``flux`` entry computed as .. math:: F = N / E = \sum{N_i} / \sum{E_i} as well as a ``flux3`` entry computed as .. math:: F = \sum{F_i} = \sum{\left(N_i / E_i\right)} where ``F`` is flux, ``N`` is excess and ``E`` is exposure. The output ``flux`` units are the inverse of the input ``exposure`` units, e.g. * ``exposure`` in cm^2 s -> ``flux`` in cm^-2 s^-1 * ``exposure`` in cm^2 s TeV -> ``flux`` in cm^-2 s^-1 TeV-1 The output ``surface_brightness`` units in addition depend on the ``solid_angle`` units, e.g. * ``exposure`` in cm^2 s and ``solid_angle`` in deg^2 -> ``surface_brightness`` in cm^-2 s^-1 deg^-2 TODOs: * integrate this with the `Stats` class. * add statistical errors on excess, flux, surface brightness Parameters ---------- counts, exposure : array_like Input arrays background, solid_angle, mask : array_like Optional input arrays Returns ------- result : dict Dictionary of total stats (for now, see the code for which entries it has). """ counts = np.asanyarray(counts) exposure = np.asanyarray(exposure) if solid_angle is None: background = np.zeros_like(counts) else: background = np.asanyarray(background) if solid_angle is None: solid_angle = np.ones_like(counts) else: solid_angle = np.asanyarray(solid_angle) if mask is None: mask = np.ones_like(counts) else: mask = np.asanyarray(mask) t = dict() t["n_pix_map"] = mask.size t["n_pix_mask"] = mask.sum() t["n_pix_fraction"] = t["n_pix_mask"] / float(t["n_pix_map"]) t["counts"] = counts[mask].sum(dtype=np.float64) t["background"] = background[mask].sum(dtype=np.float64) # Note that we use mean exposure (not sum) here!!! t["exposure"] = exposure[mask].mean(dtype=np.float64) t["solid_angle"] = solid_angle[mask].sum(dtype=np.float64) excess = counts - background t["excess"] = t["counts"] - t["background"] t["excess_2"] = excess[mask].sum(dtype=np.float64) flux = excess / exposure t["flux"] = (t["excess"]) / t["exposure"] t["flux_2"] = t["excess_2"] / t["exposure"] t["flux_3"] = flux[mask].sum(dtype=np.float64) surface_brightness = flux / solid_angle t["surface_brightness"] = t["flux"] / t["solid_angle"] t["surface_brightness_2"] = t["flux_2"] / t["solid_angle"] # Note that we use mean exposure (not sum) here!!! t["surface_brightness_3"] = surface_brightness[mask].mean(dtype=np.float64) return t