# Understanding the Bayesian Grid Approximation Probabilities

I am working through the survivability analysis shown here:

https://www.r-bloggers.com/2020/01/survival-analysis-fitting-weibull-models-for-improving-device-reliability-in-r/

My specific question is regarding the grid approximation method:

# 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)? ## 2 Answers 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. • Thank you for explaining. Can the prior be described in terms of discrete probabilities, such as those generated in a Markov chain? Feb 16 at 20:37 • If you cut a space into a grid, and if the grid covers the space where everything out of that space is below the computational precision of your computer, then the grid is a close, discrete approximation of the posterior. The grid has to cover the area where any addition would impact the computation. Computers do not let you add$10^{-1}+10^{-100}$. The answer would just be 0.1. The additional value would be effectively zero. The grid has to cover everywhere that is not effectively zero. Feb 16 at 20:41 • It is probably not good to think in terms of Markov chains, even though they are used in other types of numerical integration. The posterior would be converted to a biased, discrete approximation of the posterior, but if you change the boundary rules, you change everything potentially. Markov Chain Monte Carlo is a different method to do the same thing. It is not biased. For most purposes, though, the grid is more than adequate. Grids do not work well as dimensions increase. Feb 16 at 20:46 • Sorry, I may be misunderstanding. I get that the grid discretizes the posterior. What I'm trying to understand is whether the prior must be a continuous function, or if it can likewise be discretized? Feb 16 at 20:46 • It can also be discretized. Indeed, it should be if the likelihood is discretized. Feb 16 at 20:47 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)