## MCMC examples (Section 4.9) ## Example 1: # Vector to store values of Markov chain: tot.iter <- 10000 X <- rep(0, tot.iter) accept <- X X[1] <- 1 # setting an initial value for the chain for (n in 1:tot.iter){ Y <- rpois(1, lambda=X[n]) + 1 i <- X[n]; j <- Y; if(log(runif(1,0,1)) < 4*log(i)-j+(i-1)*log(j)+lfactorial(j-1) -(4*log(j)-i+(j-1)*log(i)+lfactorial(i-1)) ) { X[n+1] <- j; accept[n] <- 1 }else { X[n+1] <- i } } # end loop mean(accept) # a summary of our sample: table(X) ## Example 2: # Vector to store values of Markov chain: tot.iter <- 10000 X <- rep(0, tot.iter) accept <- X sig.sq <- 1 ### proposal distribution's variance (try some different values) X[1] <- 0 # setting an initial value for the chain for (n in 1:tot.iter){ Y <- rnorm(1, mean=X[n], sd=sqrt(sig.sq) ) i <- X[n]; j <- Y; if(log(runif(1,0,1)) < -(j + exp(-j)) - (-(i + exp(-i))) ) { X[n+1] <- j; accept[n] <- 1 }else { X[n+1] <- i } } # end loop mean(accept) # a picture of our sample, with summary statistics: plot(density(X)) mean(X) var(X) # Example 3: Gibbs sampler example: tot.iter <- 10000 X <- matrix(0,nrow=tot.iter,ncol=2) # matrix to store sampled random vectors: for (n in 1:(tot.iter-1)){ choice <- sample(1:2,size=1) other <- 3-choice X[n+1,choice] <- rnorm(1,mean=0.5*X[n,other], sd=sqrt(1-(0.5^2)) ) X[n+1,other] <- X[n,other] } X.stationary <- X[3001:10000,] # discarding the first 3000 values as burn-in # a picture of our sample, with summary statistics: plot(X.stationary[,1], X.stationary[,2]) apply(X.stationary,2,mean) # component means apply(X.stationary,2,sd) # component std. deviations cor(X.stationary[,1], X.stationary[,2]) # sample correlation