Poisson Process

Proceso Poisson

Definición:

Sean \(T_1\), \(T_2\), …, \(T_n\) v.a.i.i.d con distribución exponencial de parámetro \(\lambda\), : \(\mathcal{E}xp(\lambda)\). : \(\mathcal{W}_0=0\), \(\mathcal{W}_n=T_1, T_2, ..., T_n\) para \(n\geq 1\). Definimos el proceso Poisson de parámetro o intensidad \(\lambda\) por

\[ \begin{equation} N(t) = máx \{n\geq 1, \ \ \mathcal{W}_n = T_1 + T_2 + ...+ T_n \leq t\} \end{equation} \] Las variables aleatorias \(T_n\) representan los intervalos de tiempo entre eventos sucesivos, y \(\mathcal{W}_n = T_1, T_2, ..., T_n\) representa el instante en el que ocurre el n-ésimo evento, y \(N(t)\) es el número de eventos que han ocurrido hasta el instante t.

Una definción más común, sea \(\{N(t), t \geq 0\}\) un proceso estocástico que cuenta el número de eventos hasta el tiempo \(t\). Se dice que \(\{N(t)\}\) es un proceso de Poisson de parámetro o intensidad \(\lambda > 0\) si cumple las siguientes propiedades:

\(N(0) = 0\)
Incrementos independientes: Para cualquier colección de intervalos disjuntos \([t_0, t_1], [t_2, t_3], \dots\), los incrementos \[ N(t_1)-N(t_0), \ N(t_3)-N(t_2), \dots \] son variables aleatorias independientes.
Incrementos estacionarios: La distribución de los incrementos depende solo de la longitud del intervalo: \[ N(t+s)-N(s) \sim \text{Poisson}(\lambda t), \quad \forall s,t \ge 0. \]
No ocurrencia simultánea de eventos: \[ \lim_{\Delta t \to 0} P(\text{más de un evento en } [t, t+\Delta t]) = 0. \]

Proposición:

La variable aleatoria \(N(t)\) tiene distribución Poisson con parámetro \((\lambda t)\), es decir, para cualquier \(t>0\), y para \(n=0, 1, ...\)

\[\begin{equation} P(N(t) = n) = e^{-\lambda t} \frac{(\lambda t)^n}{n!} \end{equation}\]

Su valor esperado, y varianza son

\[\mathbb{E}[N(t)] = \lambda t \]

\[\mathbb{V}ar(N(t)) = \lambda t \] Para un proceso Poisson con parámetro lambda \(\lambda\), y \(s<t\) la covarianza entre \(N(s)\) y \(N(t)\) es \[ \begin{equation} \mathbb{C}ov(N(s), N(t)) = \lambda s \end{equation} \] Distribución de los Tiempos de Llegada

El \(n\)-ésimo tiempo de llegada \(W_n\) tiene una distribución Gamma:

\[ W_n \sim \text{Gamma}(n, \lambda) \] con función de densidad de probabilidad: \[ f_{W_n}(t) = \frac{\lambda^n t^{\,n-1} e^{-\lambda t}}{(n-1)!}, \quad t > 0 \]

Aquí, \(W_n = T_1 + T_2 + \dots + T_n\) representa el instante en el que ocurre el \(n\)-ésimo evento, siendo \(T_1, T_2, \dots, T_n\) variables aleatorias independientes e idénticamente distribuidas con distribución exponencial de parámetro \(\lambda\), \(T_i \sim \mathcal{E}xp(\lambda)\).

Simulación

Este código define una clase \(StochasticSimulation()\) para simular procesos estocásticos en tiempo discreto. El método \(__init__\) inicializa el número de pasos, repeticiones (\(n\_times\)), el tiempo total \(T\) y genera un vector de tiempos \(t\) con incremento \(dt\). El método \(PoissonProcess(rate, max\_time)\) simula trayectorias. Se realizan \(n\_times\) simulaciones en paralelo, almacenando cada trayectoria en filas del arreglo. Finalmente, convierte el resultado en un DataFrame de pandas.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.cm as cm
import matplotlib.patches as patches
sns.set_theme()

class StochasticSimulation():

    def __init__(self, T=None, steps=100, n_times=1):
        self.steps = steps
        self.n_times = n_times
        # Si T no se proporciona, usar T = steps
        if T is None:
            self.T = steps
        else:
            self.T = T

        self.dt = self.T/ steps
        self.t = np.linspace(0, self.T, steps+1)

    # Procesos Poisson
    def PoissonProcess(self, rate, max_time=1000):
          """Event-based simulation"""
          time_array = np.linspace(0, max_time, self.steps + 1)
          trajectories = np.zeros((self.n_times, self.steps + 1))
          
          for i in range(self.n_times):
              events = []
              current_time = 0
              while current_time <= max_time:
                  inter_arrival = np.random.exponential(1/rate)
                  current_time += inter_arrival
                  if current_time <= max_time:
                      events.append(current_time) 
              # Contamos eventos
              event_count = 0
              for j, t in enumerate(time_array):
                  while event_count < len(events) and events[event_count] <= t:
                      event_count += 1
                  trajectories[i, j] = event_count
          # retornamos un dataframe
          return pd.DataFrame(trajectories)

En este ejemplo, primero se definen los parámetros de la simulación: se tomarán \(𝑁=10,000\) pasos, lo que permitirá generar una trayectoria muy detallada del proceso, y se simulará \(𝑀=1\), una única trayectoria.

A continuación, se crea una instancia de la clase \(StochasticSimulation()\), pasando como argumentos el número de pasos \(N\)y el número de trayectorias \(M\) a simular. Esta clase se encarga de manejar los cálculos necesarios para generar procesos estocásticos.

Finalmente, se llama al método \(PoissonProcess(rate, max\_time)\) de la instancia, recibe como argumento \(rate=0.01\) parámetro lambda, así como \(max\_time=N\), con los cuales nos genera la trayectoria simulada. El resultado se almacena en df.

Simulación de un proceso Poisson en una dimensión

Simulación de una trajectoria

# Parámetros
N = 10000        # número de pasos
M = 1            # número de trayectorias simuladas
rate= 0.01
max_time = N
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.PoissonProcess(rate, max_time)

# Colormap
cmap = cm.get_cmap("Grays", M)
# lista de t
t_max = [*range(0, N, 1)]
T = sim.T
# sdt teórica : lambda*t
std_teorica = np.sqrt(np.array(t_max)*rate)
# Media teórica
mean_B = [t*rate for t in t_max]
# Intervalos de confianza 95% t(1−t/T)
conf_upper =  mean_B + 1.96 * std_teorica
conf_lower =  mean_B - 1.96 * std_teorica
# Gráfico
plt.figure(figsize=(14,6))
for i in range(M):
    plt.step( df.columns.to_list(), df.iloc[i,:].to_list(), color=cmap(i), linewidth=0.5)
plt.plot(t_max, conf_upper, 'blue', linewidth=0.3, linestyle='--')
plt.plot(t_max, conf_lower, 'blue', linewidth=0.3, linestyle='--')
plt.fill_between(t_max, conf_lower, conf_upper, color='lightblue', alpha=0.2, label='Intervalo de confianza 95% (teórico)')
plt.title(f"Proceso Poisson[{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("Events")
plt.grid(True)
plt.legend()
plt.show()

Simulación de múltiples trajectorias

# Parámetros
N = 10000        # número de pasos
M = 1000            # número de trayectorias simuladas
rate= 0.01
max_time = N
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.PoissonProcess(rate, max_time)

# Colormap
cmap = cm.get_cmap("Grays", M*2)
# lista de t
t_max = [*range(0, N, 1)]
T = sim.T
# sdt teórica : lambda*t
std_teorica = np.sqrt(np.array(t_max)*rate)
# Media teórica
mean_B = [t*rate for t in t_max]
# Intervalos de confianza 95% t(1−t/T)
conf_upper =  mean_B + 1.96 * std_teorica
conf_lower =  mean_B - 1.96 * std_teorica
# Gráfico
plt.figure(figsize=(14,6))
for i in range(M):
    plt.step( df.columns.to_list(), df.iloc[i,:].to_list(), color=cmap(i), linewidth=0.5)
plt.plot(t_max, conf_upper, 'blue', linewidth=0.3, linestyle='--')
plt.plot(t_max, conf_lower, 'blue', linewidth=0.3, linestyle='--')
plt.fill_between(t_max, conf_lower, conf_upper, color='lightblue', alpha=0.2, label='Intervalo de confianza 95% (teórico)')
plt.title(f"Proceso Poisson[{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("Events")
plt.grid(True)
plt.legend()
plt.show()

Veamos las distribuciones a distintos tiempos

# Tiempos a graficar
steps = [10,100,300,500]
# Generamos gráfico
fig, axes = plt.subplots(2, 2, figsize=(10, 6))
axes = axes.flatten() 
for i, step in enumerate(steps):
    data = df[step]
    mean = data.mean()
    var = data.var()
    sns.histplot(data, bins=30, kde=True, ax=axes[i])
    axes[i].set_title(f'Distribución en paso {step}')
    axes[i].set_xlabel('Posición')
    axes[i].set_ylabel('Frecuencia')
    axes[i].legend([f'Media = {mean:.2f}, Varianza = {var:.2f}'])
plt.tight_layout()
plt.show()

POdemos observar que para cada valor de \(t\) los parámetros, media, y varianza, se aproximan con los valores teóricos.

Simulación de un Proceso Poisson en dos dimensiones

Finalmente, tomamos los datos de la simulación anterior y gráficamos pares de caminatas aleatorias para ver su evolución en dos dimensiones.

# Colormap (gradiente de gris)
cmap = cm.get_cmap("Greys", M*2)
centro = (50, 50)
circulo = patches.Circle(centro, 100, edgecolor='grey', facecolor='black', alpha=0.8, lw=1, ls='--')
# Crear la figura y los ejes
fig, ax = plt.subplots(figsize=(14, 10))
ax.add_patch(circulo)
# Proceso Poisson
for i in range(df.shape[0] - 1):
    if i < 500: 
        ax.plot(df.iloc[i][:], df.iloc[i + 1][:], linewidth=0.2, color=cmap(i))  # Normaliza el valor de i
    else:
        continue
plt.title('Simulación de un Proceso Poisson: {} steps and {} paths'.format(N, M))
plt.xlabel('N(t)')
plt.ylabel('N(t)')
plt.gca().set_xticklabels([])  
plt.gca().set_yticklabels([])  
plt.grid(True)
plt.show()

Referencias

Ubbo F Wiersema (2008). Brownian Motion Calculus. John Wiley & Sons Ltd.
Sheldon M. Ross(1995). Stochastic Processes. Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons.
LefebvreMario(2007). Applied Stochastics Processes. Springer, New York.