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?
devtools::install_github
. $\endgroup$range(post$prod)
we see that all values are below -1230 whileexp(-1230)
is numerically zero. $\endgroup$1.716 * max(post$prod)
. $\endgroup$