# Metropolis sampling for the model: # # x1, x2, ..., xn ~ N(mu,1/tau) # # mu ~ N(0,1) # tau ~ Exp(1) # # The function below takes the data as its first argument, followed by initial # values for mu and tau, the standard deviations for the normal proposal # distributions, and the number of iterations to do. It returns a list # containing vectors of values for mu and tau at each iteration. # # The function could be speeded up by doing things like pre-computing some # sufficient statistics but this isn't done here for simplicity. met <- function(x,mu0,tau0,s.mu,s.tau,k) { mu.vec <- numeric(k) tau.vec <- numeric(k) mu <- mu0 tau <- tau0 for (i in 1:k) { prop.mu <- rnorm(1,mu,s.mu) prop.tau <- rnorm(1,tau,s.tau) if (prop.tau>0) { current.log.posterior <- log.posterior(x,mu,tau) prop.log.posterior <- log.posterior(x,prop.mu,prop.tau) if (runif(1) < exp(prop.log.posterior-current.log.posterior)) { mu <- prop.mu tau <- prop.tau } } mu.vec[i] <- mu tau.vec[i] <- tau } list (mu=mu.vec, tau=tau.vec) } # LOG POSTERIOR DENSITY. Returns the log of the posterior probability # density of mu and tau, given data x, plus some arbitrary constant. log.posterior <- function(x,mu,tau) { log.prior <- (-tau) + (-mu^2/2) log.likelihood <- log(tau)*length(x)/2 + sum (-tau*(x-mu)^2/2) log.prior + log.likelihood }