poisson point process.py

import numpy as np
import matplotlib
matplotlib.use("TkAgg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from scipy.stats import poisson, expon
import time

# Set random seed for reproducibility
np.random.seed(42)

# Global variables for figure and axes
fig = None
axes = None


def init_figure():
    global fig, axes
    # Set up the figure with subplots
    fig, axes = plt.subplots(3, 2, figsize=(15, 14))
    # Flatten the axes array for easier indexing
    axes = axes.flatten()


# Function to plot example of a Poisson distribution
def plot_poisson_distribution():
    global axes
    ax = axes[0]

    lambda_val = 5
    x = np.arange(0, 15)
    pmf = poisson.pmf(x, lambda_val)

    ax.bar(x, pmf, alpha=0.7)
    ax.axvline(lambda_val, color='r', linestyle='--', label=f'λ = {lambda_val}')
    ax.set_title(f'Poisson Distribution (λ = {lambda_val})')
    ax.set_xlabel('Number of Events (k)')
    ax.set_ylabel('Probability P(X = k)')
    ax.legend()
    ax.grid(alpha=0.3)


# Function to simulate and plot homogeneous Poisson process on a line (1D)
def simulate_homogeneous_1d():
    global axes
    ax = axes[1]

    # Parameters
    lambda_val = 1.5  # Rate parameter
    T = 10  # Time interval [0, T]

    # Step 1: Generate number of points
    N = np.random.poisson(lambda_val * T)

    # Step 2: Generate random positions (uniform)
    points = np.random.uniform(0, T, N)
    points.sort()  # Sort to get the event times in order

    # Plot the points as a counting process
    counts = np.arange(1, N + 1)

    # Plot the points as impulses
    for p in points:
        ax.axvline(x=p, ymin=0, ymax=0.05, color='b', alpha=0.7)

    # Plot the counting process N(t)
    ax.step(np.concatenate(([0], points)), np.concatenate(([0], counts)), 'r-', where='post', label='N(t)')

    ax.set_title(f'Homogeneous Poisson Process (λ = {lambda_val})')
    ax.set_xlabel('Time (t)')
    ax.set_ylabel('Count N(t)')
    ax.set_xlim([0, T])
    ax.set_ylim([0, N + 1])
    ax.grid(alpha=0.3)
    ax.legend()

    # Return the points for other demonstrations
    return points


# Function to demonstrate exponential interarrival times
def plot_interarrival_times(points):
    global axes
    ax = axes[2]

    # Calculate interarrival times
    interarrival_times = np.diff(np.concatenate(([0], points)))

    # Plot histogram of interarrival times
    ax.hist(interarrival_times, bins=15, density=True, alpha=0.7)

    # Plot theoretical exponential PDF
    lambda_val = 1 / np.mean(interarrival_times)
    x = np.linspace(0, max(interarrival_times), 1000)
    pdf = expon.pdf(x, scale=1 / lambda_val)
    ax.plot(x, pdf, 'r-', linewidth=2, label=f'Exponential PDF (λ = {lambda_val:.2f})')

    ax.set_title('Interarrival Times Distribution')
    ax.set_xlabel('Interarrival Time')
    ax.set_ylabel('Density')
    ax.grid(alpha=0.3)
    ax.legend()


# Function to simulate and plot homogeneous Poisson process in 2D
def simulate_homogeneous_2d():
    global axes
    ax = axes[3]

    # Parameters
    lambda_val = 50  # Intensity parameter (points per unit area)
    xmin, xmax = 0, 10
    ymin, ymax = 0, 10
    area = (xmax - xmin) * (ymax - ymin)

    # Step 1: Generate number of points
    N = np.random.poisson(lambda_val * area)

    # Step 2: Generate random positions (uniform in both coordinates)
    x_points = np.random.uniform(xmin, xmax, N)
    y_points = np.random.uniform(ymin, ymax, N)

    # Plot the points
    ax.scatter(x_points, y_points, s=10, alpha=0.7)
    ax.set_title(f'Homogeneous Spatial Poisson Process (λ = {lambda_val})')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_xlim([xmin, xmax])
    ax.set_ylim([ymin, ymax])
    ax.grid(alpha=0.3)

    # Create a proper ScalarMappable with a plot
    sm = plt.cm.ScalarMappable(plt.Normalize(lambda_val, lambda_val), plt.cm.Blues)
    fig.colorbar(sm, ax=ax, label='Intensity λ (points per unit area)')


# Function to simulate and plot inhomogeneous Poisson process in 1D
def simulate_inhomogeneous_1d():
    global axes
    ax = axes[4]

    # Parameters
    T = 10

    # Intensity function λ(t) = 2 + sin(t)
    def lambda_t(t):
        return 2 + np.sin(t)

    # Maximum value of λ(t) for rejection sampling
    lambda_max = 3

    # Step 1: Generate points using rejection sampling
    points = []
    t = 0
    while t < T:
        # Generate next exponential waiting time with rate lambda_max
        dt = np.random.exponential(1 / lambda_max)
        t += dt

        if t < T:
            # Generate a uniform random number for acceptance/rejection
            u = np.random.uniform(0, 1)

            # Accept the point with probability λ(t)/λ_max
            if u <= lambda_t(t) / lambda_max:
                points.append(t)

    points = np.array(points)

    # Plot the intensity function
    t_values = np.linspace(0, T, 1000)
    lambda_values = lambda_t(t_values)
    ax.plot(t_values, lambda_values, 'g-', label='λ(t) = 2 + sin(t)')

    # Plot the events
    for p in points:
        ax.axvline(x=p, ymin=0, ymax=0.05, color='r', alpha=0.7)

    # Plot the counting process
    counts = np.arange(1, len(points) + 1)
    ax.step(np.concatenate(([0], points)), np.concatenate(([0], counts)),
            'b-', where='post', label='N(t)', alpha=0.7)

    ax.set_title('Inhomogeneous Poisson Process')
    ax.set_xlabel('Time (t)')
    ax.set_ylabel('Intensity λ(t) / Count N(t)')
    ax.set_xlim([0, T])
    ax.legend()
    ax.grid(alpha=0.3)


# Function to simulate and plot inhomogeneous Poisson process in 2D
def simulate_inhomogeneous_2d():
    global axes
    ax = axes[5]

    # Parameters
    xmin, xmax = -5, 5
    ymin, ymax = -5, 5

    # Intensity function λ(x,y) = exp(-(x^2 + y^2))
    def lambda_xy(x, y):
        return np.exp(-(x ** 2 + y ** 2))

    # Maximum value of λ(x,y) for rejection sampling
    lambda_max = 1

    # Generate a large number of candidate points
    N_candidates = 1000
    x_candidates = np.random.uniform(xmin, xmax, N_candidates)
    y_candidates = np.random.uniform(ymin, ymax, N_candidates)

    # Apply thinning through rejection sampling
    accepted_indices = []
    for i in range(N_candidates):
        x, y = x_candidates[i], y_candidates[i]
        intensity = lambda_xy(x, y)

        # Accept with probability λ(x,y)/λ_max
        if np.random.uniform(0, 1) <= intensity / lambda_max:
            accepted_indices.append(i)

    # Get the accepted points
    x_points = x_candidates[accepted_indices]
    y_points = y_candidates[accepted_indices]

    # Create a grid for displaying the intensity function
    x_grid = np.linspace(xmin, xmax, 100)
    y_grid = np.linspace(ymin, ymax, 100)
    X, Y = np.meshgrid(x_grid, y_grid)
    Z = lambda_xy(X, Y)

    # Plot the intensity function as a heatmap
    contour = ax.contourf(X, Y, Z, levels=50, cmap='viridis', alpha=0.5)
    fig.colorbar(contour, ax=ax, label='Intensity λ(x,y)')

    # Plot the points
    ax.scatter(x_points, y_points, color='r', s=10, alpha=0.7)

    ax.set_title('Inhomogeneous Spatial Poisson Process\nλ(x,y) = exp(-(x² + y²))')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_xlim([xmin, xmax])
    ax.set_ylim([ymin, ymax])
    ax.grid(alpha=0.3)


def main():
    # Initialize the figure and axes
    init_figure()

    # Plot examples of each concept
    plot_poisson_distribution()
    points = simulate_homogeneous_1d()
    plot_interarrival_times(points)
    simulate_homogeneous_2d()
    simulate_inhomogeneous_1d()
    simulate_inhomogeneous_2d()

    plt.tight_layout()
    plt.savefig('poisson_processes_demo.png', dpi=300)
    plt.show()


if __name__ == "__main__":
    main()