#!/usr/bin/env python
# coding: utf-8

# > This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.
# 

# # 13.2. Simulating a Poisson process

# 1. Let's import NumPy and matplotlib.

# In[ ]:


import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')


# 2. Let's specify the `rate`: the average number of events per second.

# In[ ]:


rate = 20.  # average number of events per second


# 3. First, we will simulate the process using small time bins of 1 millisecond.

# In[ ]:


dt = .001  # time step
n = int(1./dt)  # number of time steps


# 4. On every time bin, the probability that an event occurs is about $\textrm{rate} \times dt$, if $dt$ is small enough. Besides, since the Poisson process has no memory, the occurrence of an event is independent from one bin to another. Therefore, we can sample Bernoulli random variables in a vectorized way in order to simulate our process.

# In[ ]:


x = np.zeros(n)
x[np.random.rand(n) <= rate*dt] = 1


# The vector `x` contains zeros and ones on all time bins, *one* corresponding to the occurrence of an event.

# In[ ]:


x[:10]


# 5. Let's display the simulated process. We draw a vertical line on each non-zero time bin.

# In[ ]:


plt.figure(figsize=(6,2));
plt.vlines(np.nonzero(x)[0], 0, 1);
plt.xticks([]); plt.yticks([]);


# 6. Another way of representing that same object consists in considering the associated **counting process** $N(t)$: the number of events that have occurred until time $t$. Here, we can display this process using the function `cumsum`.

# In[ ]:


plt.figure(figsize=(6,4));
plt.plot(np.linspace(0., 1., n), np.cumsum(x));
plt.xlabel("Time");
plt.ylabel("Counting process");


# 7. The other (and more efficient) way of simulating the homogeneous Poisson process is to use the property that the time interval between two successive events is an exponential random variable. Furthermore, these intervals are independent, so that we can sample these intervals in a vectorized way. Finally, we get our process by summing cumulatively all those intervals.

# In[ ]:


y = np.cumsum(np.random.exponential(1./rate, size=int(rate)))


# The vector `y` contains another realization of our Poisson process, but the data structure is different. Every component of the vector is the time of an event.

# In[ ]:


y[:10]


# 8. Finally, let's display the simulated process.

# In[ ]:


plt.figure(figsize=(6,2));
plt.vlines(y, 0, 1);
plt.xticks([]); plt.yticks([]);


# > You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).
# 
# > [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages).