#
# Simulating departure times from an M(lambda)/U(beta1,beta2)/1 queue
#
dsim <- function(n,lambda,beta1,beta2){
  u <- rexp(n,lambda)
  a <- cumsum(u)
  v <- runif(n,beta1,beta2)
  d <- 0
  for (i in 1:n){
    d<- c(d,v[i]+max(a[i],d[length(d)]))
  }
  d <- d[2:(n+1)]
  return(d)
}
#
# Generate an ABC sample from the posterior of beta2, assuming lambda and beta1 are known and a U(beta1,beta2max) prior for
# beta2.
#
dabc <- function(d,lambda,beta1,beta2max,iters=1000,errprop=0.01){
  niters <- ceiling(iters/errprop)
  beta2 <- runif(niters,beta1,beta2max)
  ddist <- rep(0,niters)
  for (i in 1:niters){
    simd <- dsim(length(d),lambda,beta1,beta2[i])
    ddist[i] <- sum((d-simd)^2)
  }
  e <- sort(ddist,index.return=TRUE)
  err <- e$x
  e <- e$ix
  beta2 <- beta2[e]
  err <- err[iters]
  beta2 <- beta2[1:iters]
  return(list("beta2"=beta2,"err"=err))
}

#
# Example run
#
# Generates 5 departure times from M(1)/U(0,3).
#
d <- dsim(5,1,0,3)
#
# Simulates 1000 values from the approximate posterior for the model M(1)/U(0,beta2) with a U(0,5) prior for beta2. 
#
c <- dabc(d,1,0,5,1000,0.001)
beta2 <- c$beta2
err <- c$err
hist(beta2,probability=TRUE)