# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from ...utils.random import get_random_state
__all__ = ["GeneralRandom"]
[docs]class GeneralRandom(object):
"""Fast random number generation with an arbitrary pdf
of a continuous variable x.
Linear interpolation is applied between points pdf(x)
at which the pdf is specified.
I started with the recipy 576556, removed some unnecessary stuff
and added some useful stuff.
Recipe 576556: Generating random numbers with arbitrary distribution
http://code.activestate.com/recipes/576556/
Note: This class can only handle 1D distributions.
Note: Should it be required the cdf could be deleted
after computing to inversecdf to free memory since it is
not required for random number generation."""
def __init__(self, pdf, min_range, max_range, ninversecdf=None, ran_res=1e3):
"""Initialize the lookup table
Inputs:
x: random number values
pdf: probability density profile at that point
ninversecdf: number of reverse lookup values
Lookup is computed and stored in:
cdf: cumulative pdf
inversecdf: the inverse lookup table
delta_inversecdf: difference of inversecdf
ran_res: Resolution of the PDF
"""
self.ran_res = ran_res
x = np.linspace(min_range, max_range, ran_res)
# This is a good default for the number of reverse
# lookups to not loose much information in the pdf
if ninversecdf is None:
ninversecdf = 5 * x.size
self.nx = x.size
self.x = x
self.pdf = pdf(x)
self.cdf = np.empty(self.nx, dtype=float)
self.cdf[0] = 0
for i in range(1, self.nx):
temp = (self.pdf[i] + self.pdf[i - 1]) * (self.x[i] - self.x[i - 1]) / 2
self.cdf[i] = self.cdf[i - 1] + temp
self.pdf = self.pdf / self.cdf.max() # normalize pdf
self.cdf = self.cdf / self.cdf.max() # normalize cdf
self.ninversecdf = ninversecdf
y = np.arange(ninversecdf) / float(ninversecdf)
self.inversecdf = np.empty(ninversecdf)
self.inversecdf[0] = self.x[0]
cdf_idx = 0
for n in range(1, self.ninversecdf):
while self.cdf[cdf_idx] < y[n] and cdf_idx < ninversecdf:
cdf_idx += 1
self.inversecdf[n] = self.x[cdf_idx - 1] + (
self.x[cdf_idx] - self.x[cdf_idx - 1]
) * (y[n] - self.cdf[cdf_idx - 1]) / (
self.cdf[cdf_idx] - self.cdf[cdf_idx - 1]
)
if cdf_idx >= ninversecdf:
break
self.delta_inversecdf = np.concatenate((np.diff(self.inversecdf), [0]))
[docs] def draw(self, N=1000, random_state="random-seed"):
"""Returns an array of random numbers with the requested distribution.
The random numbers x are generated using the lookups
inversecdf and delta_inversecdf.
Parameters
----------
N : int
array length
random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
Defines random number generator initialisation.
Passed to `~gammapy.utils.random.get_random_state`.
Returns
-------
x : `~numpy.ndarray`
random numbers
"""
random_state = get_random_state(random_state)
# Generate uniform random float index in range [0, ninversecdf-1]
idx_f = random_state.uniform(size=N, high=self.ninversecdf - 1)
# Round down to next integer
idx = np.array(idx_f, "i")
# Use the inversecdf lookup to get the corresponding x
# and the delta_inversecdf lookup for linear interpolation
x = self.inversecdf[idx] + (idx_f - idx) * self.delta_inversecdf[idx]
return x
[docs] def make_plots(self, N=1e5):
"""Plot the pdf, cdf and inversecdf
and a random distribution of sample size N.
Useful for illustrating the interpolation and debugging."""
import matplotlib.pyplot as plt
# Plot the cdf
plt.figure()
plt.plot(self.x, self.cdf)
plt.title("cdf(x)")
# Plot the inverse cdf
plt.figure()
y = np.arange(self.ninversecdf) / float(self.ninversecdf)
plt.plot(y, self.inversecdf)
plt.title("inversecdf(y)")
# Plot the pdf and a random sample distribution
plt.figure()
x = self.draw(N)
# Use the same binning as self.x
binedges = self.x
plt.hist(x, bins=binedges, normed=True)
# x1 = 0.5*(edges[0:-1] + edges[1:])
# plot(x1, p1/p1.sum(),label='hist of random draws')
plt.plot(self.x, self.pdf, label="pdf")