# Numerically stable transformation of log-likelihoods to probability [closed]

I have the following problem: I have log-likelihoods that need to be transformed to probabilities. One thing I have attempted is the following. Define $$\kappa_{s} := \log \int p (\theta | X_s) \text{d} \theta$$ and $$\nu_{j,s} := \log p (\theta_j | X_s$$ for data-generating processes $$X_s$$ and parameters $$\theta_j$$.

The log-likelihood $$\kappa_{s} \in [-\infty,\infty]$$. I have found $$\frac{\exp( \kappa_s - b )}{\sum_{s} \exp( \kappa_s - b)}$$ for $$b= \underset{s}{ \max} \kappa_{s}$$.

However, my log-likelihoods are still fairly spread out so that I wind up with lots of probabilities that are very close to zero with the maximum occupying a one. Is there a type of outlier-robust transformation or something slightly more sophisticated than the scheme described above?

For more context, the problem is more carefully described by

$$\sum_{j} \exp(\nu_{s,j}) \lambda_{j} \geq \exp (\kappa_{s})$$

which appears as a constraint of a linear programme. The log-likelihoods are fairly large and dispersed.

The reason why I want to transform them is that the LHS of the inequality approximates the integral $$\sum_{j} \exp(\nu_{s,j}) \lambda_{j} \approx \int p_{\nu} \, \text{d} \lambda$$.

The $$\text{d} \lambda$$ weights are the variable I am solving for, so logsumexp sadly does not seem to be applicable here.

This is a likelihood ratio test constraint.

Thank you!

• In what sense are those quantities "log-likelihoods" and why do you need to transform them to probabilities ? Commented May 5, 2023 at 18:29
• @J.Delaney thank you for asking! I added more context. Commented May 5, 2023 at 19:21
• I am afraid your edit made the question even less comprehensible. You introduced a bunch of symbols without explaining anything about what they represent or what you want to do with them. We cannot guess what any of this is supposed to mean. Commented May 5, 2023 at 19:38
• @J.Delaney thank you for your feedback, I added more context? The densities are fairly difficult to write in closed form but they give rise to misbehaving log-likelihoods. Commented May 5, 2023 at 19:47
– whuber
Commented May 5, 2023 at 19:58

If I understand your problem correctly (which I may not have), you want to compute the following quantity in a numerically stable way:

$$F_s(\boldsymbol{\kappa}, b) \equiv \frac{\exp( \kappa_s - b )}{\sum_{s} \exp( \kappa_s - b)}.$$

To do this you need to do computation in log-space, which requires you to manipulate your function into a form that does not involve any intermediate computation outside of log-space. For the type of function you are dealing with (which is a form of the softmax function) we can write:

\begin{align} \log F_s(\boldsymbol{\kappa}, b) &= \log \bigg(\frac{\exp( \kappa_s - b )}{\sum_{s} \exp( \kappa_s - b)} \bigg) \\[6pt] &= \log \Big( \exp( \kappa_s - b ) \Big) - \log \Big( \sum_{s} \exp( \kappa_s - b) \Big) \\[8pt] &= \kappa_s - b - \text{logsumexp} (\boldsymbol{\kappa} - b). \\[6pt] \end{align}

(In the last line we take $$\boldsymbol{\kappa} - b$$ to refer to the vector of values $$\kappa_s - b$$ taken over all indices $$s$$.) This form puts things in terms of the logsumexp function, which can itself be computed in log-space in a numerically stable way (see this related question). Consequently, we can compute the original function of interest in log-space (converting back to regular space at the end) as:

$$F_s(b, \boldsymbol{\kappa}) = \exp(\kappa_s - b - \text{logsumexp} (\boldsymbol{\kappa} - b)).$$

Coding the function: The softmax function is already coded in all relevant mathematical and statistical software, so you should be able to find it if you look around in the software you are using. If you want to code it from scratch in R you can do it like this:

#Create softmax function
softmax <- function(x, b = 0, log = FALSE) {

#Check inputs
if (!is.vector(x))          stop('Error: Input x should be a vector')
if (!is.numeric(x))         stop('Error: Input x should be a numeric vector')
if (!is.vector(b))          stop('Error: Input b should be a vector')
if (!is.numeric(b))         stop('Error: Input b should be a numeric vector')
if (length(b) != 1)         stop('Error: Input b should be a single numeric value')
if (!is.vector(log))        stop('Error: Input log should be a vector')
if (!is.logical(log))       stop('Error: Input log should be a logical vector')
if (length(log) != 1)       stop('Error: Input log should be a single logical value')

LOGS <- (x-b) - matrixStats::logSumExp(x-b)
if (log) { LOGS } else { exp(LOGS) } }

• Now my probabilities are much less degenerate thank you so much Ben The linked questions helped too but my question is sufficiently different. @whuber et al. can you please reopen it? Ben's answer is on point and solved my dilemma but I don't want the closed status to deter future readers. Ben's answer teaches facility with computations in log space which I was not used to. (To theorists, log space computations seem redundant and are thus less intuitive.) Commented May 6, 2023 at 11:17
• I didn't vote to close this question, Keynes, pending some clarification, but it looks like a duplicate of others about computing with logarithms. You can find even more threads like it by searching for logsumexp.
– whuber
Commented May 6, 2023 at 14:37
• In fairness to the moderators who closed your question, it certainly wasn't clear to me if I was interpreting it correctly, and the clarity of the question is below what we would usually expect. I recommend you edit to strip your question back to essentials (i.e., wanting to compute softmax function with a stable method). A lot of the material where you mention your likelihood is confusing.
– Ben
Commented May 6, 2023 at 23:10