# # Written by: # -- # John L. Weatherwax 2009-04-21 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # #----- library(lattice) library(zoo) source('chap_3_prob_12.R') # read in the definition of the function "Markov" set.seed(0) Pb = matrix( c(0.6,0.2,0.2,0.2,0.4,0.4,0.4,0.3,0.3), nrow=3, ncol=3, byrow=T ) # Part (a): # N = 1000 print( sprintf("Simulating the given Markov chain for %d timesteps",N) ) chain = Markov( N=N, initial.value=2, Pb ) # initial.value in {0,1,2} <=> {Sun,Cloud,Rain} print( table( chain ) / N ) # Part (b): # plotmarkov = function(n=10000, start=0, window=100, transition=Pb, npanels=5 ){ xc2 = Markov(n, start, transition) mav0 = rollmean(as.integer(xc2==0), window) mav1 = rollmean(as.integer(xc2==1), window) npanel = cut(1:length(mav0), breaks=seq(from=1, to=length(mav0), length=npanels+1), include.lowest=TRUE) df = data.frame( av0=mav0, av1=mav1, x=1:length(mav0), gp=npanel ) print( xyplot( av0+av1 ~ x | gp, data=df, layout=c(1,npanels), type="l", par.strip.text=list(cex=0.65), scales=list(x=list(relation="free")) ) ) } plotmarkov()