# M/M/1 queue simulator 
# experiment with different durations for simulation
lambda <- 2 # arrival rate
mu <- 2.5 # service rate
duration <- 1000 # total duration of the simulation
t <- 0 # current time
queue <- 0 # start with empty queue
s <- 0 # running sum for computing average queue length

# first arrival to start process
T1 <-rexp(1,rate=lambda)
currentqueue <- 1
eventsTime <- T1
t <- T1
nEvents <- 1 # total number of events that have occurred

while (t<duration) {
  nEvents <- nEvents+1
  if(currentqueue>0) {
    T1 <-rexp(1,rate=lambda+mu) # time until next event
    # is event an arrival or departure?
    p <- runif(1,0,1) 
    queue[nEvents]<-currentqueue # how many have been in queue before this new event
    currentqueue<-ifelse(p<lambda/(lambda+mu),
                         currentqueue+1, # arrival
                         currentqueue-1) # departure
  } else {
    T1 <-rexp(1,rate=lambda)
    queue[nEvents] <- currentqueue
    currentqueue <- 1
  }
  t <- t+T1 # time of next arrival
  eventsTime[nEvents] <- T1 # inter-event time
  s <- s+T1*queue[nEvents] # spent T1 at previous queue value
}

plot(cumsum(eventsTime),queue,type="s", xlab="Time",ylab="Queue length",
     main=paste("M/M/1 Simulation, mean queue length=",s/t))

# compare expected values to observed
# expected value of queue length
lambda/(mu-lambda) # theoretical 
s/t # simulation value

# compare to equilibrium distribution
pibar <-(1-lambda/mu)*(lambda/mu)^(1:20)
nubar <- matrix(0,1,20) 
for (k in 1:20) nubar[k]<-sum(eventsTime[which(queue==k)])/t
plot(1:20,pibar,type='l',xlab="Queue length",ylab="Proportion of time")
points(1:20,nubar)

# L2 norm of difference (note not finite state space here)
sqrt(sum((pibar-nubar)^2))