# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install --upgrade iqplot watermark"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
# ------------------------------We have previously modeled spiking using a noisy leaky integrate-and-fire model. We saw that in some limits, the noisy LIF model gives Poissonian spiking such that the interspike intervals are Exponentially distributed and the number of spikes in a given time interval are Poisson distributed. Poisson process-based models of spiking are widely used to good effect.
However, though well-modelled as Poisson, the rate of spiking can vary over time in what we call a nonhomogeneous Poisson process. Here, we derive the PDF for ISIs in this context and show how to sample ISIs.
Consider the waiting time for an arrival of a Poisson process that arrives at rate \(\beta\). We know that the waiting time is Exponentially distributed.
\[\begin{align} f(t\mid \beta) = \beta \mathrm{e}^{-\beta t}. \end{align} \tag{16.1}\]
The probability that a Poisson process has not arrived at or before time \(t\) is given by the complementary cumulative distribution function (CCDF) of the Exponential distribution.
\[\begin{align} P(\text{not arrived by time }t) = 1 - P(\text{arrived by time }t) = 1 - \int_0^t\mathrm{d}t'\,f(t\mid\beta) = 1 - (1 - \mathrm{e}^{-\beta t}) = \mathrm{e}^{-\beta t}. \end{align} \tag{16.2}\]
Now, consider the case where the arrival rate varies with time; \(\beta = \beta(t)\). Consider some small interval, \([0, \Delta t]\). The probability that no Poisson process has arrived in time \(\Delta t\) is
\[\begin{align} P(\text{not arrived by time } \Delta t) \approx \mathrm{e}^{-\beta(\Delta t) \Delta t}, \end{align} \tag{16.3}\]
where \(\Delta t\) is small so that \(\beta(t)\) does not change appreciably over the time interval and \(\beta(t) \approx \beta(t+\Delta t)\). We will soon take the \(\Delta t\to 0\) limit, so this approximation is justified.
Now consider another interval of length \(\Delta t\) right after the one we just considered. The probability that no process arrives in that interval is
\[\begin{align} P(\text{not arrived in time window } [\Delta t, 2\Delta t]) \approx \mathrm{e}^{-\beta(2\Delta t) \Delta t}. \end{align} \tag{16.4}\]
Therefore, owing to the memorylessness of Poisson processes, the probability that no processes arrived in the interval \([0, 2\Delta t]\) is the product of the probabilities of not arriving in the respective subintervals,
\[\begin{align} P(\text{not arrived in time } 2\Delta t) \approx \mathrm{e}^{-\beta(\Delta t) \Delta t}\,\mathrm{e}^{-\beta(2\Delta t) \Delta t}. \end{align} \tag{16.5}\]
We can consider \(m\) such intervals.
\[\begin{align} P(\text{not arrived between in time } m\Delta t) \approx \prod_{k=1}^m \mathrm{e}^{-\beta(k\Delta t) \Delta t} = \exp\left[-\sum_{k=1}^m \Delta t \,\beta(k\Delta t)\right]. \end{align} \tag{16.6}\]
In the limit of \(\Delta t\to 0\), the Riemann sum becomes an integral. Taking \(t = m\Delta t\), we have
\[\begin{align} P(\text{not arrived by time } t) = \exp\left[- \int_0^t\mathrm{d}t'\,\beta(t')\right]. \end{align} \tag{16.7}\]
This is the CCDF of the distribution desribing the arrival of a Poisson process with variable rate \(\beta(t)\). The CDF is then \(1 - \text{CCDF}\),
\[\begin{align} F(t\mid \beta(t)) = 1 - \exp\left[- \int_0^t\mathrm{d}t'\,\beta(t')\right]. \end{align} \tag{16.8}\]
The probability density function is then
\[\begin{align} f(t\mid \beta(t)) &= \frac{\mathrm{d}F}{\mathrm{d}t} = -\frac{\mathrm{d}}{\mathrm{d}t}\,\exp\left[- \int_0^t\mathrm{d}t'\,\beta(t')\right] = -\exp\left[- \int_0^t\mathrm{d}t'\,\beta(t')\right]\,\frac{\mathrm{d}}{\mathrm{d}t}\left(- \int_0^t\mathrm{d}t'\,\beta(t')\right) \nonumber \\[1em] &= \beta(t)\,\exp\left[- \int_0^t\mathrm{d}t'\,\beta(t')\right]. \end{align} \tag{16.9}\]
Say we observe arrivals of a Poisson process one after another. If we have \(n\) arrivals of a Poisson process with a variable rate occuring at times \(\mathbf{t} = \{t_1, t_2, \ldots, t_n\}\), where the times are ordered, the PDF is given by the product of the PDFs for each waiting time.
\[\begin{align} f(\mathbf{t} \mid \beta(t)) = \prod_{i=1}^n\beta(t_i)\exp\left[- \int_{t_{i-1}}^{t_i}\mathrm{d}t'\,\beta(t')\right] = \left(\prod_{i=1}^n\beta(t_i)\right)\exp\left[- \sum_{i=1}^n\int_{t_{i-1}}^{t_i}\mathrm{d}t'\,\beta(t')\right], \end{align} \tag{16.10}\]
where \(t_0\) is when we started observing (it is not a time of an arrival of a Poisson process), and we take \(t_0 = 0\). The integrals add together, giving
\[\begin{align} f(\mathbf{t} \mid \beta(t)) = \left(\prod_{i=1}^n\beta(t_i)\right)\exp\left[-\int_0^{t_n}\mathrm{d}t'\,\beta(t')\right]. \end{align} \tag{16.11}\]
In practice, we start observing at time \(t = 0\) and end at time \(t = T\). While \(t_n\) is the time of the last Poisson process we saw arrive, we should also take into account that we continued watching for time \(T - t_n\) and saw no arrivals during that time. We therefore have
\[\begin{align} f(\mathbf{t}, T;\beta(t)) &= f(\mathbf{t} \mid \beta(t)) \times P(\text{no arrivals between }t_n\text{ and }T) \nonumber \\[1em] &= \left(\prod_{i=1}^n\beta(t_i)\right)\exp\left[-\int_0^{t_n}\mathrm{d}t'\,\beta(t')\right]\,\exp\left[\int_{t_n}^T\mathrm{d}t'\,\beta(t')\right] \nonumber \\[1em] &= \left(\prod_{i=1}^n\beta(t_i)\right)\exp\left[-\int_0^{T}\mathrm{d}t'\,\beta(t')\right]. \end{align} \tag{16.12}\]
Note that in the special case of a constant rate, \(\beta(t) = \beta\), we get
\[\begin{align} f(\mathbf{t}, T;\beta(t)) &= \beta^n\,\mathrm{e}^{-\beta T}, \end{align} \tag{16.13}\]
In this case, the PDF is not dependent on the times, but only on the number of arrivals in the time interval. We could have derived by noting that the interrarrival times are Exponentially distributed.
\[\begin{align} f(\mathbf{t}, T;\beta(t)) &= \beta\,\mathrm{e}^{-\beta t_1}\,\left(\prod_{i=2}^{n-1}\beta\,\mathrm{e}^{-\beta (t_{i+1}-t_i)}\right)\mathrm{e}^{-\beta(T-t_n)} = \beta^n\,\mathrm{e}^{-\beta T}. \end{align} \tag{16.14}\]
Consider now a time interval from \(0\) to \(T\). We break this time interval into \(N\) small intervals, all of which are small enough such that, as before, \(\beta(t)\) does not change appreciably across a small interval and furthermore that the probability of more than one spike landing in a given small interval is zero. Let \(\Delta t_i = t_{i+1}-t_i\) be the width of small interval \(i\). For reasons that will become clear momentarily, we choose the widths of the small intervals such that the small interval width times the rate is the same for all intervals (which we will define as \(\xi\) for convenience), or
\[\begin{align} \beta(t_{i+1}) \,\Delta t_i = \xi \text{ for all }i. \end{align} \tag{16.15}\]
Since every \(\beta(t_{i+1}) \,\Delta t_i\) is the same, we can write
\[\begin{align} \xi = \frac{1}{N}\sum_{i=0}^N\,\beta(t_{i+1}) \,\Delta t_i = \frac{1}{N}\int_0^T\mathrm{d}t'\,\beta(t') \equiv \frac{\Lambda(T)}{N}, \end{align} \tag{16.16}\]
where we have considered small \(\Delta t_i\) and evaluated the Riemann sum, which is valid for unequal intervals provided the intervals are small enough. For convenience, we have defined
\[\begin{align} \Lambda(T) = \int_0^T\mathrm{d}t'\,\beta(t'). \end{align} \tag{16.17}\]
The probability of a spike not arriving in a small interval spanning from \(t_i\) to \(t_{i+1}\), as we have already worked out, is
\[\begin{align} P(\text{not arrived between } t_i \text{ and } t_{i+1}) \approx \mathrm{e}^{-\beta(t_{i+1})\,\Delta t_i} = \mathrm{e}^{-\xi}, \end{align} \tag{16.18}\]
such that the probability a spike does arrive in that interval is
\[\begin{align} P(\text{arrived between } t_i \text{ and } t_{i+1}) \approx 1 - \mathrm{e}^{-\beta(t_{i+1}) \Delta t_i} = 1-\mathrm{e}^{-\xi}. \end{align} \tag{16.19}\]
Again owing the memorylessness of spiking, the probability of getting a spike in a given small interval constitutes a Bernoulli trial with probability \(\theta\) of success. Then, the probability of getting \(n\) spikes in the interval \([0, T]\) is Binomially distributed,
\[\begin{align} f(n;N,\xi) \approx \frac{N!}{n!(N-n)!}\,(1 - \mathrm{e}^{-\xi})^n\,(\mathrm{e}^{-\xi})^{N-n}. \end{align} \tag{16.20}\]
We now consider the limit of large \(N\). In this limit,
\[\begin{align} \mathrm{e}^{-\xi} \approx 1 - \xi, \end{align} \tag{16.21}\]
such that
\[\begin{align} f(n;N,\xi) \approx \frac{N!}{n!(N-n)!}\,\xi^n\,(1-\xi)^{N-n}. \end{align} \tag{16.22}\]
Now, in the limit of large \(N\), \(N \gg n\) such that \(N - n \approx N\). Further, we have that \(N = \Lambda(T) / \xi\), such that
\[\begin{align} (1-\xi)^{N-n} \approx (1-\xi)^{N} = (1-\xi)^{\Lambda(T)/\xi}. \end{align} \tag{16.23}\]
Because
\[\begin{align} \lim_{\xi\to 0}(1-\xi)^{1/\xi} = \mathrm{e}^{-1}, \end{align} \tag{16.24}\]
we have
\[\begin{align} (1-\xi)^{N-n} \approx (1-\xi)^{\Lambda(T)/\xi} \approx \mathrm{e}^{-\Lambda(T)}. \end{align} \tag{16.25}\]
Finally, also have that in the limit of large \(N\),
\[\begin{align} N!/(N-n)! \approx N^n = \left(\frac{\Lambda(T)}{\xi}\right)^n. \end{align} \tag{16.26}\]
Putting everything together, we have
\[\begin{align} f(n;\beta(t), T) = \frac{(\Lambda(T))^n}{n!}\,\mathrm{e}^{-\Lambda(T)}, \end{align} \tag{16.27}\]
where
\[\begin{align} \Lambda(T) = \int_0^T\mathrm{d}t'\,\beta(t'), \end{align} \tag{16.28}\]
giving the PMF for the number of arrivals of a nonhomogenous Poisson process in an interval of length \(T\).
We can easily sample arrival times of a homogeneous Poisson process by noting that the interarrival times are Exponentially distributed. We keep drawing interarrival times from an Exponential distribution and cumulatively sum them to get the arrival times.
rng = np.random.default_rng()
def sample_homogeneous_poisson_process(beta, T):
"""Draw samples of arrival times of a homogeneous Poisson process."""
# Intialize arrival times
arrival_times = []
# Typical arrival time (Numpy likes this)
tau = 1 / beta
# Draw interarrival times and keep adding
t = rng.exponential(tau)
while t < T:
arrival_times.append(t)
t += rng.exponential(tau)
return np.array(arrival_times)
# Give it a whirl
sample_homogeneous_poisson_process(1, 20)array([ 0.30127478, 2.12937407, 3.36913533, 3.42966043, 3.75863601,
4.47142917, 4.5084683 , 4.55892601, 4.83689062, 5.61171643,
5.79705461, 5.94193015, 9.54332889, 9.80945717, 10.56075264,
12.96824828, 14.4488708 , 15.07726874, 16.50660407, 16.88488701,
18.12896431, 18.43671618])However, sampling arrival times of a nonhomogeneous Poisson is more challenging. The arrival time varies, so we cannot directly draw out of an Exponential distribution. To perform the sampling, we can use a thinning method. This works as follows. We set a fast rate \(\beta_u\) that is greater than \(\beta(t)\) everywhere on the interval \([0, T]\). We draw a time \(t_1\) out of an Exponential distribution with rate \(\beta_u\). We accept \(t_1\) as a sample out of the nonhomogeneous distribution with probability \(\beta(t_1)/\beta_u\). We then draw an interarrival time \(\delta t_2\) from an Exponential distribution parametrized with rate \(\beta_u\). We accept \(t_1 + \delta t_2\) as a draw from the nonhomogeneous distribution with probability \(\beta(t_1 + \delta t_2) / \beta_u\). We continue in this manner until the total time of the draws meets or exceeds \(T\). This is implemented in the function below.
def sample_nhpp(beta, beta_u, T, beta_args=()):
"""Draw arrival times of a nonhomogeneous Poisson process.
Parameters
----------
beta : function, call signature beta(t, *beta_args)
The function of the rate of arrivals as a function of time.
beta_u : float
A value of beta that is greater than beta(t) for all time.
T : float
The maximum time of observation.
beta_args : tuple, default ()
Any other arguments passed to the function `beta()`.
Returns
-------
output : Numpy array
Times for arrivals of the nonhomogeneous Poisson process.
Notes
-----
.. This is an implementation of the algorithm on page 85 of
Simulation, 5th Ed., by Sheldon Ross, Academic Press, 2013.
"""
samples = []
tau = 1 / beta_u
t = rng.exponential(tau)
while t < T:
r = beta(t, *beta_args) * tau
if r > 1:
raise RuntimeError('beta_u is less than beta; sampling invalid.')
if np.random.uniform() <= r:
samples.append(t)
t += rng.exponential(tau)
return np.array(samples)Let’s take this for a spin. We will have \(\beta(t)\) given by a sinusoid.
beta = lambda t: 1.15 + np.sin(t/10)
arrival_times = sample_nhpp(beta, 2.15, 500)
# Plot the arrivals
p = iqplot.strip(arrival_times, marker='dash', marker_kwargs=dict(alpha=0.3))
bokeh.io.show(p)Nice!
If we want to speed things up, we can break the interval up into smaller segments and for each interval we make sure that we choose a \(\beta_u\) specific to that interval that is close to the maximum \(\beta(t)\) on that interval. The code (which also includes other speed enhancements) to achieve this is below.
def sample_nhpp(beta, beta_j, t_j, beta_args=()):
"""Draw arrival times of a nonhomogeneous Poisson process.
Parameters
----------
beta : function, call signature beta(t, *beta_args)
The function of the rate of arrivals as a function of time.
beta_j : scalar or array_like
If scalar, a value of beta that is greater than beta(t)
for all time. If array_like, then beta_j[j] > beta(t) for
all times between t_j[j-1] and t_j[j].
t_j : scalar or array_like
If scalar, the maximum time of observation. If array_like, must
be the same length of `beta_j`. beta_j[j] is the value
of the the upper bound of the rate for the interval between
t[j-1] and t[j].
beta_args : tuple, default ()
Any other arguments passed to beta.
Returns
-------
output : Numpy array
Times for arrivals of the nonhomogeneous Poisson process.
Notes
-----
.. This is an implementation of the algorithm on page 86 of
Simulation, 5th Ed., by Sheldon Ross, Academic Press, 2013.
"""
# Convert scalar inputs to arrays
if np.isscalar(beta_j):
beta_j = np.array([beta_j])
if np.isscalar(t_j):
t_j = np.array([t_j])
# Make sure dimensions match
if len(beta_j) != len(t_j):
raise RuntimeError(
f'`beta_j` is length {len(beta_j)} '
f'and `t_j` is length {len(t_j)}. '
'They must have the same length.'
)
return _sample_nhpp(beta, beta_j, t_j, beta_args)
def _sample_nhpp(beta, beta_j, t_j, beta_args=()):
"""
Parameters
----------
beta : function, call signature beta(t, *beta_args)
The function of the rate of arrivals as a function of time.
beta_j : Numpy array
If scalar, a value of beta that is greater than beta(t)
for all time. If array_like, then beta_j[j] > beta(t) for
all times between t_j[j-1] and t_j[j].
t_j : Numpy array
Must be the same length of `beta_j`. beta_j[j] is the value
of the the upper bound of the rate for the interval between
t[j-1] and t[j].
beta_args : tuple, default ()
Any other arguments passed to beta.
Returns
-------
output : Numpy array
Times for arrivals of the nonhomogeneous Poisson process.
Notes
-----
.. This is an implementation of the algorithm on page 86 of
Simulation, 5th Ed., by Sheldon Ross, Academic Press, 2013.
"""
# Number of samples to take before concatenating arrays
n_samples = 1000
# Initializations
t = 0.0 # time
j = 0 # Index in beta_j and t_j arrays
i = 0 # index in sample array
n = 0 # total number of samples drawn
x_from_boundary = False # If we've hit a boundary of
# Array for storing subtrajectory
samples = np.empty(n_samples, dtype=float)
# Output array for all samples
samples_output = np.array([], dtype=float)
# Loop until done (we exceed final time point)
not_done = True
while not_done:
# Take samples until we fill array
# We do it this way for speed to avoid list append operations
while not_done and i < n_samples:
if x_from_boundary:
x = (x - t_j[j] + t) * beta_j[j] / beta_j[j+1]
t = t_j[j]
j += 1
else:
x = np.random.exponential(1 / beta_j[j])
if t + x > t_j[j]:
# If we got here, we went past the edge of this interval
if j == len(t_j) - 1:
not_done = False
else:
x_from_boundary = True
else:
t += x
x_from_bounday = False
if np.random.uniform() <= beta(t, *beta_args) / beta_j[j]:
samples[i] = t
i += 1
n += 1
samples_output = np.concatenate((samples_output, samples))
i = 0
return np.array(samples_output[:n])%load_ext watermark
%watermark -v -p numpy,iqplot,bokeh,jupyterlabPython implementation: CPython
Python version : 3.12.9
IPython version : 9.1.0
numpy : 2.1.3
iqplot : 0.3.7
bokeh : 3.6.2
jupyterlab: 4.4.2