Sample multinomial distribution from the log-probabilities

Question

I have three numbers:

lp_1 <- -316480
lp_2 <- -85041
lp_3 <- -1197042

The goal is to sample a multinomial distribution with probabilities proportional to the exp(lp_i). With these numbers, one gets 0 if one takes the exponential, and even if we set one number to 0, we will not manage to get some probabilities proportional to exp(lp_i) because we will get either 0 or Inf when taking the exponential of the two other numbers.

Do you see any way to achieve this goal nevertheless? I think it's impossible, but I'd like to be sure.

score 4 · Accepted Answer · answered Aug 01 '23 at 00:26

Normally we could do this without loss of generality by subtracting the maximum value (equivalent to dividing through by the largest probability) but this isn't going to help enough here:

logprobs  <- c(-316480, -85041, -1197042)
lp2 <- logprobs - max(logprobs)

Let's see how big these probabilities would be relative to the max:

lp2/log(10)
[1] -100512.7       0.0 -482935.9

So the second-most-probable outcome is many orders of magnitude smaller than the most probable outcome (larger than a googol, although less than a googolplex). There is no way, before the sun burns out, that you're going to sample anything but the most probable outcome (i.e., category 2).

Double-checking with a simple example to check my math:

probs <- log(c(0.0001, 0.01, 0.000001))
> probs
[1]  -9.21034  -4.60517 -13.81551
> (probs - max(probs))/log(10)
[1] -2  0 -4

This correctly concludes that the first and third outcomes are two and four orders of magnitude less probable than the middle.

You're right, thanks. I had the idea to subtract one value but I thought this was a dead way. — Stéphane Laurent, Aug 01 '23 at 05:52

score 1 · Answer 2 · answered Aug 01 '23 at 07:16

1

I think @Ben Bolker has figured out the possibility when with your logprobs.

If you have some relative "close" logprobs (in terms of the order of magnitude), probably you can optimize like below to find a factor alpha like below to have the desired probability p, e.g.,

logprob <- c(-1, -8.5, -100.9)

f <- function(x) abs(sum(exp(logprob * x)) - 1)

alpha <- optimise(f, c(0, 1), tol = .Machine$double.eps)

p <- exp(logprob * alpha$minimum)

and you will see

> p
[1] 8.182418e-01 1.817582e-01 1.620909e-09

> sum(p)
[1] 1

answered Aug 01 '23 at 07:16

ThomasIsCoding

96,636
9
24
81

Thanks. But we don't need to normalize with `rmultinom`. – Stéphane Laurent Aug 01 '23 at 07:33
@StéphaneLaurent Yes, I know you need the proportional relationship only for that. I was just used `=1` to set up the objective function `f`, we actually can use any other positive number. – ThomasIsCoding Aug 01 '23 at 08:27

Sample multinomial distribution from the log-probabilities

2 Answers2