Source code for gammapy.time.simulate

# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
from astropy.time import TimeDelta
from ..utils.random import get_random_state

__all__ = ["random_times"]


[docs]def random_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 (dimension: 1 / TIME) dead_time : `~astropy.units.Quantity` or `~astropy.time.TimeDelta`, optional Dead time after event (dimension: TIME) return_diff : bool Return time difference between events? (default: no, return absolute times) random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. 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.time import random_times >>> rate = Quantity(10, 'Hz') >>> times = random_times(size=100, rate=rate, random_state=0) >>> times[-1] <TimeDelta object: scale='None' format='sec' value=9.186484131475076> """ random_state = get_random_state(random_state) dead_time = TimeDelta(dead_time) scale = (1 / rate).to("s").value time_delta = random_state.exponential(scale=scale, size=size) time_delta += dead_time.to("s").value if return_diff: return TimeDelta(time_delta, format="sec") else: time = time_delta.cumsum() return TimeDelta(time, format="sec")