Weighted Bootstrapping in obtaining posterior distribution

[GaryBAYLOR]

An example for binomial (n, p) likelihood and beta (a, b) prior.

wb.bin <- function(n, x, a = 1, b = 1, simulation = 5e3, k = 1e3) {
    random.prior <- rbeta(simulation, a, b)
    likelihood <- dbinom(x, n, random.prior)
    random.posterior <- sample(random.prior, k, replace = TRUE, prob = likelihood)
    actual <- (a + x)/(a + b + n)
    simulated <- mean(random.posterior)
    return(list(actual = actual, simulated = simulated))
}

When running this function we get

> wb.bin(20, 3)
$actual
[1] 0.1818182

$simulated
[1] 0.1844998

> wb.bin(12, 4)
$actual
[1] 0.3571429

$simulated
[1] 0.3570672

This method is pretty easy to understand and to implement. Thus it is great!