Toy MCMC

toy-mcmc.R

data <- c(120, 130, 110, 105, 121)
n.iter <- 10000
mu.prior <- function(x) dnorm(x, mean=0, sd=10)
sd.prior <- function(x) dnorm(x, mean=0, sd=10)
posterior <- function(x, mu, sigma) sum(c(
                                         dnorm(x, mu, abs(sigma), log=TRUE),
                                         mu.prior(x), sd.prior(x)
                                     ))
proposal <- function(x) x + rnorm(1, 0, sd=.25)
mu <- 90
sigma <-5

output <- matrix(NA, nrow=n.iter+1, ncol=3)
colnames(output) <- c("mu", "sigma", "posterior")
output[1,] <- c(mu, sigma, posterior(data, mu, sigma))

for(i in 1:n.iter){
    new.mu <- proposal(mu)
    new.sigma <- proposal(sigma)
    post <- posterior(data, mu, sigma)
    new.post <- posterior(data, new.mu, new.sigma)
    if(runif(1) <= 10^(new.post-post)){
        mu <- new.mu
        sigma <- new.sigma
        post <- new.post
    }
    output[i+1,] <- c(mu, sigma, post)
}