## 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