Source code for gammapy.stats.data

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

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


[docs]class Stats(object): """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(['%s = %s' % (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 ------- """ 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: {0}".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: {0}".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). See also -------- gammapy.image.profile.FluxProfile.compute """ 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