Understanding the Bayesian Grid Approximation Probabilities

Question

I am working through the survivability analysis shown here. My specific question is regarding the grid approximation method, implemented in the following R code:

# function to get log-likelihood of the data for a given shape & scale pair
grid_function <- function(shape, scale) {
  dweibull(data_tbl$fatigue_duration, shape = shape, scale = scale, log = T) %>%
    sum()
}

# set up grid of possible shape, scale parameters
n <- 100
shape_grid <- seq(1, 3, length.out = n)
scale_grid <- seq(60, 130, length.out = n)
two_param_grid <- expand_grid(shape_grid, scale_grid)

# map the grid_function over all candidate parameter pairs
# multiply LL by prior and convert to probability
full_tbl <- two_param_grid %>%
  mutate(log_likelihood = map2(shape_grid, scale_grid, grid_function)) %>%
  unnest() %>%
  mutate(
    shape_prior = 1,
    scale_prior = 1
  ) %>%
  mutate(product = log_likelihood + shape_prior + scale_prior) %>%
  mutate(probability = exp(product - max(product)))

I don't understand the last two lines creating the product and probability values in the table.

1. Why are we adding the shape/scale prior to the PDF value?
2. Why do we use product - max(product) in the probability calculation?

My understanding is that we use the log values to avoid datatype precision errors, but it's not clear to me why we don't just take the PDF log-likelihood values directly and instead use these operations. Does this approach change if we use different priors (e.g., non-flat priors)?

Answer 1

I don't understand the last two lines creating the product and probability values in the table.

1)Why are we adding the shape/scale prior to the PDF value?

You are adding them because the software is using the logarithmic transformation. You are really multiplying. It is not a pdf value, it is a likelihood.

With a pdf, the data is random. With a likelihood, the parameters are random.

An example of a pdf is $$\pi^{-1}[1+(x-\mu)^2]^{-1},\forall{x}\in\Re,$$ and example of a likelihood function would be $$\pi^{-1}[1+(x-\mu)^2]^{-1},\forall\mu\in\Re.$$

2)Why do we use product - max(product) in the probability calculation?

You are converting everything to a relative density. When $product=\max(product)$, then the relative density is 1 at the maximum value and all other values are expressed as a percentage of the maximum. That allows you to avoid having everything scaled at some problematic value such as 2.31 E -631. It helps reduce computational issues down the line.

My understanding is that we use the log values to avoid datatype precision errors, but it's not clear to me why we don't just take the PDF log-likelihood values directly and instead use these operations.

You will underflow with just a few values, or you will zero out. It can be challenging to work with sometimes even with logarithms. I was comparing two models and one model had relative odds of approximately $$10^{-65,000,000}:1.$$ Of course, that is zero, but the data set was enormous so the relative likelihoods were very small in both cases.

Does this approach change if we use different priors (e.g., non-flat priors)?

Yes. The prior would become a function and it may not be independent. You may have $p(\mu|\sigma)p(\sigma)$, for example. A prior is not a function of the data, so you would add the value of the log-prior based on the function that creates the prior.

Answer 2

Why do we use product - max(product) in the probability calculation?

I actually find this part to be a bit silly. The other answers here say it is to get a "relative probability" taken with respect to the maximally probable element, but that doesn't make sense to me. (Why would you want to know that?) A much better way to do this would be to scale things to get a proper probability distribution where the probabilities add up to one. To do this in log-space you would compute $\mathbf{l} - \text{logsumexp}(\mathbf{l})$ instead of $\mathbf{l} - \max(\mathbf{l})$. This can be done with the following alternate code for the last part:

# map the grid_function over all candidate parameter pairs
# multiply LL by prior and convert to probability
full_tbl <- two_param_grid %>%
  mutate(log_likelihood = map2(shape_grid, scale_grid, grid_function)) %>%
  unnest() %>%
  mutate(
    shape_prior = 1,
    scale_prior = 1
  ) %>%
  mutate(product = log_likelihood + shape_prior + scale_prior) %>%
  mutate(probability = exp(product - matrixStats::logSumExp(product)))

Answer 3

Here's my current understanding (hoping someone will correct it if it's off).

According to Bayes rule, the posterior distribution is proportional to the product of the prior and the likelihood.

$p(\theta|D) = p(D|\theta) p(\theta) / p(D)$

where D = data. I.e.,

$posterior = (likelihood)(prior)/(evidence)$

The priors defined above are flat, uninformed priors. In other words, we don't have information about the distribution of the parameters, so the prior beliefs are set to 100% probability:

shape_prior = 1,
scale_prior = 1

If we instead had information to adjust the strength of this belief, these priors would be adjusted to fit those new beliefs in accordance with how confident we are about the model assumptions. The likelihood probability is computed from the grid_function() in the link using the R dweibull() function to get the densities:

grid_function <- function(shape, scale) {
  dweibull(data_tbl$fatigue_duration, shape = shape, scale = scale, log = T) %>%
    sum()
}

So following Bayes rule, the posterior is product of the prior and likelihood. Since we're actually using the log-likelihood probabilities to avoid precision errors (log = T in the grid_function() above), this product becomes additive:

product = log_likelihood + shape_prior + scale_prior)

To ensure the conjoint probabilities (i.e. the probabilities across all possible parameter/data combinations) are relative with the max probability = 1, we divide all the products by the max product. Or, since we use the log-probabilities, we subtract them:

probability = exp(product - max(product)