Jim — Sep 22, 2013, 10:00 PM
# Brute-force method of doing Bayesian inference
# for a single parameter
# Suppose y is binomial(p)
# The proportion is given a normal(mu, sigma) prior
# Want to simulate from posterior of p
# Set up a grid of p values from .01 to .99 in steps of 0.02
p <- seq(0.01, 0.99, 0.02)
# Suppose mu = 0.3, sigma = 0.1
# We take a sample of 15 and observe 9 successes
mu <- 0.3; sigma <- 0.1
n <- 15; y <- 9
# Compute the product of the likelihood and prior on the grid
product <- dbinom(y, size=n, prob=p) * dnorm(p, mean=mu, sd=sigma)
# Normalize the products
posterior <- product / sum(product)
# Plot posterior
plot(p, posterior, type="l")
# Simulate from posterior
sim.p <- sample(p, size=1000, replace=TRUE, prob=posterior)
# Construct density estimate from simulated values
plot(density(sim.p))
