metropolis <- function(init, logp, M = 1000, stepsize = 1.0, ...) { D <- length(init) samples <- matrix(rep(NA, M * D), ncol = D) samples[1,] <- init ar <- 0 for (m in 2:M) { prev <- samples[m-1,] lp <- logp(prev, ...) prop <- rnorm(D, prev, stepsize) lp_prop <- logp(prop, ...) r <- lp_prop - lp a <- log(runif(1)) < min(r, 0) samples[m,] <- a * prop + (1 - a) * prev ar <- ar + (a - ar) / (m-1) } list(samples = samples, acceptance_rate = ar) } # example 1: 1D fake data ---- x <- rnorm(100, 3.14, 1) lp <- function(mu, data) { sum(dnorm(data, mu, 1, log=TRUE)) } init <- rexp(1) o <- metropolis(init, lp, stepsize = 0.1, data = x) mu_samples <- o$samples ar <- o$acceptance_rate plot(mu_samples, type = "l", main = "Trace plot") abline(h=3.14, lty = 2) hist(mu_samples, breaks = "FD", main = "Histogram") # example 2: 2D fake data ---- x <- rnorm(100, 3.14, sqrt(2)) lp <- function(theta, data) { mu <- theta[1] sd <- exp(theta[2]) sum(dnorm(data, mu, sd, log=TRUE)) } init <- rnorm(2) o <- metropolis(init, lp, stepsize = 0.1, M = 10000, data = x) samples <- o$samples tail(samples) plot(samples[, 1], type = "l", main = "Trace plot, mu") abline(h = 3.14, lty = 2) hist(samples[,1], breaks = "FD", main = "Histogram, mu") plot(exp(samples[, 2]), type = "l", main = "Trace plot, sigma") abline(h = sqrt(2), lty = 2) hist(exp(samples[,2]), breaks = "FD", main = "Histogram, sigma") # example 3: real data ---- df <- read.csv("https://roualdes.sfo3.digitaloceanspaces.com/data/abalone.csv") lp <- function(theta, data) { mu <- theta[1] sd <- exp(theta[2]) sum(dnorm(data, mu, sd, log=TRUE)) + dnorm(mu, log=TRUE) + dexp(sd, log=TRUE) } init <- rnorm(2) M <- 10000 o <- metropolis(init, lp, stepsize = 0.05, M = M, data = df$age) samples <- o$samples tail(samples) warmup <- 2000 wdx <- warmup:10000 plot(samples[wdx, 1], type = "l", main = "Trace plot, mu") hist(samples[wdx,1], breaks = "FD", main = "Histogram, mu") plot(exp(samples[wdx, 2]), type = "l", main = "Trace plot, sigma") hist(exp(samples[wdx,2]), breaks = "FD", main = "Histogram, sigma") # example 4: real data ---- lp <- function(theta, data) { beta0 <- theta[1] beta1 <- theta[2] sd <- exp(theta[3]) sum(dnorm(data$age, beta0 + beta1*df$shell_weight, sd, log=TRUE)) + dnorm(beta0, log=TRUE) + dnorm(beta1, log=TRUE) + dexp(sd, log=TRUE) } init <- rnorm(3) M <- 100000 o <- metropolis(init, lp, stepsize = 0.1, M = M, data = df) samples <- o$samples tail(samples) warmup <- 2000 wdx <- warmup:M plot(samples[wdx, 1], type = "l", main = "Trace plot, beta0") hist(samples[wdx,1], breaks = "FD", main = "Histogram, beta0") plot(samples[wdx, 2], type = "l", main = "Trace plot, beta1") hist(exp(samples[wdx,2]), breaks = "FD", main = "Histogram, beta1") plot(exp(samples[wdx, 3]), type = "l", main = "Trace plot, sigma") hist(exp(samples[wdx, 3]), breaks = "FD", main = "Histogram, sigma") plot(samples[wdx, 1], samples[wdx, 2], pch = 20, xlab="beta0", ylab="beta1")