# Understanding bayesian model code from chapter 4 of “Statistical rethinking”

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?

• Might that simply be normalizing to prevent overflow/underflow when computing exp? – Glen_b Oct 12 '18 at 8:31
• Since rethinking is not a CRAN package, you could add a call to devtools::install_github. – Yves Oct 12 '18 at 8:35
• @Glen_b I don't think that this would be the issue here. The formula is quite clear on what we are calculating. Where did the integral (corresponding to $Pr(data)$) go? – January Oct 12 '18 at 8:59
• Glen_b gave the answer. Looking at ̀range(post$prod) we see that all values are below -1230 while exp(-1230) is numerically zero. – Yves Oct 12 '18 at 14:46 • Um, but that means that we might as well subtract 1.716 * max(post$prod). – January Oct 12 '18 at 19:32

## 2 Answers

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.

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.

• I am aware that you do not always need the marginal likelihood, e.g. when calculating the Bayes factor. However, the book explicitly states that we are calculating the posterior distribution. If you use an arbitrary value in the denominator for purely technical reasons, maybe you are calculating something useful, but that thing is not the posterior distribution. Am I wrong? Thus, I edited my question. – January Oct 12 '18 at 19:28
• If the posterior is required, it might be that it is normalized later. In many cases, if the posterior is of simple form, the exact normalizing constant is readily calculated. – Glen_b Oct 13 '18 at 1:02