-2
$\begingroup$

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?

enter image description here

$\endgroup$
2
$\begingroup$

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

sample poisson process

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

Poisson variables vs. count from Poisson processes

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

enter image description here

$\endgroup$
  • 2
    $\begingroup$ 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$ $\endgroup$ – StackG Aug 10 '20 at 1:15
  • $\begingroup$ Thank you! But how to make the poisson trajectories be piecewise constant instead of given as points? $\endgroup$ – Math122 Aug 10 '20 at 14:12
  • 1
    $\begingroup$ 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 $\endgroup$ – StackG Aug 10 '20 at 14:20
  • 1
    $\begingroup$ I add drawstyle='steps-post' and works :) Thanks for help $\endgroup$ – Math122 Aug 10 '20 at 14:32
  • 1
    $\begingroup$ Cant you just add a point at (0,0)? $\endgroup$ – StackG Aug 10 '20 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.