# SIMULATE A QUEUE. One person at the front of the queue is served at each # (integer) time, after which some number of new people arrive, with the number # being Poisson distributed with the given mean. The queue can hold # a maximum of max people; any beyond that have to go away (and never return). # The queue is simulated (starting with an empty queue) for the given number # of steps. The result returned is a list with an element called "length" # that is a vector of size steps containing the length of the queue after # each step, and a element called "away" that is a vector of size steps # containing the number of people who went away because the queue was full # during each time step. simulate.queue = function (steps, mean, max) { # Allocate space for vectors that will contain results. length = rep(0,steps) away = rep(0,steps) # Set the length of the queue to zero to start. qlen = 0 # Loop for the requested number of steps. for (t in 1:steps) { # Take one person out of the queue, if there is one. if (qlen>0) { qlen = qlen - 1 } # Increase the number in the queue by the number of new arrivals. qlen = qlen + rpois(1,mean) # If there are too many people in the queue, force some to go away. if (qlen>max) { away[t] = qlen-max qlen = max } # Record how long the queue is after this step. length[t] = qlen } # Return the list with the results. list (length=length, away=away) } # FIND THE DISTRIBUTION OF QUEUE LENGTH USING MATRIX OPERATIONS. Assumes # that the queue starts out empty. Returns a vector of state probabilities # for the distribution after "steps" steps. queue.length.distribution = function (steps, mean, max) { # Create the matrix of one-step transition probabilities. P.1 = matrix(0,max+1,max+1) for (j in 0:max) { P.1[1,j+1] = dpois(j,mean) } P.1[1,max+1] = P.1[1,max+1] + (1-ppois(max,mean)) for (i in 1:max) { for (j in (i-1):max) { P.1[i+1,j+1] = dpois(j-(i-1),mean) } P.1[i+1,max+1] = P.1[i+1,max+1] + (1-ppois(max-(i-1),mean)) } # Find the "steps" power of the one-step transition matrix. P.steps = diag(max+1) # Start with (max+1) by (max+1) identity matrix for (t in 1:steps) { P.steps = P.steps %*% P.1 } # Return the first row of P.steps, giving the distribution after # "steps" transitions starting with an empty queue. P.steps[1,] }