# How to simulate Poisson and Compound Poisson process

someone knows, maybe websites / blogs where I can find tips (preferably ready codes) to simulate the trajectory of processes? So far I only need the Poisson process and the compound Poisson process but I would like to learn to simulate more advanced processes as well so I will be grateful for every source

EDIT: I edit code from anserws a little bit:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import poisson

mu = 3
num_events_simulated = 5

def generate_poisson_process(mu, num_events):
time_intervals = -np.log(np.random.random(num_events)) / mu
total_events = time_intervals.cumsum()
events = pd.DataFrame(np.ones(num_events), index=total_events)
events[0] = events[0].cumsum()

return events

plt.plot(generate_poisson_process(mu, num_events_simulated), marker='o', drawstyle='steps-post')
plt.title("Sample Poisson Process")
plt.xlabel("time")
plt.ylabel("events")


but the trajectory starts from the first jump. How to add to the graph this fragment in which the trajectory starts from zero and is 0 until 1 jump?

The time between two events in a poisson distribution has an exponential distribution, so the easiest thing to do is simulate a sequence of exponentially distributed variables and use these as the times between events, as discussed in this primer.

To simulate variables given a uniform RNG, we need the reverse CDF of the distribution, which maps uniform distributions to our distribution of choice

For the exponential distribution this is just \begin{align} F^{-1}(x) = {\frac {-\log(1-x)} {\lambda}} \end{align} where $$\lambda$$ is the poisson parameter

We can simulate that using numpy as follows

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import poisson

mu = 15
num_events_simulated = 50

def generate_poisson_process(mu, num_events):
time_intervals = -np.log(np.random.random(num_events)) / mu
total_events = time_intervals.cumsum()
events = pd.DataFrame(np.ones(num_events), index=total_events)
events[0] = events[0].cumsum()

return events

plt.plot(generate_poisson_process(mu, num_events_simulated), marker='o', linestyle='none')
plt.title("Sample Poisson Process")
plt.xlabel("time")
plt.ylabel("events")
plt.legend()


which generates

So far, so good, but how can we demonstrate that this is a poisson process? Well, for a poisson process, the number of events in a period of time $$\tau$$ is distributed as $$Poi(\tau\lambda)$$, and numpy can generate these directly for us to compare.

Here, we generate 100,000 processes and compare the number of events in $$\tau = 1$$ to the results of the variables generated directly:

results = []
for x in range(100000):
process = generate_poisson_process(mu, num_events_simulated)
results.append(process[:1][0].iloc[-1])

plt.hist(results, bins=np.linspace(0, 35, 36), alpha=0.5, label='simulated poisson', ec='black')

r = poisson.rvs(mu, size=100000)
plt.hist(r, bins=np.linspace(0, 35, 36), alpha=0.5, label='counting process', ec='black')

plt.title("Poisson-Distributed Variables")
plt.ylabel("Count")
plt.xlabel("X")
plt.legend()


which generates

Match is very good, so we're happy with the path generation process!

The compound poisson process is a simple extension of this as long as you know which secondary distribution you want to use, and can generate variables distributed according to it. Here is an example of a compound poisson process generated using binomial $$X ~ B(10,0.5)$$ as the seocndary distribution:

def binomial_generator(num_events):
return np.random.binomial(10, 0.5, num_events)

def generate_compound_poisson_process(mu, num_events, generator):
time_intervals = -np.log(np.random.random(num_events)) / mu
total_events = time_intervals.cumsum()
events = pd.DataFrame(generator(num_events), index=total_events)
events[0] = events[0].cumsum()

return events

plt.plot(generate_compound_poisson_process(mu, num_events_simulated, binomial_generator), marker='o', linestyle='none')
plt.title("Sample Compound Poisson Process (Binomial)")
plt.xlabel("time")
plt.ylabel("events")
plt.legend()


which generates

• Another approach is discussed here: columbia.edu/~ks20/4703-Sigman/4703-07-Notes-PP-NSPP.pdf where instead of creating exponential variables, you just create a single poisson-distributed variable for time $\tau$, call this $N$, and then generate $N$ uniform variables $u_0$ to $u_N$, which you line up from smallest to largest and then the event times are just $\tau * u_i$ for each $i$ – StackG Aug 10 '20 at 1:15
• Thank you! But how to make the poisson trajectories be piecewise constant instead of given as points? – Math122 Aug 10 '20 at 14:12
• I've plotted them here as points, but if you want to know the value at any given time then just slice to that time in the dataframe coming out of generate_poisson_process, and the most recent value in the column 0 is the value that you need – StackG Aug 10 '20 at 14:20
• I add drawstyle='steps-post' and works :) Thanks for help – Math122 Aug 10 '20 at 14:32
• Cant you just add a point at (0,0)? – StackG Aug 10 '20 at 22:22