Sampling from a grid-approximate posterior

The problem is a binomial problem, the Globe Tossing problem. Pick a point randomly in the globe and it will fall in water W or in land L. Supposse the data WLWWWLWLW (n=9, w=6), and the question: What is the most likely fraction of water on the surface of the globe? we don’t know, but we know the most likely ratio that accounts for the given data WLWWWLWLW, is 6/9 = 0.666.

The Likelihood

The same value can be obtained in a more complicated way as the maximum of the binomial distribution considered as a function of p, the probablity, given that 6 out of 9 samples were W. We can estimate this in a graph.

# define grid
p_grid <- seq(from=0, to=1, length.out = 1000)
# compute likelihood at each value in grid
likelihood <- dbinom(6, size = 9, prob = p_grid)
plot(p_grid, likelihood, type = 'l', xlab="probability of water")
mtext("1000 points")

You can appreciate the maximum occurs around 0.66. The maximum can be derived with calculus making the derivative with respect to p equal to zero.(DOIT). Also the Likelihood is not a distribution function, it does not add up to one for example.(DOIT).

The Posterior

A way to turn the likelihood into a probability density is by conditioning the likelihood to the prior. This is the posterior probability that we obtain using the Bayes formula.

# define an uninformative prior
prior <- rep(1, length(p_grid))
# compute the product of the likelihood and prior
unstd.posterior <- likelihood * prior
# standardize the posteriour, so it sums to 1
posterior <- unstd.posterior / sum(unstd.posterior)
plot(p_grid, posterior, type='l', xlab="probability of water")
mtext("1000 points")

The maximum likelihood now correstponds to the maximum aposteriory or MAP but only because the prior is uniform.

Sampling from a grid-approximated posterior

Now we wish to draw 10,000 samples from a probability distribution given as a vector of probability weights which is our posteriour.

samples <- sample(p_grid, size=1e4, replace=TRUE, prob=posterior)
plot(samples)

We now plot the density estimate computed from these samples using the convenient interface for plotting densities defined in the library “rethinking”

dens(samples)