library(lattice) #---------------------------------------------------------------- # Drift model for Dynamic Alternative Routing erlang <- function(rho,C) dpois(C,rho)/ppois(C,rho) drift <- function(b,rho,C) erlang(rho*(1+2*b*(1-b)),C) - b # A plot of the drift function B <- seq(0,.3,by=.01) Bdrift <- drift(B, rho=4750,C=5000) xyplot(Bdrift~B, type='l', panel=function(...) { panel.abline(h=0,col='grey60'); panel.xyplot(...) }) # Drift diagram rescale <- function(x) sign(x)*sqrt(abs(x)) xyplot(B~0:20, panel=function(x,...) { for (xx in x) larrows(xx,B,xx,B+rescale(Bdrift*.004),length=.2,unit='native', type='closed',col='grey70',fill='grey70') } ) # Finding the fixed points # I'll update my estimate for B 1000 times, and I'll do the process simulateneously # for a range of values of B. # Answer: we end up with two fixed points, 8.81e-6 and 2.11e-1. B <- seq(0,.3,by=.01) for (i in 1:1000) B <- B + drift(B, rho=4750,C=5000) B