# Define real pars mu and sigma, sample 100x | |
trueMu <- 5 | |
trueSig <- 2 | |
set.seed(100) | |
randomSample <- rnorm(100, trueMu, trueSig) | |
# Grid approximation, mu %in% [0, 10] and sigma %in% [1, 3] | |
grid <- expand.grid(mu = seq(0, 10, length.out = 200), | |
sigma = seq(1, 3, length.out = 200)) | |
# Compute likelihood | |
lik <- sapply(1:nrow(grid), function(x){ | |
sum(dnorm(x = randomSample, mean = grid$mu[x], | |
sd = grid$sigma[x], log = T)) | |
}) | |
# Multiply (sum logs) likelihood and priors | |
prod <- lik + dnorm(grid$mu, mean = 0, sd = 5, log = T) + | |
dexp(grid$sigma, 1, log = T) | |
# Standardize the lik x prior products to sum up to 1, recover unit | |
prob <- exp(prod - max(prod)) | |
# Sample from posterior dist of mu and sigma, plot | |
postSample <- sample(1:nrow(grid), size = 1e3, prob = prob) | |
plot(grid$mu[postSample], grid$sigma[postSample], | |
xlab = "Mu", ylab = "Sigma", pch = 16, col = rgb(0,0,0,.2)) | |
abline(v = trueMu, h = trueSig, col = "red", lty = 2) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment