A simple MC (example 1 of notes2)

MC1 <- function(n,a) {
 if( a >= 1 ) stop( "a must < 1" )
 if( a <= 0 ) stop( "a must > 0" )
 N <- as.integer(n)

vec <- rep(9, N)
vec[1] <- a
for(i in 2:N) vec[i] <- runif(1, min= 1-vec[i-1], max =1)
return(vec)
}













A simple Metropolis MCMC




Let's look at simulating from a normal density with zero mean and unit variance using a 
Metropolis algorithm with uniform as proposal distribution. 

Implementation in R
A function for the Metropolis sampler for this problem is given below. The chain 
is initialised at zero, and at each stage a U(-alpha, alpha) innovation is 
proposed. 

MetroNorm <- function (n, a) 
{
        vec <- vector("numeric", n)
        x <- 0
        vec[1] <- x
        for (i in 2:n) {
                can <- x + runif(1, min= -a, max=a)
                aprob <- dnorm(can)/dnorm(x)
                u <- runif(1)
                if (u < aprob) 
                        x <- can
                vec[i] <- x
        }
        return(vec)
}


Bayesian Example:  Casella & Berger example 7.2.16

 xsample <-  rnorm(10, mean=3, sd= 2)

We assume mean above is unknown parameter, needs to be estimated. While sd=2 is known.

We took The prior distribution of the parameter
normal ( mean = 2, sd= 4 )

fprior <- function(mu) { dnorm(mu, mean=2, sd=4) }

The likelihood function

lik <- function(mu) { prod( dnorm(xsample, mean=mu, sd=2) )  }

The posterior is then proportional to (need xsample to be available)

 
fpost <- function(mu){ dnorm(mu, mean=2, sd=4) * prod( dnorm(xsample, mean=mu, sd=2) ) }

Use this fpost( ) to replace dnorm( ) in the MetroNorm( ) function.