Nonstationary Poisson ProcessesDiscrete-Event Simulation:A First CourseSection 7.5: Nonstationary Poisson ProcessesSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesSection 7.5: Nonstationary Poisson ProcessesSuppose we want the arrival rate λ to change over time: λ(t)Recall the algorithm to generate a stationary Poisson process:Stationary Poisson Processa0= 0.0;n = 0;while (an< τ) {an+1= an+ Exponential(1 / λ);n++;}return a1, a2, ..., an−1;Above algorithm generates a stationary Poisson processTime interval is 0 ≤ t < τEvent times are a1, a2, a3, . . .Process has constant rate λSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesIncorrect AlgorithmChange constant λ to function:Incorrect Algorithma0= 0.0;n = 0;while (an< τ) {an+1= an+ Exponential(1 / λ(an));n++;}return a1, a2, ..., an−1;Incorrect: ignores future evolution of λ(t) after t = anIf λ(an) < λ(an+1) then an+1− anwill tend to be toolargeIf λ(an) > λ(an+1) then an+1− anwill tend to be toosmall“Inertia error” will be small only if λ(t) varies slowly with tSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesExample 7.5.1Piecewise-constant rate functionλ(t) =1 0 ≤ t < 152 15 ≤ t < 351 35 ≤ t < 50Simulated using incorrect algorithm for τ = 50Process replicated 10000 timesPartitioned interval 0 ≤ t ≤ 50 into 50 binsCounted number of events for each bin, divided by 10000Result isˆλ(t), an estimate of λ(t)Section 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesIncorrect process generation0 15 35 500.51.01.52.02.5ˆλ(t)tλ(t) · · ·ˆλ(t) —Incorrect algorithm has inertia error:ˆλ(t) under-estimates λ(t) after the rate increaseˆλ(t) over-estimates λ(t) after the rate decreaseUse one of 2 correct algorithms insteadSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesThinning methodDue to Lewis and Shedler, 1979Uses an upper bound λmax≥ λ(t) for 0 ≤ t < τGenerates a stationary Poisson process with rate λmaxDiscards (thins) some events, probabilisticallyEvent at time s is kept with probability λ(s)/λmaxans0λmaxλ(t)t.......................................................................................................................................................................................................................................................................................................................................................................................................................................................•λ(s)Efficiency depends on λmaxbeing a tight boundSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesAlgorithm 7.5.1Algorithm 7.5.1a0= 0.0;n = 0;while (an< τ) {s = an;do {s = s + Exponential(1 / λmax);u = Uniform(0, λmax);} while (u > λ(s));an+1= s;n++;}return a1, a2, ..., an−1;When λ(s) is low,The event at time s is more likely to be discardedThe number of loop iterations is more likely to be largeSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesExample 7.5.2The thinning method was applied to Example 7.5.1, usingλmax= 2Computation time increased by a factor of about 2.2The algorithm is not synchronizedEven if a separate stream is used for Uniform0 15 35 500.51.01.52.02.5ˆλ(t)λ(t) · · ·ˆλ(t) —tSection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesInversion MethodDue to C¸inlar, 1975Similar to inversion for random variate generationRequires only one call to Random per eventBased upon the cumulative event rate function:Λ(t) =Zt0λ(s) ds 0 ≤ t < τΛ(t) represents the expected number of events in interval[0, t)If λ(t) > 0 thenΛ(·) is strictly monotone increasingThere exists an inverse Λ−1(·)Section 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesAlgorithm 7.5.2: ideaGenerates a stationary “unit” Poisson process u1, u2, u3, . . .Equivalent to n random points in interval 0 < ui< Λ(τ)Each uiis transformed into aiusing Λ−1(·)a1a2a3a40u1u2u3u4Λ(t)t.........................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................Section 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesAlgorithm 7.5.2: detailsAlgorithm 7.5.2a0= 0.0;u0= 0.0;n = 0;while (an< τ) {un+1= un+ Exponential(1.0);an+1= Λ−1(un+1);n++;} return a1, a2, ..., an−1;The algorithm is synchronizedUseful when Λ−1(·) can be evaluated efficientlySection 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesExample 7.5.3Use the rate function from Example 7.5.2λ(t) =1 0 ≤ t < 152 15 ≤ t < 351 35 ≤ t < 50By integration we obtainΛ(t) =t 0 ≤ t < 152t − 15 15 ≤ t < 35t + 20 35 ≤ t < 50Solving u = Λ(t) for t we obtainΛ−1(u) =u 0 ≤ u < 15(u + 15)/2 15 ≤ u < 35u − 20 35 ≤ u < 50Section 7.5: Nonstationary Poisson Processes Discrete-Event Simulationc2006 Pearson Ed., Inc. 0-13-142917-5Nonstationary Poisson ProcessesExample 7.5.3 results0 15 35 500.51.01.52.02.5ˆλ(t)λ(t) · · ·ˆλ(t) —tGeneration time using inversion is similar to the incorrectalgorithmFor this example, Λ(t) was easily invertedIf Λ(t) cannot be inverted in closed formUse thinning if λmaxcan be foundUse numerical methods to invert