Understanding bayesian model code from chapter 4 of "Statistical rethinking"

Question

I'm trying to learn bayesian statistics from "Statistical rethinking" by Richard McElreath. In chapter 4, a model with Gaussian distribution of heights is introduced:

$h_i \sim N(\mu, \sigma)$

$\mu \sim N(178, 20)$

$\sigma \sim U(0, 50)$

To calculate the posterior probability, we need to calculate this:

$Pr(\mu, \sigma|h) = \frac{Pr(h|\mu, \sigma) Pr(\mu, \sigma)}{Pr(h)} = \frac{\prod_i N(h_i|\mu, \sigma) N(\mu|178,20) U(\sigma|0, 50)}{\int\int \prod_i N(h_i|\mu, \sigma) N(\mu|178,20) U(\sigma|0, 50)d\mu d\sigma}$

So far, so good. McElreath gives the following code to calculate the grid approximation of the posterior distribution:

## load the library and prepare the data
library(devtools)
devtools::install_github("rmcelreath/rethinking")
library(rethinking)
data(Howell1)
d2 <- Howell1[ Howell1$age >= 18, ]

mu.list <- seq(140, 160, length.out=200)
sigma.list <- seq(4, 9, length.out=200)
post <- expand.grid(mu=mu.list, sigma=sigma.list)
post$LL <- sapply(1:nrow(post), function(i) sum(dnorm(d2$height, 
             mean=post$mu[i], sd=post$sigma[i], log=TRUE)))

OK, this is clear enough, that is the $\prod_i N(h_i|\mu, \sigma)$ part, logarithmized so we calculate the sum rather than product.

post$prod <- post$LL + dnorm(post$mu, 178, 20, TRUE) + dunif(post$sigma, 0, 500, T)

Fine, we now multiply the likelihood by the prior $N(\mu|178,20) U(\sigma|0, 50)$

post$prob <- exp(post$prod - max(post$prod))

Aaaand I'm lost. Instead of calculating the integral

$\int\int \prod_i N(h_i|\mu, \sigma) N(\mu|178,20) U(\sigma|0, 50)d\mu d\sigma$

we now divide by the maximum. Why? I can't wrap my head around that. I am not sure how it should look like (should we sum the products?) but it certainly does not look like the maximum to me.

EDIT. Assuming that it is just a hack to avoid over/underflow in computations (as suggested in comments), how should the correct calculation (i.e. calculating precisely the marginal probability) look like?

Answer 1

I think I have the answer having played around with it for a while. Firstly, it is important to remember that the denominator is just a scaling constant and really all it does is re-scale the area under the posterior to 1.

McElreath doesn't bother calculating this marginal likelihood and instead re-scales. However, he doesn't seem to scale properly and instead re-scales using the max value. The reason he uses this max value is due to numerical limitations. Calculate anything above exp(-745)=4.94e-324 and the answer in R is 0. As the values in post$prod go from -3822 to -1228 if they were just exponentiated they would all go to 0! Instead in log-space he subtracts the max value to rescale them all (now from -2594 to 0) so all values in the main mass of probability will have values greater than 0 when exponentiated. However, he doesn't complete the final step as sum(post$prob)=301.83 which is incorrect as it should sum to 1. Instead to be formally correct he should have added this final line:

post$prob <- post$prob/sum(post$prob)

It makes little practical difference if it doesn't sum to 1 as for plotting it is the relative values in the posterior distribution that matter and not the absolute values when viewing the shape of the posterior. In fairness, he does say in the book that he won't bother going into the technical details of the code but it is a little misleading.

So, that final line is only a re-scaling due numerical constraints and not a calculation of the marginal likelihood. It should really have a proper scaling to 1 when back on the linear scale.

Answer 2

The denominator in the expression for the posterior is called the marginal likelihood or the evidence. One way to think about it is as a normalizing constant that ensures the posterior integrates to 1. This means the numerator is proportional to the posterior. In some applications (e.g. inference using MAP estimation or some MCMC methods) that's all we need, and the normalizing constant can be ignored. This is convenient because the integral is often inconvenient (or outright intractable) to compute.

I think this is also what's happening in your example. The numerator is calculated at a grid of points, then re-scaled so that the maximum value is 1. The goal is probably just to inspect/plot values that are proportional to the posterior.

The log of the numerator is calculated at a grid of parameter values. The max is then subtracted before exponentiating. The reason to calculate things this way is to avoid numerical overflow/underflow, as Glen_b mentioned. As above, this is probably just for plotting purposes--for inference, we can often work directly with the log and no re-scaling is needed.