#### Code by Matt Asher. Published at StatisticsBlog.com #### #### CONFIG #### # Where do you want the animation stored? setwd("path/to/directory/") # Number of slots to fill numbSlots = 5 # Total time intervals in the simulation intervals = 250 # Probability that a new person will show up during an interval # Note, a maximum of one new person can show up during an interval p = 0.1 # Average time each person takes at the teller, using discretized exponential distribution # Times will be augmented by one, so that everyone takes at least 1 interval to serve meanServiceTime = 60 #### INITIALIZATION #### queueLengths = rep(0, intervals) slots = list() waitTimes = c() leavingTimes = c() queue = list() arrivalTimes = c() frontOfLineWaits = c() id = 1 #### Libraries #### # Use the proto library to treat people like objects in traditional oop library("proto") #### Functions #### # R is missing a nice way to do ++, so we use this inc <- function(x) { eval.parent(substitute(x <- x + 1)) } # Main object, really a "proto" function # Changed based on suggestion by Johann (http://www.statisticsblog.com/2011/10/waiting-in-line-waiting-on-r/comment-page-1/#comment-12272) person <- proto( new = function (.,id) { # Set id and how muct teller time this person will need .\$proto(intervalsNeeded = floor(rexp(1, 1/meanServiceTime)) + 1,personId = id) }, intervalArrived = 0, intervalAttended = NULL, intervalsWaited = 0, intervalsWaitedAtHeadOfQueue = 0, ) # Fill our slots with Nullos (don't Google that term, NSFL) nullPerson = person\$new(0) for(k in 1:numbSlots) { slots[[k]] = nullPerson } #### Main loop #### for(i in 1:intervals) { # Check if anyone is leaving the slots for(j in 1:numbSlots) { if(slots[[j]]\$personId) { if( (slots[[j]]\$intervalAttended + slots[[j]]\$intervalsNeeded) == i) { # They are leaving the queue, slot to 0 slots[[j]] = nullPerson leavingTimes = c(leavingTimes, i) } } } # See if a new person is to be added to the queue if(runif(1) < p) { newPerson = person\$new(id) inc(id) newPerson\$intervalArrived = i queue = c(queue, newPerson) arrivalTimes = c(arrivalTimes, i) } # Can we place someone into a slot? for(j in 1:numbSlots) { # Is this slot free if(!slots[[j]]\$personId) { if(length(queue) > 0) { slots[[j]] = queue[[1]] slots[[j]]\$intervalAttended = i waitTimes = c(waitTimes, queue[[1]]\$intervalsWaited) # Only interested in these if person waited 1 or more intevals at front of line if(queue[[1]]\$intervalsWaitedAtHeadOfQueue) { frontOfLineWaits = c(frontOfLineWaits, queue[[1]]\$intervalsWaitedAtHeadOfQueue) } # Remove placed person from queue queue[[1]] = NULL } } } # Everyone left in the queue has now waited one more interval to be served if(length(queue)) { for(j in 1:length(queue)) { inc(queue[[j]]\$intervalsWaited) } # The (possibley new) person at the front of the queue has had to wait there one more interval. inc(queue[[1]]\$intervalsWaitedAtHeadOfQueue) } # End of the interval, what is the state of things queueLengths[i] = length(queue); # Create a snapshot of this interval. png(file=paste("snapshot-",paste(sprintf( "%04d", i),sep="",collapse=""),".png",sep=""), width=550, height=550) plot.new() plot(xlim=c(0, numbSlots), ylim=c(0, 20), 1, 1, pch=18, col=(slots[[1]]\$personId %% 153), cex=sqrt((slots[[1]]\$intervalsNeeded) ), main=paste("Interval ",i,sep=""), xlab="Queueing area", ylab="Slots") # Add the rest of the slots if(numbSlots > 1) { for(sl in 1:numbSlots) { points(sl,1, pch=18, col=(slots[[sl]]\$personId %% 153), cex=log(slots[[sl]]\$intervalsNeeded+1)) } } # Show the line, if there is one if(length(queue)) { for(j in 1:length(queue)) { points((numbSlots/2),(j+1), pch=18, col=(queue[[j]]\$personId %% 153), cex=log(queue[[j]]\$intervalsNeeded+1)) } } dev.off() } # Create an animated gif of all of this. You will need ImageMagick installed for this system('"path/to/convert" -delay 15 snapshot-*.png animation.gif') #### Output #### # par(ask=T) # So we can see plots one by one plot(queueLengths, type="o", col="blue", pch=20, main="Queue lengths over time", xlab="Interval", ylab="Queue length") #plot(waitTimes, type="o", col="blue", pch=20, main="Wait times", xlab="Person", ylab="Wait time") #if(length(frontOfLineWaits)) { hist(frontOfLineWaits, col="blue", breaks=10) } # Optional removal of all snapshot files. Careful with this. # file.remove(list.files(pattern=".png"))