Source code for gammapy.utils.random.inverse_cdf

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import html
import numpy as np
from .utils import get_random_state

__all__ = ["InverseCDFSampler"]


[docs] class InverseCDFSampler: """Inverse CDF sampler. It determines a set of random numbers and calculate the cumulative distribution function. Parameters ---------- pdf : `~gammapy.maps.Map` Map of the predicted source counts. axis : int Axis along which sampling the indexes. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. """ def __init__(self, pdf, axis=None, random_state=0): self.random_state = get_random_state(random_state) self.axis = axis if axis is not None: self.cdf = np.cumsum(pdf, axis=self.axis) self.cdf /= self.cdf[:, [-1]] else: self.pdf_shape = pdf.shape pdf = pdf.ravel() / pdf.sum() self.sortindex = np.argsort(pdf, axis=None) self.pdf = pdf[self.sortindex] self.cdf = np.cumsum(self.pdf) def _repr_html_(self): try: return self.to_html() except AttributeError: return f"<pre>{html.escape(str(self))}</pre>"
[docs] def sample_axis(self): """Sample along a given axis. Returns ------- index : tuple of `~numpy.ndarray` Coordinates of the drawn sample. """ choices = self.random_state.uniform(high=1, size=len(self.cdf)) shape_cdf = self.cdf.shape cdf_all = np.insert(self.cdf, 0, 0, axis=1) edges = np.arange(shape_cdf[1] + 1) - 0.5 pix_coords = [] for cdf, choice in zip(cdf_all, choices): pix = np.interp(choice, cdf, edges) pix_coords.append(pix) return np.array(pix_coords)
[docs] def sample(self, size): """Draw sample from the given PDF. Parameters ---------- size : int Number of samples to draw. Returns ------- index : tuple of `~numpy.ndarray` Coordinates of the drawn sample. """ # pick numbers which are uniformly random over the cumulative distribution function choice = self.random_state.uniform(high=1, size=size) # find the indices corresponding to this point on the CDF index = np.searchsorted(self.cdf, choice) index = self.sortindex[index] # map back to multi-dimensional indexing index = np.unravel_index(index, self.pdf_shape) index = np.vstack(index) index = index + self.random_state.uniform(low=-0.5, high=0.5, size=index.shape) return index