Skip to content

Instantly share code, notes, and snippets.

@willpearse
Created October 24, 2021 13:43
Show Gist options
  • Save willpearse/56b58c917f5aa19932ea46830232dc20 to your computer and use it in GitHub Desktop.
Save willpearse/56b58c917f5aa19932ea46830232dc20 to your computer and use it in GitHub Desktop.
Toy MCMC
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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment