#----------------------------------------------------------------
# Drift model for Erlang link

# First, load in the library for solving differential equations.
# If it's not already installed on your system, you need to
# pick "Package | Install package" from the menu and install it.
library(deSolve)
# and the library for plotting
library(lattice)

# Define the drift model.
# From state y, at time t, what is the drift in each component of y?
# In our case, there is only one component we're interested in, called 'n'.
driftmodel <- function(t,y,...,lambda,mu,C) {
  n <- y['n']
  drift <- list(dn=lambda-n*mu)
  drift
}
initstate <- c(n=0)

# Solve the differential equation
theorydata <- ode(initstate, times=seq(0,8,.1), driftmodel, lambda=150,mu=1,C=200)
theorydata <- as.data.frame(theorydata)

xyplot(n~time, data=theorydata, type='l')


# Add in a simulation
lambda <- 150; mu <- 1; C <- 200
trace <- list()
t <- 0; calls <- 0
while (t<8) {
  trace <- c(trace, list( c(t,calls) ))
  t <- t + rexp(1, lambda+mu*calls)
  calls <- calls + if(runif(1)<=lambda/(lambda+mu*calls)) 1 else -1 }
simdata <- data.frame(do.call('rbind',trace))
names(simdata) <- c('time','n')

df <- rbind(cbind(theorydata,type='driftmodel'), cbind(simdata,type='sim'))
xyplot(n~time, groups=type, data=df, type='l', auto.key=TRUE)