# Section 2.3. Fixed point models. lambda0 <- 3 lambda1 <- 5 lambda2 <- 7 C1 <- 10 C2 <- 12 m <- 1 # A useful way of writing the Erlang function is # E(rho,C) = Prob(X=C)/Prob(X<=C) where X~Poisson(rho) # R has built-in functions for Prob(X=C) and Prob(X<=C) that we can use. erlang <- function(rho,C) dpois(C,rho)/ppois(C,rho) # Input: B, a vector consisting of two elements, B1 and B2 # Output: A vector of two elements, the new estimates for B1 and B2 updateB <- function(B) { B1old <- B[1] B2old <- B[2] B1new <- erlang((lambda1+lambda0*(1-B2old))*m,C1) B2new <- erlang((lambda2+lambda0*(1-B1new))*m,C2) c(B1new,B2new) } # I'll store estimates for (B1,B2) in a matrix. Row i contains the ith estimate. # Start off with some initial guess e.g. (.5,.5) # Repeatedly take the last guess (the bottom row of the matrix) and update it, # and attach this new estimate to the bottom of the matrix. bvals <- matrix(c(.5,.5),nrow=1) for (i in 1:10) { bold <- bvals[nrow(bvals),] bnew <- updateB(bold) bvals <- rbind(bvals,bnew) } # Look at the estimates we got. The estimates very quickly converge to # B1=0.1073, B2=0.1076 bvals # (ii) Plot a 2-dimensional drift diagram, B1 on x-axis, B2 on y-axis. drift <- function(B) updateB(B) - B # What is the drift, at different points? drift(c(.5,.5)) drift(c(0,1)) drift(c(1,0)) drift(c(1,1)) drift(c(0,0)) # Here is some more advanced R code, which plots many arrows showing the drift. # It also superimposes the successive estimates we got, i.e. bvals. # Notice that the arrows all point towards the fixed point, and that # these iterations very quickly took us there. library(lattice) # set up a data frame with all points at which we want to draw arrows drifts <- expand.grid(b1=seq(0,1,.05),b2=seq(0,1,.05)) # calculate the drifts at each of these points u <- t(sapply(1:nrow(drifts), function(i) updateB(c(drifts[i,'b1'],drifts[i,'b2'])))) drifts$b1new <- u[,1] drifts$b2new <- u[,2] # decide how long to draw each arrow (normalizing, so the arrows look OK) drifts$len <- ((drifts$b1new-drifts$b1)^2+(drifts$b2new-drifts$b2)^2)^0.2 drifts$d1 <- drifts$b1+(drifts$b1new-drifts$b1)*drifts$len/max(drifts$len)*0.1 drifts$d2 <- drifts$b2+(drifts$b2new-drifts$b2)*drifts$len/max(drifts$len)*0.1 # plot the iterations from the fixed point method, # and superimpose the drift arrows xyplot(bvals[,2]~bvals[,1], type='b', xlim=c(0,.55),ylim=c(0,.55), panel=function(...) { larrows(drifts$b1,drifts$b2,drifts$d1,drifts$d2, length=sqrt(drifts$len*.02),col='grey50') panel.xyplot(...,lwd=3,cex=3) }, xlab="b1",ylab="b2") # The same, but with baby steps eps <- .1 bvals <- matrix(c(.2,.5),nrow=1) for (i in 1:100) { bold <- bvals[nrow(bvals),] bnew <- updateB(bold) bvals <- rbind(bvals,eps*bnew+(1-eps)*bold) }