Stochastic Processes for Computer Science: Poisson, Brownian Motion, Queueing and Reliability

The Poisson process is the universal model of random arrivals. Whether you are modeling packet arrivals at a router, job submissions to a scheduler, failures in a distributed storage system, or mutations in a DNA sequence, the Poisson process is almost certainly the right first approximation. Its ubiquity derives from the Poisson limit theorem: the superposition of many independent, sparse point processes converges to a Poisson process. This is the Palm-Khintchine theorem, and it explains why the Poisson process appears everywhere from telephony (Erlang, 1909) to cloud computing (capacity planning, 2024).

But the Poisson process is just the beginning. Continuous-time Markov chains (CTMCs) extend the Markov property to continuous time, providing the mathematical framework for queueing theory (birth-death processes, M/M/1 queues, Jackson networks). Brownian motion—the scaling limit of random walks—provides the heavy-traffic approximations that let us reason about congested systems. And the reflection principle, a beautiful symmetry argument, gives exact hitting-time distributions for Brownian motion, which translate into bounds on buffer overflow probabilities, response time tail latencies, and time-to-failure in reliability models.

1. The Poisson Process

Definition 1.1. A Poisson process with rate (\lambda > 0) is a counting process ((N(t))_{t \geq 0}) satisfying:

1. (N(0) = 0).
2. Independent increments: For (0 \leq t_1 < t_2 \leq t_3 < t_4), (N(t_2) - N(t_1)) is independent of (N(t_4) - N(t_3)).
3. Stationary increments: (N(t + s) - N(t) \sim \mathrm{Poisson}(\lambda s)) for all (t, s \geq 0).

Equivalently, the inter-arrival times (T_1, T_2, \ldots) between consecutive events are independent (\mathrm{Exponential}(\lambda)) random variables.

Theorem 1.1 (Memoryless Property). The exponential distribution is the only continuous distribution with the memoryless property: for (s, t \geq 0),

[ \mathbb{P}(T > s + t \mid T > s) = \mathbb{P}(T > t) ]

Proof. (\mathbb{P}(T > s + t \mid T > s) = e^{-\lambda(s+t)} / e^{-\lambda s} = e^{-\lambda t} = \mathbb{P}(T > t)). Uniqueness follows from solving the functional equation (\bar{F}(s+t) = \bar{F}(s)\bar{F}(t)), whose only continuous solution is (\bar{F}(t) = e^{-\lambda t}). ∎

1.1 Properties of the Poisson Process

Superposition: The sum of independent Poisson processes with rates (\lambda_1, \ldots, \lambda_k) is a Poisson process with rate (\sum \lambda_i). This is why aggregate arrivals at a data center (with thousands of independent users) look Poisson.

Thinning: Given a Poisson process with rate (\lambda), independently classify each event as type (i) with probability (p_i). The resulting type-(i) process is Poisson with rate (\lambda p_i), and the processes for different types are independent. This models, for instance, a router that classifies packets into different QoS classes.

PASTA property: “Poisson Arrivals See Time Averages.” For a queueing system in equilibrium, arrivals from a Poisson process see the system in its time-averaged state distribution—they are not biased toward busy or idle periods. This is a consequence of the independent-increments property and is essential for analyzing M/G/1 queues.

2. Continuous-Time Markov Chains

A continuous-time Markov chain (CTMC) ((X_t){t \geq 0}) on a countable state space (\Omega) is defined by transition rates (q(x, y)) for (x \neq y) (the rate of transitioning from (x) to (y)). The total rate out of (x) is (q(x) = \sum{y \neq x} q(x, y)), and the holding time in state (x) is (\mathrm{Exponential}(q(x))). Upon leaving (x), the chain jumps to (y) with probability (q(x, y) / q(x)).

The generator matrix (or Q-matrix) (Q) has entries (Q(x, y) = q(x, y)) for (x \neq y) and (Q(x, x) = -q(x)). The transition probabilities (P_t(x, y) = \mathbb{P}(X_t = y \mid X_0 = x)) satisfy the Kolmogorov forward equations:

[ \frac{d}{dt} P_t = P_t Q ]

and the backward equations (\frac{d}{dt} P_t = Q P_t). The solution is the matrix exponential (P_t = e^{tQ} = \sum_{n=0}^\infty (tQ)^n / n!).

Definition 2.1. A distribution (\pi) is stationary for the CTMC if (\pi Q = 0), i.e., for all (y),

[ \sum_{x \neq y} \pi(x) q(x, y) = \pi(y) q(y) ]

This is the balance equation: the total rate of transitions into (y) equals the total rate out of (y).

2.1 Birth-Death Processes

A birth-death process is a CTMC on (\mathbb{N}) where from state (n), the only possible transitions are to (n+1) (birth, rate (\lambda_n)) or to (n-1) (death, rate (\mu_n)). The generator is tridiagonal. The stationary distribution (if it exists) satisfies the detailed balance equations:

[ \pi(n) \lambda_n = \pi(n+1) \mu_{n+1} ]

Solving recursively:

[ \pi(n) = \pi(0) \prod_{k=0}^{n-1} \frac{\lambda_k}{\mu_{k+1}}, \quad \pi(0) = \left(1 + \sum_{n=1}^\infty \prod_{k=0}^{n-1} \frac{\lambda_k}{\mu_{k+1}}\right)^{-1} ]

The M/M/1 queue is a birth-death process with (\lambda_n = \lambda) (constant arrival rate) and (\mu_n = \mu) (constant service rate). The stationary distribution is geometric: (\pi(n) = (1 - \rho) \rho^n), where (\rho = \lambda / \mu < 1) for stability.

3. Brownian Motion and the Reflection Principle

Brownian motion (Wiener process) is the fundamental continuous-time, continuous-state stochastic process. It arises as the scaling limit of random walks: if (S_n = \sum_{i=1}^n X_i) with (\mathbb{E}[X_i] = 0), (\mathrm{Var}(X_i) = 1), then (S_{\lfloor nt \rfloor} / \sqrt{n} \Rightarrow B_t), where ((B_t)_{t \geq 0}) is standard Brownian motion.

Definition 3.1. Standard Brownian motion ((B_t)_{t \geq 0}) is characterized by:

1. (B_0 = 0) almost surely.
2. Independent increments: For (0 \leq t_1 < t_2 \leq t_3 < t_4), (B_{t_2} - B_{t_1}) is independent of (B_{t_4} - B_{t_3}).
3. Gaussian increments: (B_t - B_s \sim \mathcal{N}(0, t - s)) for (0 \leq s < t).
4. Continuous paths: (t \mapsto B_t) is continuous almost surely.

3.1 The Reflection Principle

The reflection principle is one of the most elegant tools for computing probabilities involving Brownian motion and its running maximum.

Theorem 3.1 (Reflection Principle). For (a > 0) and (b \leq a),

[ \mathbb{P}\left(\max_{0 \leq s \leq t} B_s \geq a, B_t \leq b\right) = \mathbb{P}(B_t \geq 2a - b) ]

Proof sketch. Consider the first hitting time (\tau_a = \inf{s \geq 0 : B_s = a}). By the strong Markov property, the process after (\tau_a),

[ \tilde{B}s = B{\tau_a + s} - B_{\tau_a} ]

is again a Brownian motion, independent of the past. Reflect the path after (\tau_a): define (B’_s = B_s) for (s \leq \tau_a) and (B’s = 2a - B_s) for (s > \tau_a). By symmetry of Brownian motion, (B’) is also a Brownian motion. The event ({\max{s \leq t} B_s \geq a, B_t \leq b}) for (B) corresponds to ({B’_t \geq 2a - b}) for (B’). By equality in distribution, the probabilities are equal. ∎

Corollary 3.2 (Distribution of the Running Maximum).

[ \mathbb{P}\left(\max_{0 \leq s \leq t} B_s \geq a\right) = 2\mathbb{P}(B_t \geq a) = 2\left(1 - \Phi\left(\frac{a}{\sqrt{t}}\right)\right) ]

where (\Phi) is the standard normal CDF.

3.2 Brownian Motion with Drift and the Girsanov Theorem

For queueing applications, we often need Brownian motion with drift: (X_t = \mu t + \sigma B_t), representing the net input rate ((\mu) > 0 for input exceeding output, growing queue). The Girsanov theorem provides a change of measure that removes the drift, enabling exact computations via the reflection principle.

Theorem 3.3 (Girsanov). Let (X_t = B_t + \int_0^t \theta_s ds) under (\mathbb{P}). Then under (\mathbb{Q}) defined by the Radon-Nikodym derivative (\frac{d\mathbb{Q}}{d\mathbb{P}} = \exp\left(-\int_0^T \theta_s dB_s - \frac{1}{2}\int_0^T \theta_s^2 ds\right)), the process (X_t) is a standard Brownian motion.

For constant drift (\theta_s \equiv \theta), this gives (\mathbb{Q}(A) = \mathbb{E}_{\mathbb{P}}[e^{-\theta B_T - \theta^2 T/2} \mathbf{1}_A]), and the reflected Brownian motion computations carry over.

4. Queueing Theory

4.1 The M/M/1 Queue

The M/M/1 queue is the simplest queueing model: Poisson arrivals at rate (\lambda), exponential service times at rate (\mu), single server, infinite buffer. The number in system (N(t)) is a birth-death process with stationary distribution (\pi(n) = (1 - \rho)\rho^n).

Key performance measures (at stationarity):

• Mean number in system: (\mathbb{E}[N] = \rho / (1 - \rho)).
• Mean sojourn time (Little’s Law): (\mathbb{E}[W] = \mathbb{E}[N] / \lambda = 1 / (\mu - \lambda)).
• Tail probability: (\mathbb{P}(N \geq k) = \rho^k).

As (\rho \to 1) (heavy traffic), (\mathbb{E}[N] \to \infty) with a pole at (\rho = 1). This is the critical phenomenon: near saturation, small changes in load cause large changes in performance.

4.2 Jackson Networks

A Jackson network is a network of (K) M/M/1-like queues with probabilistic routing. External arrivals to node (i) are Poisson with rate (\alpha_i). After service at node (i), a job is routed to node (j) with probability (p_{ij}) or leaves the system with probability (1 - \sum_j p_{ij}).

Jackson’s theorem (1963): the stationary distribution of the network is product-form:

[ \pi(n_1, \ldots, n_K) = \prod_{i=1}^K (1 - \rho_i) \rho_i^{n_i} ]

where (\rho_i = \lambda_i / \mu_i) and the effective arrival rates (\lambda_i) satisfy the traffic equations:

[ \lambda_i = \alpha_i + \sum_{j=1}^K \lambda_j p_{ji} ]

This is a remarkable result: the nodes behave independently at stationarity, despite the complex routing. The product-form property enables efficient analysis of large distributed systems and is the foundation for capacity planning in cloud computing.

5. Reliability and Failure Models

Stochastic processes provide the mathematical foundation for reliability engineering.

5.1 The Failure Rate Function

The failure rate (or hazard rate) of a nonnegative random variable (T) with density (f) and CDF (F) is:

[ h(t) = \frac{f(t)}{1 - F(t)} = -\frac{d}{dt} \log(1 - F(t)) ]

Intuitively, (h(t) dt = \mathbb{P}(t \leq T < t + dt \mid T \geq t)). The exponential distribution has constant failure rate (h(t) = \lambda). The Weibull distribution ((h(t) = \alpha \beta t^{\beta-1})) captures increasing ((\beta > 1), aging/wear-out) or decreasing ((\beta < 1), infant mortality) failure rates.

5.2 Renewal Processes and Alternating Renewal

A renewal process counts the number of renewals (e.g., repairs) of a system that fails and is repaired repeatedly. If the time to failure is (T_i) and the repair time is (R_i), the availability (fraction of time the system is operational) is:

[ A = \frac{\mathbb{E}[T]}{\mathbb{E}[T] + \mathbb{E}[R]} ]

by the renewal reward theorem. For a system with multiple components, Markov models (CTMCs) capture the state-dependent failure and repair behavior, enabling the computation of system-level availability and mean time to failure (MTTF).

6. Martingales and Stopping Times: Optional Stopping in Queueing

A martingale is a stochastic process ((M_t){t \geq 0}) satisfying (\mathbb{E}[|M_t|] < \infty) and (\mathbb{E}[M{t+s} \mid \mathcal{F}_t] = M_t) for all (s, t \geq 0). The martingale property captures the notion of a “fair game”: your expected future fortune, given the past, equals your current fortune. Brownian motion is a martingale. So is the compensated Poisson process (N(t) - \lambda t).

Theorem 6.1 (Optional Stopping Theorem). Let (M_t) be a martingale and (\tau) a bounded stopping time. Then (\mathbb{E}[M_\tau] = \mathbb{E}[M_0]). If (\tau) is unbounded but (M_{t \wedge \tau}) is uniformly integrable, the same conclusion holds.

Application: Queue Busy Periods. Consider an M/G/1 queue. The workload process (W_t) (remaining work in the system) satisfies:

[ W_t = W_0 + \sum_{i=1}^{A(t)} S_i - t ]

where (A(t)) is the number of arrivals by time (t) and (S_i) are service times. The process (W_t - (\rho - 1)t) is a martingale (where (\rho = \lambda \mathbb{E}[S])). For the busy period duration (B) (time until (W_t) hits 0), optional stopping gives:

[ \mathbb{E}[B] = \frac{\mathbb{E}[W_0]}{1 - \rho} ]

This elegant formula is the foundation for analyzing task completion times, server utilization, and buffer dynamics in storage systems.

6.1 Wald’s Equation and Its Martingale Proof

Wald’s equation states that for i.i.d. (X_i) with finite mean and a stopping time (\tau) with finite expectation that is independent of the future,

[ \mathbb{E}\left[\sum_{i=1}^\tau X_i\right] = \mathbb{E}[\tau] \cdot \mathbb{E}[X_1] ]

This follows from the optional stopping theorem applied to the martingale (M_n = \sum_{i=1}^n (X_i - \mu)). Wald’s equation is fundamental to the analysis of sequential algorithms, randomized data structures (skip lists, treaps), and the expected runtime of Las Vegas algorithms.

7. Stochastic Differential Equations and Itô Calculus for Network Dynamics

Brownian motion is nondifferentiable almost everywhere, so ordinary calculus fails. Itô calculus provides the correct framework for integrating with respect to Brownian motion and solving stochastic differential equations (SDEs).

Definition 7.1 (Itô Integral). The Itô integral of a process (H_t) with respect to Brownian motion is defined as the (L^2) limit:

[ \int_0^T H_s dB_s = \lim_{|\Pi| \to 0} \sum_{i} H_{t_i} (B_{t_{i+1}} - B_{t_i}) ]

Crucially, the integrand is evaluated at the left endpoint of each interval—the Itô convention. This makes the integral a martingale and gives it zero expectation: (\mathbb{E}[\int H_s dB_s] = 0).

Theorem 7.1 (Itô’s Lemma). For (f \in C^2) and (dX_t = \mu_t dt + \sigma_t dB_t),

[ df(X_t) = \left(f’(X_t)\mu_t + \frac{1}{2} f’’(X_t)\sigma_t^2\right)dt + f’(X_t)\sigma_t dB_t ]

The correction term (\frac{1}{2} f’’(X_t)\sigma_t^2) is the “Itô correction”—absent in ordinary calculus, essential for stochastic processes.

7.1 Geometric Brownian Motion and TCP Throughput Models

Geometric Brownian motion (GBM): (dS_t = \mu S_t dt + \sigma S_t dB_t). Solution by Itô’s lemma:

[ S_t = S_0 \exp\left((\mu - \sigma^2/2)t + \sigma B_t\right) ]

GBM models congestion window evolution in TCP under random loss (Mathis et al., 1997). The TCP sending rate (X_t) satisfies roughly (dX_t = (1/RTT^2) dt - (X_t/2) dN_t), where (N_t) is a loss-indication Poisson process. In the heavy-traffic, many-flows limit, the aggregate throughput follows a reflected GBM, and the stationary distribution gives the throughput formula:

[ \mathbb{E}[\text{Throughput}] \approx \frac{MSS}{RTT} \cdot \frac{1}{\sqrt{p}} ]

where (p) is the packet loss probability—a classic result derived via stochastic calculus.

8. Large Deviations of Queueing Systems: The Effective Bandwidth Principle

The effective bandwidth of a traffic source with cumulative arrivals (A(t)) is:

[ \alpha(\theta) = \lim_{t \to \infty} \frac{1}{\theta t} \log \mathbb{E}[e^{\theta A(t)}] ]

This function encapsulates the burstiness characteristics relevant to buffer overflow. For a Poisson process of rate (\lambda), (\alpha(\theta) = \lambda(e^\theta - 1)/\theta). For Gaussian traffic (fractional Brownian motion with Hurst parameter (H)), (\alpha(\theta) = m + \theta \sigma^2 t^{2H-1}/2).

Theorem 8.1 (Gibson & Sen, 1991). For a queue fed by (N) independent sources, each with effective bandwidth (\alpha(\theta)), served at rate (C), the overflow probability satisfies:

[ \lim_{B \to \infty} \frac{1}{B} \log \mathbb{P}(Q \geq B) = -\theta^* ]

where (\theta^) solves (\alpha(\theta^) = C) (if all sources are identical and always active). This principle—the effective bandwidth approximation—is the foundation for Connection Admission Control (CAC) in ATM and MPLS networks, and for buffer dimensioning in data center switches.

8.1 Many-Sources Asymptotic and the Bahadur-Rao Refinement

For (N) sources sharing a buffer of size (B = N b), the many-sources asymptotic (Botvich and Duffield, 1995) gives:

[ \lim_{N \to \infty} \frac{1}{N} \log \mathbb{P}(Q_N \geq N b) = -I(b) ]

where (I(b)) is a rate function obtained via the Gärtner-Ellis theorem applied to the empirical arrival process. The Bahadur-Rao refinement adds a (\Theta(1/\sqrt{N})) polynomial prefactor for accurate dimensioning at moderate scale.

9. Markov-Modulated Processes and the Matrix-Geometric Method

Many real systems exhibit burstiness that a simple Poisson model cannot capture. A Markov-modulated Poisson process (MMPP) has an underlying CTMC (the “environment”) modulating the arrival rate. When the environment is in state (i), arrivals follow a Poisson process with rate (\lambda_i).

Definition 9.1 (Quasi-Birth-Death Process). A QBD is a CTMC on ({(i, j) : i \geq 0, 1 \leq j \leq m}) where the “level” (i) changes by at most (\pm 1) per transition. The generator is block-tridiagonal:

[ Q = \begin{pmatrix} B_0 & A_0 & 0 & \cdots \ A_2 & A_1 & A_0 & \cdots \ 0 & A_2 & A_1 & \cdots \ \vdots & \vdots & \vdots & \ddots \end{pmatrix} ]

Theorem 9.1 (Matrix-Geometric Solution, Neuts, 1981). If the QBD is positive recurrent, the stationary distribution has the matrix-geometric form:

[ \pi_i = \pi_0 R^i, \quad i \geq 0 ]

where (R) (the rate matrix) is the minimal nonnegative solution to the matrix quadratic equation:

[ A_0 + R A_1 + R^2 A_2 = 0 ]

The matrix (R) can be computed efficiently via iterative schemes (logarithmic reduction, cyclic reduction).

The MMPP/M/1 queue is a QBD, and its stationary queue length distribution is matrix-geometric—a tractable generalization of the M/M/1 geometric distribution that captures burstiness through the environmental Markov chain. This framework has been applied to model wireless channel fading, disk I/O with seek delays, and cloud VM provisioning under workload spikes.

MMPP/M/1 queue block structure:

        ┌─────────────────────┐
Level:  0  1  2  3  ...
        │  B0 A0              │
     Q =│  A2 A1 A0           │
        │     A2 A1 A0        │
        │        A2 A1 ...    │
        └─────────────────────┘

Solution: π_i = π_0 R^i, where R solves A_0 + R A_1 + R^2 A_2 = 0

11. Stochastic Optimal Control and the Hamilton-Jacobi-Bellman Equation

Many problems in computer systems involve making decisions under uncertainty: admission control for queues, dynamic power management, and rate adaptation in wireless networks. Stochastic optimal control provides the mathematical framework for these problems, with the Hamilton-Jacobi-Bellman (HJB) equation as its centerpiece.

11.1 The Controlled Diffusion Process

Consider a system whose state X_t evolves according to a controlled stochastic differential equation:

dX_t = b(X_t, u_t) dt + sigma(X_t) dB_t

where u_t is a control process chosen from an admissible set U. The objective is to minimize (or maximize) an expected cost functional:

J(x, u) = E [ int_0^T f(X_t, u_t) dt + g(X_T) | X_0 = x ]

where f is the running cost and g is the terminal cost.

Theorem 11.1 (Hamilton-Jacobi-Bellman Equation). Let V(x, t) = inf_u J(x, u) be the value function (optimal cost-to-go). If V is sufficiently smooth, then V satisfies the HJB partial differential equation:

-V_t = inf_{u in U} [ f(x, u) + b(x, u) . nabla V + (1/2) Tr(sigma(x) sigma(x)^T nabla^2 V) ]

with terminal condition V(x, T) = g(x). The optimal control at state (x, t) is the argmin of the Hamiltonian:

u*(x, t) = argmin_{u in U} [ f(x, u) + b(x, u) . nabla V(x, t) ]

Proof sketch. The proof uses the dynamic programming principle: V(x, t) = inf_u E[ int_t^{t+dt} f ds + V(X_{t+dt}, t+dt) ]. Expanding V(X_{t+dt}, t+dt) via Ito’s lemma and taking the limit dt -> 0 yields the HJB equation. The verification theorem ensures that any smooth solution of the HJB equation equals the value function. This is the stochastic analog of the Bellman equation for discrete-time MDPs. ∎

11.2 Application: Optimal Admission Control in Queueing Systems

Consider a single server queue with Poisson arrivals (rate lambda) and exponential service (rate mu). The controller can accept or reject each arriving job. Accepted jobs pay a reward R; rejected jobs incur a penalty C. The state is the queue length Q_t, and the control u_t in {accept, reject} is chosen at each arrival.

The HJB equation for the value function V(q) (in the infinite-horizon discounted case with discount rate beta) becomes:

beta V(q) = lambda max{ R + V(min(q+1, B)), V(q) } + mu(V(max(q-1, 0)) - V(q))

for a buffer of capacity B. This is a system of linear equations in the unknown V(0), …, V(B), solvable by value iteration or policy iteration. The optimal policy has a threshold structure: accept if q < q* and reject if q >= q*, for some optimal threshold q* determined by the relative values of R, C, and the traffic intensity rho = lambda/mu.

11.3 The Linear-Quadratic-Gaussian (LQG) Regulator

When the dynamics are linear and the cost is quadratic, the HJB equation has a closed-form solution. For:

dX_t = (A X_t + B u_t) dt + sigma dB_t J = E[ int_0^T (X_t^T Q X_t + u_t^T R u_t) dt + X_T^T S X_T ]

the value function is quadratic: V(x, t) = x^T P(t) x + c(t), where P(t) solves the Riccati differential equation:

-P_t = A^T P + P A - P B R^{-1} B^T P + Q

with P(T) = S. The optimal control is linear state feedback: u*(x) = -R^{-1} B^T P(t) x. This is the foundation of linear-quadratic control in continuous time, with applications in congestion control (TCP Vegas/AIMD as approximate LQG) and in thermal management of data centers.

12. Levy Processes and Heavy-Tailed Phenomena in Computer Systems

While Poisson processes and Brownian motion capture many system behaviors, modern computer workloads often exhibit heavy-tailed characteristics: file sizes on web servers, CPU bursts in operating systems, and inter-arrival times of packets in self-similar network traffic. Levy processes generalize Brownian motion to include jumps, providing a flexible modelling framework.

12.1 The Levy-Khintchine Representation

Definition 12.1 (Levy Process). A Levy process X_t is a stochastically continuous process with stationary independent increments and X_0 = 0. Every Levy process can be decomposed as the sum of a deterministic drift, a Brownian motion, and a pure-jump process:

X_t = gamma t + sigma B_t + int_{|x|<1} x (N_t(dx) - t nu(dx)) + sum_{s <= t} Delta X_s 1_{|Delta X_s| >= 1}

Theorem 12.1 (Levy-Khintchine Formula). The characteristic function of a Levy process is:

E[e^{i theta X_t}] = exp( t [ i gamma theta - (1/2) sigma^2 theta^2 + int_{R{0}} (e^{i theta x} - 1 - i theta x 1_{|x|<1}) nu(dx) ] )

where (gamma, sigma^2, nu) is the Levy triplet: gamma is the drift, sigma^2 the Brownian variance, and nu the Levy measure encoding the intensity of jumps of different sizes.

12.2 Alpha-Stable Processes and Self-Similar Traffic

An alpha-stable Levy process has the scaling property X_{ct} = c^{1/alpha} X_t in distribution. When alpha < 2, the variance is infinite (heavy tails), and the process exhibits self-similarity: aggregated traffic looks statistically identical at different timescales. This is precisely the behavior observed in Ethernet LAN traffic (Leland, Taqqu, Willinger, and Wilson, 1994).

The parameter alpha controls the tail heaviness: P(|X_t| > x) ~ C x^{-alpha} as x -> infinity. For network traffic, empirical measurements give alpha in the range 1.1 to 1.7, indicating infinite variance and the breakdown of classical queueing formulas based on exponential decay of buffer overflow probabilities. Instead, buffer overflow probabilities decay polynomially: P(Q > B) ~ B^{-alpha+1}, much slower than the exponential decay of light-tailed systems.

12.3 Subordinators and Time-Changed Processes

A subordinator is a non-decreasing Levy process (e.g., the Gamma process, inverse Gaussian process). Subordinators can be used to time-change Brownian motion to create stochastic volatility models:

Y_t = mu * S_t + sigma B_{S_t}

where S_t is a subordinator. For the variance gamma process (Madan and Seneta, 1990), S_t is a Gamma process, yielding a model with both skewness and excess kurtosis. In computer systems, time-changed Brownian motion models CPU utilization where the “business time” S_t represents the cumulative CPU demand, which itself evolves randomly due to varying workload intensity.

13. Point Processes, Palm Calculus, and the Feller Coupling

Beyond Poisson processes, the general theory of point processes provides a unified framework for modelling discrete events in continuous time: packet arrivals, page faults, sensor readings, and user session starts.

13.1 The Campbell-Mecke Formula

Definition 13.1 (Point Process). A point process Phi on a space E is a random counting measure: Phi = sum_i delta_{X_i} where X_i are random points in E. The intensity measure Lambda(A) = E[Phi(A)] gives the expected number of points in A.

Theorem 13.1 (Campbell-Mecke Formula). For a point process Phi with intensity measure Lambda and any non-negative measurable function f:

E[ sum_{x in Phi} f(x, Phi) ] = int_E E[ f(x, Phi + delta_x) ] Lambda(dx)

This is the Palm calculus version of the Campbell formula. It decomposes the expectation over the point process into an integral over the intensity measure, where the integrand is the expected value of f when a point is artificially added at x. This is the continuous-time analogue of the “law of total probability for random sums” and is fundamental for the analysis of queueing networks with general arrival processes.

13.2 The Feller Coupling and the Poisson-Dirichlet Distribution

The Feller coupling provides a remarkable connection between the Poisson process and combinatorial structures. Consider n i.i.d. uniform points on [0,1]. Their spacings (ordered differences) follow the Dirichlet distribution. In the limit n -> infinity, viewed through the lens of a Poisson process on [0, infinity) with rate 1/n -> 0, the normalized spacings converge to the Poisson-Dirichlet distribution, which governs the asymptotic distribution of cycle lengths in random permutations, the sizes of components in the random graph near the critical window, and the allele frequencies in population genetics (Ewens sampling formula).

Theorem 13.2 (GEM Distribution via Size-Biased Permutation). Let (V_1, V_2, …) be i.i.d. Beta(1, theta) random variables, and define P_1 = V_1, P_k = V_k prod_{i=1}^{k-1} (1 - V_i). Then (P_1, P_2, …) has the Griffiths-Engen-McCloskey (GEM) distribution, which is a size-biased permutation of the Poisson-Dirichlet distribution with parameter theta. The GEM distribution arises naturally from the stick-breaking construction of the Dirichlet process, which is fundamental to Bayesian nonparametrics and the Chinese Restaurant Process for clustering.

13.3 Application: Cache Replacement and the Working Set Model

Consider a cache that stores items requested by a point process. The working set at time t is the set of distinct items requested in [t - T, t], where T is the window size. Under a Poisson request process with rate lambda and item popularity following a Zipf distribution, the expected working set size is:

E[ W(t) ] = sum_{i=1}^{infinity} (1 - exp(-lambda p_i T))

where p_i is the probability of requesting item i. This formula, derived via the Campbell-Mecke formula applied to the Poisson process of requests for each item, is the theoretical basis for LRU (Least Recently Used) cache performance analysis and for the design of the Adaptive Replacement Cache (ARC) algorithm. The analysis extends to time-varying request rates via the Cox process (doubly stochastic Poisson process), where lambda(t) is itself a stochastic process.

13.4 Determinantal Point Processes and Diversity in Recommendation Systems

A determinantal point process (DPP) is a point process where the correlation functions are given by determinants of a kernel matrix K. DPPs exhibit repulsion: points tend not to cluster together. This makes them ideal for modelling diversity in recommendation systems: given a kernel K(x, y) measuring similarity between items x and y, a DPP with kernel K selects a subset of items that are simultaneously relevant (high K(x, x)) and diverse (low K(x, y) for y != x).

The probability of selecting a subset S is proportional to det(K_S), where K_S is the submatrix of K indexed by S. Computing the MAP (maximum a posteriori) subset under a DPP is NP-hard in general, but for certain kernels (e.g., when K is the reproducing kernel of a low-dimensional feature space), it reduces to a determinant maximization problem solvable by greedy algorithms with approximation guarantees. This connects stochastic geometry (point processes) to modern machine learning (recommendation diversity).

The theory of point processes and Palm calculus bridges the gap between the microscopic randomness of individual events and the macroscopic behavior of systems. For the computer scientist, these tools provide rigorous methods for analyzing the performance of caches, the reliability of distributed protocols, and the diversity of recommendation outputs. The Campbell-Mecke formula, in particular, is the continuous-time analogue of the law of large numbers for random sums and is indispensable for any analysis involving random measures in continuous time.

13.5 The Fermat-Weber Problem and Spatial Poisson Processes in Facility Location

Consider the problem of placing K servers to minimize the expected distance to the nearest server for users distributed according to a spatial Poisson process with intensity lambda(x) over a region D. This is the continuous K-median problem. If the users form a homogeneous Poisson process of intensity lambda, the optimal server locations are the K-medians of the region, and the expected cost (average distance) scales as Theta(1/sqrt(K)) for planar regions.

More generally, for a Poisson process with intensity lambda(x), the expected distance to the k-th nearest point is:

E[ R_k ] = Gamma(k + 2/d) / Gamma(k) * (1 / (lambda * V_d))^{1/d}

where V_d is the volume of the unit ball in d dimensions and Gamma is the gamma function. This formula governs the performance of nearest-neighbor search in spatial databases, the connectivity range in wireless sensor networks, and the expected latency in geographically distributed content delivery networks. The analysis uses the fact that for a homogeneous Poisson process, the number of points in a ball of radius r is Poisson distributed with mean lambda * V_d * r^d, and the distance to the k-th nearest neighbor follows a generalized Gamma distribution.

These spatial Poisson process models form the mathematical backbone of facility location theory, wireless network planning, and the analysis of spatial databases—demonstrating once again how the simplest stochastic process yields deep insights into complex system design problems.

The interplay between Poisson-driven spatial randomness and deterministic optimization creates a rich mathematical landscape with direct practical implications for infrastructure design at scale.

14. Summary

Stochastic processes—Poisson, Markov, Brownian—are the mathematical language of randomness in time. The Poisson process models arrivals. CTMCs model state transitions. Brownian motion models heavy-traffic limits and provides exact hitting-time distributions via the reflection principle. Queueing theory translates these processes into performance measures: latency, throughput, buffer occupancy. Reliability theory translates them into availability, MTTF, and failure rates.

For the computer scientist building distributed systems, these tools are indispensable. They enable capacity planning (how many servers do we need to meet an SLO?), tail latency analysis (what is the 99.9th percentile response time?), and reliability engineering (what is the probability of data loss in a given year?). The mathematics of stochastic processes is not just elegant—it is the engine of evidence-based system design.

To go deeper, Ross’s Stochastic Processes is the standard reference. Kleinrock’s Queueing Systems is the classic two-volume treatment. Harrison’s Brownian Motion and Stochastic Flow Systems develops the heavy-traffic theory. And for reliability, Trivedi’s Probability and Statistics with Reliability, Queuing, and Computer Science Applications bridges theory and practice.