Arrival Architect Slides Arrival Architect
Poisson Process: Axioms and Simulation
Graduate Level Continuous-Time Models
The Essential Question
"If we simulate the arrival of emails to a server, does the timing of the last email affect how long we wait for the next one?"
Theory
Independence and memorylessness.
Computation
Generating random interarrival times.
The Axiomatic Definition
Axiom 1: Independent Increments
Numbers of events in non-overlapping time intervals are independent random variables.
Axiom 2: Stationary Increments
The distribution of events in an interval depends only on the length of the interval, not its position.
Axiom 3: Rare Events
As \( h \to 0 \), \( P(N(h)=1) = \lambda h + o(h) \) and \( P(N(h) \ge 2) = o(h) \).
Simulating Arrivals
Method A: Interarrival Gaps
The time between events \( X_i \) follows an Exponential distribution with rate \( \lambda \).
interarrivals = rexp(n, rate = lambda)
Method B: Time-Slicing
Divide time \( T \) into \( n \) small intervals. In each, an event occurs with probability \( p = \lambda \frac{T}{n} \).
# Binomial Approximation
The Memoryless Property
\[ P(X > s + t \mid X > s) = P(X > t) \]
Crucial implication for simulation: The time until the next event is independent of the time elapsed since the last event. In our Poisson world, the past is irrelevant to the future waiting time.
Does the real world work like this? Consider bus arrivals vs. email arrivals.
Arrival Architect Lab Arrival Architect Lab
Lesson 1: Poisson Process Simulation
Student Name / Date
Objective
In this computational lab, you will validate the Poisson axioms through simulation. You will compare the two primary methods of generating Poisson events—interarrival gap generation and the binomial limit—and empirically test the memoryless property of the Exponential distribution.
1
The Interarrival Gap Method
Given a rate parameter \( \lambda = 5 \) (events per hour), simulate a Poisson process for \( T = 100 \) hours. Generate the sequence of interarrival times \( X_i \sim \text{Exp}(\lambda) \) and calculate the cumulative arrival times \( S_n = \sum_{i=1}^n X_i \).
A. Pseudocode / Logic
B. Simulation Result Sketch
[Sketch a plot of N(t) vs t]
Reflect: How do you determine the total number of events \( n \) to generate if you want to cover a fixed time \( T \)?
2
The Binomial Limit
Divide the interval \( [0, 1] \) into \( n \) sub-intervals. Let \( Y_i \sim \text{Bernoulli}(p) \) where \( p = \lambda / n \).
# Parameters: lambda = 10, T = 1
What happens to the distribution of the total count \( N(1) = \sum Y_i \) as \( n \to \infty \) while holding \( \lambda \) constant? Why does this confirm the Poisson distribution of \( N(t) \)?
3
The Memoryless Empirical Test
"Given that we have already waited 10 minutes for an email, what is the probability we wait another 5?"
Write a simulation script to generate \( 10,000 \) samples from \( X \sim \text{Exp}(\lambda=0.5) \).
Calculate the mean of all samples \( X \).
Filter the samples to keep only those where \( X > 2 \). Subtract 2 from these remaining samples (i.e., \( X_{new} = X - 2 \mid X > 2 \)).
Calculate the mean of the new samples \( X_{new} \).
Expected Theoretical Result
Simulated Result Hypothesis
Synthesis
If the interarrival times are independent and exponentially distributed, explain why the number of events in non-overlapping intervals must be independent (Axiom 1).
Computational Stochastic Processes - Advanced Topics in Statistics
Stochastic Foundations Teacher Guide Teacher Implementation Guide
Lesson 1: Arrival Axioms
Lesson Time 90 Min
Instructional Goal
Students will move beyond simple counting distributions to understanding the temporal structure of the Poisson process. The bridge between the discrete Binomial and the continuous Poisson is the critical conceptual hurdle here.
Key Mathematical Notation
\( N(t) \): Number of arrivals in \( [0, t] \)
\( X_i \): \( i^{th} \) interarrival time, \( X_i \sim \text{Exp}(\lambda) \)
\( S_n = \sum_{j=1}^n X_j \): Arrival time of the \( n^{th} \) event
Required Toolkit
R, Python, or Julia
`rexp`, `rpois`, `rbinom`
Cumulative Sum functions
Pacing & Facilitation
Phase Mins Teacher Actions & Discussion Prompts Hook & Theory 20 Present the "Email Server" scenario. Ask: "Does waiting longer make an arrival more likely?" Contrast Poisson with periodic systems. Axiom Deep Dive 25 Walk through the three axioms. Focus heavily on Axiom 3 (\( o(h) \)). Show how Axiom 3 prevents "simultaneous" arrivals. Lab: Part 1 & 2 30 Students code interarrival cumulative sums. Circulate and check: Are they using \( \lambda \) or \( 1/\lambda \) correctly in their code? Debrief 15 Discuss the Memoryless empirical result. Ask: "Is it counter-intuitive?" Discuss how it simplifies modeling vs. real-world complexity.
Common Graduate Pitfalls
1. Parametrization Confusion
In R, `rexp` takes `rate = lambda`. In some Python libraries (like SciPy), it may take `scale = 1/lambda`. Students will often get this inverted.
2. Boundary Effects
When simulating over a fixed time \( T \), the last interarrival usually exceeds \( T \). Students must decide if they truncate or only count arrivals up to \( T \).
Extension Questions for High-Performers
"How would you modify the interarrival simulation if the rate \( \lambda \) was a function of time \( \lambda(t) \)?" (Leads to Non-Homogeneous Poisson Processes).
"Can you prove that \( S_n \) follows a Gamma distribution using the sum of exponentials property?"
Exponential Risks Slides Exponential Risks
Memoryless Property & Competing Risks
Lesson 2 Reliability Theory
The "Holding Time"
In continuous-time models, the time a system spends in a specific state is the holding time.
Why must it be Exponential?
If the future of the process only depends on the current state (Markov Property), then the distribution of the remaining time in that state must be independent of how long we've already been there.
Competing Risks
The Setup
Imagine multiple independent failure processes happening simultaneously:
\( X_1 \sim \text{Exp}(\lambda_1) \)
\( X_2 \sim \text{Exp}(\lambda_2) \)
\( X_n \sim \text{Exp}(\lambda_n) \)
The Race
\( T = \min(X_1, X_2, \dots, X_n) \)
The Winning Rate
\[ \lambda_{total} = \sum_{i=1}^n \lambda_i \]
Who Fails First?
What is the probability that component \( j \) is the one that triggers the transition?
The Formula
\[ P(X_j = \min_i X_i) = \frac{\lambda_j}{\sum_{k=1}^n \lambda_k} \]
"Probability is proportional to the rate."
Visualizing Failure Rates
Constant Hazard Rate
Unlike humans, components in this model do not "wear out." They are just as likely to fail at age 100 as they are at age 1.
Simulation Check
If we simulate 10,000 components, the histogram of failure times should always look the same, even if we truncate the first \( T \) years.
Declining PDF vs. Constant Hazard
Race to Failure Problem Set Race to Failure
Problem Set: Competing Risks & Reliability
Student Identifier
Theoretical Context
In system reliability, a system "fails" when any one of its critical components fails. If components have independent exponential lifetimes with rates \( \lambda_1, \lambda_2, \dots, \lambda_n \), the system lifetime follows \( \text{Exp}(\sum \lambda_i) \). This "race" determines which state transition occurs in a Continuous-Time Markov Chain.
01 The Server Trio
A server cluster has three independent hardware components that can fail. Component A has an MTBF (Mean Time Between Failures) of 1,000 hours. Component B has an MTBF of 2,500 hours. Component C has an MTBF of 5,000 hours.
A. Calculate the failure rate \( \lambda \) for each component (failures per hour).
B. What is the probability that the cluster survives at least 500 hours?
C. Calculate the probability that Component B is the first to fail.
02 Parallel Redundancy
Consider a system where two identical components work in parallel. The system only fails when both have failed. Each component has a failure rate \( \lambda = 0.01 \).
Explain why the time until the first failure is Exponential, but the time until the system fails is NOT Exponential.
Hint: Consider the memoryless property and the state space of the system.
03 Empirical Verification
If you were to simulate the "Server Trio" 10,000 times, describe the steps you would take to verify your answer to 1C.
# Sketch code logic here...
Random Processes in Statistics Lesson 2: Reliability & Competing Risks
Race to Failure Mastery Key Mastery Key: Race to Failure
Instructor Solution Guide & Derivations
01
The Server Trio
A. Rate Calculations
MTBF is the reciprocal of the rate for an Exponential distribution.
B. Survival Probability
\( P(T > 500) = e^{-\lambda_{total} \times 500} \)
C. Prob(B fails first)
\( P(X_B = \min) = \frac{\lambda_B}{\lambda_A + \lambda_B + \lambda_C} = \frac{0.0004}{0.0016} = \frac{1}{4} = 0.25 \).
02
Parallel Redundancy
Conceptual Derivation
Let \( T \) be the time to system failure. \( T = \max(X_1, X_2) \).
03
Empirical Verification
# R Simulation Example
# Generate failures
# Find which was min in each row
Process Merging Slides Process Merging
Thinning and Superposition of Poisson Processes
Lesson 3 Stream Operations
Superposition
When we combine multiple independent Poisson streams into one, the result is still a Poisson process.
\( \lambda_{merged} = \sum_{i=1}^k \lambda_i \)
Example: Server Traffic
Mobile users arrive at 10 req/s.
Desktop users arrive at 5 req/s.
Total System Traffic: 15 req/s
Thinning (Decomposition)
If we split a Poisson process \( N(t) \) into two streams based on a probability \( p \), the resulting streams are independent Poisson processes.
Source
\( \lambda \)
Stream 1
\( \lambda p \)
Stream 2
\( \lambda(1-p) \)
The "Surprise" of Independence
Crucially, the resulting streams are statistically independent.
"Knowing that Stream 1 received a lot of arrivals in the last minute tells you absolutely nothing about how many Stream 2 received."
This only holds if the source is Poisson.
Coding Thinning
Algorithm
Simulate arrivals \( S_n \) with rate \( \lambda \).
For each arrival, flip a biased coin with \( P(\text{Heads}) = p \).
Assign arrival to Stream 1 if Heads, else Stream 2.
# R Snippet
n <- length(arrivals)
selector <- runif(n) < p
stream1 <- arrivals[selector]
stream2 <- arrivals[!selector]
Traffic Flow Analysis Activity Traffic Flow Analysis
Simulation Activity: Thinning & Superposition
Team Members / Lab Section
Scenario
A network router receives two types of packets: Data (\( \lambda_D = 40 \) per ms) and Voice (\( \lambda_V = 20 \) per ms). The router thins the Data stream, dropping \( 5\% \) of packets due to congestion. We want to investigate the statistical properties of the resulting merged and thinned streams.
A
The Superposition Test
Simulate the Data and Voice arrival streams independently for \( T = 1000 \) ms. Combine them into a single sorted stream of timestamps.
1. Theoretical Arrival Rate
Show calculation here
2. Interarrival Distribution
What specific distribution should the merged interarrival times follow?
3. Variance Comparison
In a Poisson process, \( E[N(t)] = Var[N(t)] \). Describe how you would verify this for the merged stream using your simulation results.
B
The Independence Paradox
Take your simulated Data stream (\( \lambda = 40 \)). For each packet, "drop" it with probability \( p = 0.05 \). This creates two streams: Accepted and Dropped .
# Logic Check
Provocation:
"If I tell you that 50 packets were dropped in the first 10ms, does that change your expectation of how many were accepted in that same 10ms?"
Verification Task:
Explain how you would use a correlation coefficient (\( \rho \)) on binned counts to prove or disprove independence between the Accepted and Dropped streams.
Computational Stochastic Processes - Module 3: Stream Operations
Stream Operations Teacher Guide Implementation Notes: Stream Operations
Lesson 3: Process Merging & Decomposing
Pedagogical Strategy
The "Independence Paradox" in thinning is the most frequent point of confusion. Students intuitively feel that if you have a fixed number of arrivals, knowing Stream A has many arrivals should mean Stream B has fewer.
The Key: Remind them that in a Poisson process, the total number of arrivals \( N(t) \) is itself a random variable. This variability allows both sub-streams to be independent.
Common Misconceptions
Misconception: Deterministic thinning (e.g., dropping every 2nd packet) also results in a Poisson process.
Reality: Deterministic thinning results in an Erlang-2 distribution for interarrivals, which is NOT Poisson (memory is introduced).
Workshop Guide: Traffic Flow Analysis
Phase Duration Discussion Prompts / Tips
Introduction 15m "Why do we care about merging? Think about a cloud load balancer receiving requests from millions of clients."
Lab Part A 30m Encourage students to use `density` plots to compare the histogram of merged interarrivals to the theoretical Exponential curve.
The Paradox 20m Ask students to calculate Cov(N1, N2). They will be shocked when it comes out to zero. Prepare a slide for the proof.
Closure 15m Summarize: Thinning = Bernoulli Choice. Superposition = Simple Addition of Rates.
The Proof Hint
If you need to show the independence of thinned processes to the class, use the joint probability:
\[ P(N_1=n, N_2=m) = P(N_1=n, N_2=m \mid N=n+m) P(N=n+m) \] \[ = \binom{n+m}{n} p^n (1-p)^m \times \frac{e^{-\lambda} \lambda^{n+m}}{(n+m)!} \] \[ = \left( \frac{e^{-\lambda p} (\lambda p)^n}{n!} \right) \left( \frac{e^{-\lambda(1-p)} (\lambda(1-p))^m}{m!} \right) \]
The factorization into two independent Poisson PMFs is the "smoking gun."
Matrix Metamorphosis Slides Matrix Metamorphosis
Continuous-Time Markov Chains & Q-Matrices
Lesson 4 Infinitesimal Generators
The Infinitesimal View
In Discrete Time (DTMC), we use a transition matrix P.
In Continuous Time (CTMC), transitions can happen at any moment. We define rates instead of probabilities.
"The probability of moving from \( i \) to \( j \) in a tiny window \( h \) is \( q_{ij}h + o(h) \)."
The Evolution Equation
\[ P(t) = e^{Qt} \]
Where \( Q \) is the Generator Matrix.
Anatomy of \( Q \)
The Rules
Off-diagonal entries \( q_{ij} \) are non-negative rates.
Row sums must be Zero .
Diagonal entries are \( q_{ii} = -\sum_{j \neq i} q_{ij} \).
A 3-State Example
\[ Q = \begin{pmatrix} -3 & 1 & 2 \\ 0 & -2 & 2 \\ 4 & 0 & -4 \end{pmatrix} \]
Note: All rows sum to 0.
Kolmogorov Differential Equations
How do transition probabilities change over time?
Backward Equation
\[ \frac{d}{dt} P(t) = Q P(t) \]
Conditioning on the first transition.
Forward Equation
\[ \frac{d}{dt} P(t) = P(t) Q \]
Conditioning on the last transition.
The Jump Chain
A CTMC is really just two things:
Exponential Holding Times
The "Embedded" Discrete Chain
To simulate: Sample how long you stay, then flip a coin to see where you jump.
Embedded Transition Probabilities
\[ p_{ij} = \frac{q_{ij}}{-q_{ii}} \]
Generator Workshop Handout Generator Workshop
Lesson 4: Building and Solving Q-Matrices
Student Name
State Duration
\( \text{Exp}(-q_{ii}) \)
Jump Prob
\( q_{ij} / |q_{ii}| \)
Steady State
\( \pi Q = 0 \)
01 Diagram to Matrix
A system has three states: Idle (1) , Processing (2) , and Error (3) . Transitions occur at the following rates:
Idle to Processing: \( \alpha = 5 \)
Processing to Idle (Success): \( \beta = 10 \)
Processing to Error (Failure): \( \gamma = 1 \)
Error to Idle (Reset): \( \delta = 0.5 \)
[Draw your state transition diagram here]
Construct the Generator Matrix \( Q \):
\( Q = \dots \)
02 Long-Run Behavior
Set up the system of linear equations to solve for the stationary distribution \( \pi = (\pi_1, \pi_2, \pi_3) \).
Note: Don't forget the normalization constraint \( \sum \pi_i = 1 \).
03 The Simulation Logic
Describe the Gillespie-style simulation steps for the system above. If the system is currently in the Processing state, what are the next two specific values you must sample to determine the next state and the time spent in Processing?
Sample 1: Holding Time
Distribution & Parameter:
Sample 2: Destination State
Probabilities:
Stochastic Lab Generator Matrices
Generator Mastery Key Mastery Key: Generator Workshop
Detailed Derivations for Lesson 4
01
Matrix Construction
The rates are: \( 1 \to 2 \) at 5; \( 2 \to 1 \) at 10; \( 2 \to 3 \) at 1; \( 3 \to 1 \) at 0.5.
\[ Q = \begin{pmatrix} -5 & 5 & 0 \\ 10 & -11 & 1 \\ 0.5 & 0 & -0.5 \end{pmatrix} \]
02
Stationary Equations
From \( \pi Q = 0 \):
1. \( -5\pi_1 + 10\pi_2 + 0.5\pi_3 = 0 \)
Solving these (using Eq 3: \( \pi_2 = 0.5\pi_3 \), Eq 2: \( \pi_1 = \frac{11}{5}\pi_2 = 1.1\pi_3 \)) results in proportions of the time spent in each state.
03
Gillespie Logic
Sample 1: Holding Time
Sample \( \Delta t \sim \text{Exp}(\lambda=11) \). This is because the total rate leaving State 2 is \( 10 + 1 = 11 \).
Sample 2: Jump Direction
Jump to State 1 with \( P = 10/11 \).
Instructor Use Only - Do Not Distribute
Queue Dynamics Slides Queue Dynamics
Birth-Death Processes & Queuing Theory
Lesson 5 Stability Analysis
The Birth-Death Model
A special CTMC where transitions only occur between adjacent states \( n \to n+1 \) or \( n \to n-1 \).
Birth Rate (\( \lambda_n \)): Rate of arrival when system has \( n \) items.
Death Rate (\( \mu_n \)): Rate of service when system has \( n \) items.
n-1
lambda mu
n
lambda mu
n+1
The M/M/1 Archetype
Stability Condition
The system is stable (doesn't grow to infinity) if and only if the arrival rate is less than the service rate.
\[ \rho = \frac{\lambda}{\mu} < 1 \]
Key Metrics
Avg Number in System: \( L = \frac{\rho}{1-\rho} \)
Avg Waiting Time: \( W = \frac{1}{\mu - \lambda} \)
As \( \rho \to 1 \), wait times explode exponentially.
Simulating the Bottleneck
The "Unstable" Regime
When \( \rho \ge 1 \), the queue length follows a Random Walk with Drift . It never reaches a steady state.
Observation:
In simulation, you will see the queue length trend upward linearly with time.
# Simulation loop
while time < T:
rate = lambda + mu if q > 0 else lambda
dt = rexp(1, rate)
if runif(1) < lambda/rate:
q += 1
else:
q -= 1
Little's Law
\[ L = \lambda W \]
"The average number of customers in a system is equal to the average arrival rate multiplied by the average time spent in the system."
This holds for almost ANY queuing system, regardless of distributions!
Queue Crashers Case Study Queue Crashers
Case Study: Stability & Performance
Student Name / Group ID
In this case study, you will analyze a technical support helpdesk modeled as an M/M/1 queue. Customers call at a rate of \( \lambda = 15 \) calls per hour. The average technician takes 3 minutes to resolve a call.
1 Steady State Baseline
A. Traffic Intensity (\( \rho \))
Calculate the utilization of the technician.
B. Average Wait Time (\( W \))
What is the average time a caller spends on the phone (waiting + service)?
2 Stability Breakdown
Management is considering a marketing campaign that will increase the arrival rate to \( \lambda = 21 \) calls per hour.
Question:
What happens to the queue stability? Describe the expected behavior of the system over a 24-hour simulation run at this new rate.
[Enter qualitative analysis here...]
3 The Multi-Server Solution
To handle the new rate (\( \lambda = 21 \)), the company adds a second technician (\( k=2 \)). Both have the same service rate \( \mu \).
Birth-Death State Space:
Write the service rate \( \mu_n \) as a function of \( n \).
\( \mu_n = \begin{cases} \dots \end{cases} \)
Simulated Insight:
How does having two servers change the idle probability (\( \pi_0 \)) compared to one faster server with rate \( 2\mu \)?
Applied Stochastic Processes - End of Sequence Case Study
Queue Crashers Mastery Key Instructor Solutions: Queue Crashers
Final Assessment Guide
Unit
Simulation & Analysis
01
Steady State Baseline
A. Traffic Intensity:
\( \lambda = 15 \) calls/hr. Service time = 3 mins \( \implies \mu = \frac{60}{3} = 20 \) calls/hr.
B. Average Wait Time (\( W \)):
\( W = \frac{1}{\mu - \lambda} = \frac{1}{20 - 15} = \frac{1}{5} \) hours = 12 minutes.
02
Stability Breakdown
# Qualitative Analysis
Since rho > 1, the system is unstable. The queue size N(t) is expected to grow
03
Multi-Server Expansion
Service Rate Function
\( \mu_n = \begin{cases} n\mu & 1 \le n < 2 \\ 2\mu & n \ge 2 \end{cases} \)
Simulated Insight
A single fast server (\( 2\mu \)) is generally more efficient at reducing wait times than two servers at rate \( \mu \), because a single server is never idle while a customer is waiting in the other's line. However, two servers provide better reliability.
Conclusion of Sequence: "Simulation and Analysis of Continuous-Time Models"
