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

Question

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!

Answer 1

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) } }