To formulate the problem, consider calculating the posterior probability given a series of weighted likelihoods. Let
\[ x_i \]
be the likelihood for observation i. We wish to calculate the posterior probability as:
\[ p_i = \frac{x_i}{\sum_j x_j} \]
As Robert Eisele very nicely discusses in his blog post, we can run into underflow issues. Suppose you have the following vector of logged likelihoods.
x <- c(-9000, -3000, -1000)
exp(x)/sum(exp(x))
# [1] NaN NaN NaNAs we can see, all the values will be zero. However, consider making the calculation in log-space.
\[\begin{align} \text{log}\left(\frac{x_i}{\sum_j x_j}\right) &= \text{log}x_i - \text{log}\sum_j x_j \\ &= \text{log}x_i - \text{log}\sum_j e^{\text{log}x_j} \\ \end{align}\]
Given the proof from Robert Eisele, we know:
\[ \text{log}\sum_i e^{y_i} = \text{log}\sum_i e^{y_i - a} + a \]
for arbitrary a. Let a be the maximum log-likelihood value, and we get:
\[\begin{align} \text{log}\left(\frac{x_i}{\sum_j x_j}\right) &= \text{log}x_i - \left(\text{log}\sum_j e^{\text{log}x_j - \operatorname{max}_j(\text{log}x_j)} + \operatorname{max}_j(\text{log}x_j) \right) \\ p_i & = \text{exp}\left\{ \text{log}x_i - \left(\text{log}\sum_j e^{\text{log}x_j - \operatorname{max}_j(\text{log}x_j)} + \operatorname{max}_j(\text{log}x_j) \right) \right\} \end{align}\]
Finally, a simple R function to do the work:
getPost <- function(l) {
## l: vector of logged-likelihood values
exp(l - (log(sum(exp(l - max(l)))) + max(l)))
}
getPost(x)
# [1] 0 0 1and we see it gives the exact same values for tractable likelihoods:
y <- c(-9, -3, -1)
exp(y)/sum(exp(y))
# [1] 0.0002953872 0.1191677110 0.8805369018
getPost(y)
# [1] 0.0002953872 0.1191677110 0.8805369018