Как вести себя с логом суммы более двух сверхмалых вероятностей

Question

По какой-то злой причине мне нужно вычислить логарифм суммы 500 сверхмалых вероятностей, каждый член которых вычисляется как

dmvnorm(X[,i], mean=rep(0,3), sigma=diag(3))

Иногда приведенные выше коды возвращают 0 из-за недостаточного значения, но при использовании логарифмов хорошо.

dmvnorm(X[,i], mean=rep(0,3), sigma=diag(3), log=TRUE)

Я знаю, что могу математически обработать два термина: log(p1 + p2) = log(p2) + log(1 + p1/p2). Но можем ли мы обобщить этот подход на большее количество терминов? Кто-нибудь с большим опытом в этом?

---

Кстати, я написал рекурсивную функцию для вычисления этого. Математически это работает. Но я не думаю, что это практично.

MESSY <- function (pv) 
{
  N <- length(pv)
  if (N==1)
    {return(pv[1])}
  else
    {w <- pv[N]
     pv <- pv[1:N-1]-w
     return(w + log(1 + exp(MESSY(pv))))
    }
}

Поскольку некоторые p очень малы, и я могу использовать только w=log(p), у нас есть log(exp(w1)+exp(w2)) = w2 + log(1+exp(w1-w2)) и log(exp(w1)+exp(w2)+exp(w3)) = w3 + log(1 + exp(w1-w3) + exp(w2-w3)) для двух и трех членов.

Answer 1

Эта функция переводится с внутренней функции logspace_add в исходном коде R здесь

logspace_add <- function(logx,logy) { 
    pmax(logx,logy) + log1p(exp(-abs(logx - logy))) 
}

Не обязательно самая эффективная, но вы должны быть в состоянии сделать это для > 2 элемента с помощью Reduce():

logspace_add_mult <- function(x) {
    Reduce(logspace_add, x)
}

Быстрый тест (с использованием значений, которые достаточно велики , а не , чтобы потерять значение, чтобы мы могли сравнить результаты регулярных вычислений и вычислений в лог-пространстве ).

x <- c(1e-4,1e-5,1e-6)
logspace_add_mult(log(x))
## [1] -9.10598
log(sum(x))
## [1] -9.10598

Насколько я знаю, это более или менее стандартный подход к добавлению пространства журналов. Преимущество использования чего-то другого, кроме этой домашней реализации, будет: (1) зрелость кода и тестирование и (2) скорость (по крайней мере, для версии logspace_add_mult; я сомневаюсь, что было бы много преимуществ от нативной C (или что угодно) реализация logspace_add). Пакет Brobdingnag использует аналогичные представления:

library(Brobdingnag)
brob(log(x))
## [1] +exp(-9.2103) +exp(-11.513) +exp(-13.816)
sum(brob(log(x)))
## [1] +exp(-9.106)
log(as.numeric(sum(brob(log(x)))))
## [1] -9.10598

В Python, numpy имеет logaddexp , но это работает только попарно: вы можете использовать functools.reduce() обобщить, как указано выше.

import numpy as np
import functools as f
x = np.array([1e-4,1e-5,1e-6])
f.reduce(np.logaddexp,np.log(x))

Это, вероятно, немного быстрее, чем Reduce() в R.