How to simulate Poisson and Compound Poisson process

Question

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?

Answer 1

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