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

По какой-то злой причине мне нужно вычислить логарифм суммы 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)) для двух и трех терминов.


r probability underflow logarithm
person Paw in Data    schedule 14.01.2020    source источник
comment
Можете ли вы показать небольшой воспроизводимый пример. Кроме того, если это три члена, какая будет формула   -  person akrun    schedule 14.01.2020
comment
Является ли log(p2) + log(1 + p1/p2) более точным, чем log(p1 + p2)? Все еще кажется, что p1/p2 будет проблематично.   -  person Gregor Thomas    schedule 14.01.2020
comment
Какова цель этого?   -  person akash87    schedule 14.01.2020
comment
Это может помочь: statmodeling.stat.columbia.edu /2016/06/11/   -  person Kent Johnson    schedule 14.01.2020
comment
Спасибо! @akrun Я только что отредактировал сообщение.   -  person Paw in Data    schedule 15.01.2020
comment
@Грегор -- восстанови Монику -- Ты прав. Это может быть проблематично. Но некоторые члены вероятности слишком малы, поэтому я должен сделать это другими способами.   -  person Paw in Data    schedule 15.01.2020
comment
@KentJohnson Большое спасибо! Это может быть самым элегантным решением! Однако я не могу найти функцию log_sum_exp( ) или logSumExp( ). Ты знаешь, из какого они пакета?   -  person Paw in Data    schedule 15.01.2020
comment
Я думаю, что они находятся в исходном коде C, лежащем в основе R, но не экспортируются... -L230" rel="nofollow noreferrer">github.com/wch/r-source/blob/ . Переводя на R, logspace_add <- function(logx,logy) { max(logx,logy) + log1p(exp(-abs(logx - logy))) } (это для двух элементов, надо бы подумать, как обобщить...)   -  person Ben Bolker    schedule 15.01.2020
comment
@BenBolker Да, это имеет смысл. Большое спасибо!   -  person Paw in Data    schedule 16.01.2020
comment
@akash87 akash87 Это для вычисления минимальной длины описания, если вы об этом спрашиваете?   -  person Paw in Data    schedule 21.01.2020

Ответы (1)


arrow_upward
2
arrow_downward

Эта функция переведена из внутренней функции 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.

person Ben Bolker    schedule 16.01.2020
comment
Большое спасибо! Этот подход работает, и это намного быстрее. Хотя я получил некоторые неожиданные результаты, поэтому я подозреваю, что для некоторых терминов все еще происходит недополнение, только я этого не вижу. Но я не уверен. Вы бы порекомендовали MATLAB или Python для таких расчетов? - person Paw in Data; 24.01.2020
comment
Я сомневаюсь, что будет очень большая разница в численной стабильности. - person Ben Bolker; 24.01.2020
comment
Я понимаю. Большое спасибо! - person Paw in Data; 24.01.2020