# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Helper functions to work with distributions."""importnumbersimportnumpyasnpimportscipy.integratefromastropy.coordinatesimportAnglefromastropy.timeimportTimeDelta__all__=["draw","get_random_state","normalize","pdf","sample_powerlaw","sample_sphere","sample_sphere_distance","sample_times",]
[docs]defnormalize(func,x_min,x_max):"""Normalize a 1D function over a given range."""deff(x):returnfunc(x)/scipy.integrate.quad(func,x_min,x_max)[0]returnf
[docs]defpdf(func):"""One-dimensional PDF of a given radial surface density."""deff(x):returnx*func(x)returnf
[docs]defdraw(low,high,size,dist,random_state="random-seed",*args,**kwargs):"""Allow drawing of random numbers from any distribution."""from.inverse_cdfimportInverseCDFSamplern=1000x=np.linspace(low,high,n)pdf=dist(x)sampler=InverseCDFSampler(pdf=pdf,random_state=random_state)idx=sampler.sample(size)x_sampled=np.interp(idx,np.arange(n),x)returnnp.squeeze(x_sampled)
[docs]defget_random_state(init):"""Get a `numpy.random.RandomState` instance. The purpose of this utility function is to have a flexible way to initialise a `~numpy.random.RandomState` instance, a.k.a. a random number generator (``rng``). See :ref:`dev_random` for usage examples and further information. Parameters ---------- init : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Available options to initialise the RandomState object: * ``int`` -- new RandomState instance seeded with this integer (calls `~numpy.random.RandomState` with ``seed=init``) * ``'random-seed'`` -- new RandomState instance seeded in a random way (calls `~numpy.random.RandomState` with ``seed=None``) * ``'global-rng'``, return the RandomState singleton used by ``numpy.random``. * `~numpy.random.RandomState` -- do nothing, return the input. Returns ------- random_state : `~numpy.random.RandomState` RandomState instance. """ifisinstance(init,(numbers.Integral,np.integer)):returnnp.random.RandomState(init)elifinit=="random-seed":returnnp.random.RandomState(None)elifinit=="global-rng":returnnp.random.mtrand._randelifisinstance(init,np.random.RandomState):returninitelse:raiseValueError("{} cannot be used to seed a numpy.random.RandomState"" instance".format(init))
[docs]defsample_sphere(size,lon_range=None,lat_range=None,random_state="random-seed"):"""Sample random points on the sphere. Reference: http://mathworld.wolfram.com/SpherePointPicking.html Parameters ---------- size : int Number of samples to generate. lon_range : `~astropy.coordinates.Angle`, optional Longitude range (min, max). lat_range : `~astropy.coordinates.Angle`, optional Latitude range (min, max). random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Default is "random-seed". Returns ------- lon, lat: `~astropy.coordinates.Angle` Longitude and latitude coordinate arrays. """random_state=get_random_state(random_state)iflon_rangeisNone:lon_range=Angle([0.0,360.0],"deg")iflat_rangeisNone:lat_range=Angle([-90.0,90.0],"deg")# Sample random longitudeu=random_state.uniform(size=size)lon=lon_range[0]+(lon_range[1]-lon_range[0])*u# Sample random latitudev=random_state.uniform(size=size)z_range=np.sin(lat_range)z=z_range[0]+(z_range[1]-z_range[0])*vlat=np.arcsin(z)returnlon,lat
[docs]defsample_sphere_distance(distance_min=0,distance_max=1,size=None,random_state="random-seed"):"""Sample random distances if the 3-dim space density is constant. This function uses inverse transform sampling (`Wikipedia <http://en.wikipedia.org/wiki/Inverse_transform_sampling>`__) to generate random distances for an observer located in a 3-dim space with constant source density in the range ``(distance_min, distance_max)``. Parameters ---------- distance_min : float, optional Minimum distance. Default is 0. distance_max : float, optional Maximum distance. Default is 1. size : int or tuple of ints, optional Output shape. Default is one sample. Passed to `numpy.random.uniform`. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Default is "random-seed". Returns ------- distance : `~numpy.ndarray` Array of samples. """random_state=get_random_state(random_state)# Since the differential distribution is dP / dr ~ r ^ 2,# we have a cumulative distribution# P(r) = a * r ^ 3 + b# with P(r_min) = 0 and P(r_max) = 1 implying# a = 1 / (r_max ^ 3 - r_min ^ 3)# b = -a * r_min ** 3a=1.0/(distance_max**3-distance_min**3)b=-a*distance_min**3# Now for inverse transform sampling we need to use the inverse of# u = a * r ^ 3 + b# which is# r = [(u - b)/ a] ^ (1 / 3)u=random_state.uniform(size=size)distance=((u-b)/a)**(1.0/3)returndistance
[docs]defsample_powerlaw(x_min,x_max,gamma,size=None,random_state="random-seed"):"""Sample random values from a power law distribution. f(x) = x ** (-gamma) in the range x_min to x_max It is assumed that *gamma* is the **differential** spectral index. Reference: http://mathworld.wolfram.com/RandomNumber.html Parameters ---------- x_min : float Minimum x range. x_max : float Maximum x range. gamma : float Power law index. size : int, optional Number of samples to generate. Default is None. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Default is "random-seed". Returns ------- x : `~numpy.ndarray` Array of samples from the distribution. """random_state=get_random_state(random_state)size=int(size)exp=-gammabase=random_state.uniform(x_min**exp,x_max**exp,size)x=base**(1/exp)returnx
[docs]defsample_times(size,rate,dead_time=TimeDelta(0,format="sec"),return_diff=False,random_state="random-seed",):"""Make random times assuming a Poisson process. This function can be used to test event time series, to have a comparison what completely random data looks like. Can be used in two ways (in either case the return type is `~astropy.time.TimeDelta`): * ``return_delta=False`` - Return absolute times, relative to zero (default) * ``return_delta=True`` - Return time differences between consecutive events. Parameters ---------- size : int Number of samples. rate : `~astropy.units.Quantity` Event rate. dead_time : `~astropy.units.Quantity` or `~astropy.time.TimeDelta`, optional Dead time after event. return_diff : bool Return time difference between events. Default False, return absolute times. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Default is "random-seed". Returns ------- time : `~astropy.time.TimeDelta` Time differences (second) after time zero. Examples -------- Example how to simulate 100 events at a rate of 10 Hz. As expected the last event occurs after about 10 seconds. >>> from astropy.units import Quantity >>> from gammapy.utils.random import sample_times >>> rate = Quantity(10, 'Hz') >>> times = sample_times(size=100, rate=rate, random_state=0) >>> times[-1] <TimeDelta object: scale='None' format='sec' value=9.186484131475074> """random_state=get_random_state(random_state)dead_time=TimeDelta(dead_time)scale=(1/rate).to("s").valuetime_delta=random_state.exponential(scale=scale,size=size)time_delta+=dead_time.to("s").valueifreturn_diff:returnTimeDelta(time_delta,format="sec")else:time=time_delta.cumsum()returnTimeDelta(time,format="sec")