This module offers various stochastic generators for point processes that can be used as spike trains.
Create an StGen object:
>>> st_gen = StGen()
This will initialize the stochastic generator and by default try to create a numpy random generator instance.
Optionally, you can also pass a random number generator instance to the constructor:
>>> import numpy
>>> st_gen = StGen(rng = numpy.random.RandomState())
You can also use random number generators from gnu scientific library (gsl):
>>> from pygsl.rng import rng
>>> st_gen_gsl = StGen(rng = rng())
If you want to seed the random number generator with a specific seed, you can do so in the constructor:
>>> st_gen = StGen(seed = 1234567)
Alternatively, you can re-seed the random number generator when the StGen object has already been created:
>>> st_gen.seed(7654321)
Using the StGen-object, you can generate point processes with inter-spike-intervals distributed according to a poisson distribution:
>>> st_gen = StGen()
>>> spike_train_poisson = st_gen.poisson_generator(rate = 100.,
tstart = 0.,
tstop = 2500.)
This generates a NeuroTools.SpikeTrain object, containing spike times with an approximate rate of 100 Hz and a duration of 2.5 seconds.
If you want a numpy array of spike times rather than a SpikeTrain object, specify the array keyword:
>>> spike_train_array = st_gen.poisson_generator(rate = 100., array = True)
StGen can also generate inhomogeneous poisson processes, i.e. spike trains with dynamically changing rates:
>>> spike_train_dyn = st_gen.poissondyn_generator(rate = [50., 80., 30.],
t = [0., 1000., 2000.],
tstop = 2.5,
array = False)
This will generate a SpikeTrain object containing spike times with an approximate rate of 50 Hz for one second, followed by 80 Hz for one second, and finally 30 Hz for half a second. Note that t[0] is used as tstart.
A collection of tools for stochastic process generation.
shotnoise_fromspikes - Convolves the provided spike train with shot decaying exponential.
gamma_hazard - Compute the hazard function for a gamma process with parameters a,b.
Generates an Orstein Ulbeck process using the forward euler method. The function returns an AnalogSignal object.
dt - the time resolution in milliseconds of th signal tau - the correlation time in milliseconds sigma - std dev of the process y0 - initial value of the process, at t_start t_start - start time in milliseconds t_stop - end time in milliseconds array - if True, the functions returns the tuple (y,t)
where y and t are the OU signal and the time bins, respectively, and are both numpy arrays.
>> stgen.OU_generator(0.1, 2, 3, 0, 0, 10000)
Generates an Orstein Ulbeck process using the forward euler method. The function returns an AnalogSignal object.
OU_generator_weave1, as opposed to OU_generator, uses scipy.weave and is thus much faster.
dt - the time resolution in milliseconds of th signal tau - the correlation time in milliseconds sigma - std dev of the process y0 - initial value of the process, at t_start t_start - start time in milliseconds t_stop - end time in milliseconds array - if True, the functions returns the tuple (y,t)
where y and t are the OU signal and the time bins, respectively, and are both numpy arrays.
>> stgen.OU_generator_weave1(0.1, 2, 3, 0, 0, 10000)
OU_generator
Returns a SpikeTrain whose spikes are a realization of a gamma process with the given shape a, b and stopping time t_stop (milliseconds). (average rate will be a*b)
Note: t_start is always 0.0, thus all realizations are as if they spiked at t=0.0, though this spike is not included in the SpikeList.
a,b - the parameters of the gamma process t_start - the beginning of the SpikeTrain (in ms) t_stop - the end of the SpikeTrain (in ms) array - if True, a numpy array of sorted spikes is returned,
rather than a SpikeTrain object.
>> gen.gamma_generator(10, 1/10., 0, 1000) >> gen.gamma_generator(20, 1/5., 5000, 10000, array=True)
inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
Returns a SpikeList whose spikes are an inhomogeneous realization (dynamic rate) of the so-called 2D adapting markov process (see references). 2D implies the process has two states, an adaptation state, and a refractory state, both of which affect its probability to spike. The implementation uses the thinning method, as presented in the references.
For the 1d implementation, with no relative refractoriness, see the inh_adaptingmarkov_generator.
- a,bq - arrays of the parameters of the hazard function where a[i] and bq[i]
- will be active on interval [t[i],t[i+1]]
tau_s - the time constant of adaptation (in milliseconds). tau_r - the time constant of refractoriness (in milliseconds). qrqs - the ratio of refractoriness conductance to adaptation conductance.
typically on the order of 200.
- t - an array specifying the time bins (in milliseconds) at which to
- specify the rate
t_stop - length of time to simulate process (in ms) array - if True, a numpy array of sorted spikes is returned,
rather than a SpikeList object.
- t_start=t[0]
- a is in units of Hz. Typical values are available in Fig. 1 of Muller et al 2007, a~5-80Hz (low to high stimulus)
- bq here is taken to be the quantity b*q_s in Muller et al 2007, is thus dimensionless, and has typical values bq~3.0-1.0 (low to high stimulus)
- qrqs is the quantity q_r/q_s in Muller et al 2007, where a value of qrqs = 3124.0nS/14.48nS = 221.96 was used.
- tau_s has typical values on the order of 100 ms
- tau_r has typical values on the order of 2 ms
Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories Neural Comput. 2007 19: 2958-3010.
Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
See source:trunk/examples/stgen/inh_2Dmarkov_psth.py
inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
Returns a SpikeList whose spikes are an inhomogeneous realization (dynamic rate) of the so-called adapting markov process (see references). The implementation uses the thinning method, as presented in the references.
This is the 1d implementation, with no relative refractoriness. For the 2d implementation with relative refractoriness, see the inh_2dadaptingmarkov_generator.
- a,bq - arrays of the parameters of the hazard function where a[i] and bq[i]
- will be active on interval [t[i],t[i+1]]
tau - the time constant of adaptation (in milliseconds). t - an array specifying the time bins (in milliseconds) at which to
specify the ratet_stop - length of time to simulate process (in ms) array - if True, a numpy array of sorted spikes is returned,
rather than a SpikeList object.
- t_start=t[0]
- a is in units of Hz. Typical values are available in Fig. 1 of Muller et al 2007, a~5-80Hz (low to high stimulus)
- bq here is taken to be the quantity b*q_s in Muller et al 2007, is thus dimensionless, and has typical values bq~3.0-1.0 (low to high stimulus)
- tau_s has typical values on the order of 100 ms
Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories Neural Comput. 2007 19: 2958-3010.
Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
See source:trunk/examples/stgen/inh_2Dmarkov_psth.py
inh_poisson_generator, inh_gamma_generator, inh_2dadaptingmarkov_generator
Returns a SpikeList whose spikes are a realization of an inhomogeneous gamma process (dynamic rate). The implementation uses the thinning method, as presented in the references.
- a,b - arrays of the parameters of the gamma PDF where a[i] and b[i]
- will be active on interval [t[i],t[i+1]]
- t - an array specifying the time bins (in milliseconds) at which to
- specify the rate
t_stop - length of time to simulate process (in ms) array - if True, a numpy array of sorted spikes is returned,
rather than a SpikeList object.
t_start=t[0] a is a dimensionless quantity > 0, but typically on the order of 2-10. a = 1 results in a poisson process. b is assumed to be in units of 1/Hz (seconds).
Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories Neural Comput. 2007 19: 2958-3010.
Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
See source:trunk/examples/stgen/inh_gamma_psth.py
inh_poisson_generator, gamma_hazard
Returns a SpikeTrain whose spikes are a realization of an inhomogeneous poisson process (dynamic rate). The implementation uses the thinning method, as presented in the references.
- rate - an array of the rates (Hz) where rate[i] is active on interval
- [t[i],t[i+1]]
- t - an array specifying the time bins (in milliseconds) at which to
- specify the rate
t_stop - length of time to simulate process (in ms) array - if True, a numpy array of sorted spikes is returned,
rather than a SpikeList object.
t_start=t[0]
Eilif Muller, Lars Buesing, Johannes Schemmel, and Karlheinz Meier Spike-Frequency Adapting Neural Ensembles: Beyond Mean Adaptation and Renewal Theories Neural Comput. 2007 19: 2958-3010.
Devroye, L. (1986). Non-uniform random variate generation. New York: Springer-Verlag.
>> time = arange(0,1000) >> stgen.inh_poisson_generator(time,sin(time), 1000)
poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
Returns a SpikeTrain whose spikes are a realization of a Poisson process with the given rate (Hz) and stopping time t_stop (milliseconds).
Note: t_start is always 0.0, thus all realizations are as if they spiked at t=0.0, though this spike is not included in the SpikeList.
rate - the rate of the discharge (in Hz) t_start - the beginning of the SpikeTrain (in ms) t_stop - the end of the SpikeTrain (in ms) array - if True, a numpy array of sorted spikes is returned,
rather than a SpikeTrain object.
>> gen.poisson_generator(50, 0, 1000) >> gen.poisson_generator(20, 5000, 10000, array=True)
inh_poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator
Compute the hazard function for a gamma process with parameters a,b where a and b are the parameters of the gamma PDF: y(t) = x^(a-1) exp(-x/b) / (Gamma(a)*b^a)
Compute the hazard function for a gamma process with parameters a,b where a and b are the parameters of the gamma PDF: y(t) = x^(a-1) exp(-x/b) / (Gamma(a)*b^a)
Convolves the provided spike train with shot decaying exponentials yielding so called shot noise if the spike train is Poisson-like. Returns an AnalogSignal if array=False, otherwise (shotnoise,t) as numpy arrays.
spike_train - a SpikeTrain object q - the shot jump for each spike tau - the shot decay time constant in milliseconds dt - the resolution of the resulting shotnoise in milliseconds t_start - start time of the resulting AnalogSignal
If unspecified, t_start of spike_train is used
- t_stop - stop time of the resulting AnalogSignal
- If unspecified, t_stop of spike_train is used
array - if True, returns (shotnoise,t) as numpy arrays, otherwise an AnalogSignal. eps - a numerical parameter indicating at what value of the shot kernal the tail is cut. The default is usually fine.
Spikes in spike_train before t_start are taken into account in the convolution.
>> stg = stgen.StGen() >> st = stg.poisson_generator(10.0,0.0,1000.0) >> g_e = shotnoise_fromspikes(st,2.0,10.0,dt=0.1)
poisson_generator, inh_gamma_generator, inh_adaptingmarkov_generator, OU_generator ...