Численные проблемы в оценке предельного правдоподобия

Question

Я пытаюсь использовать итерационную процедуру для оценки предельной вероятности в байесовской настройке для выбора модели.Если вас интересует конкретика, см., Например, эту или эту бумагу.Я уже давно борюсь с числовыми проблемами, поэтому я подумал, что могу попробовать и здесь и попытать счастья.Я попытался сделать следующее ниже как можно более коротким, воспроизводимым и общим.

Параметр

Цель состоит в том, чтобы запустить итерационный процесс в форме

, чтобы получить оценку предельной вероятности y моей модели (этот итерационный процесс довольно быстро сходится к некоторому значениюдля у ).Все остальные переменные, которые вы видите в формуле, являются скалярами или вероятностями параметров, которые уже вычислены и, следовательно, являются фиксированными значениями.Единственный член, изменяющийся в каждой итерации, это значение y .Так что в основном у меня есть векторы p_l и q_l длины L и вектор p_m и q_m длины M. Пока все хорошо.

Проблема

Значения правдоподобия, использованные в приведенной выше формуле, к сожалению, являются не логарифмическими правдоподобиями (которые я вычисляю и сохраняю в векторах), а фактическими правдоподобиями.Это источник численных проблем, которые у меня есть.Вы, вероятно, знаете старую проблему возведения в степень отрицательных логарифмических вероятностей, которые являются большими по величине.Вероятности моего журнала настолько отрицательны, что exp (вероятность журнала) почти всегда будет 0. Есть некоторые похожие вопросы с ответами, уже находящимися в сети, но ни одно из этих решений не помогло мне.

Что япробовал

Я думаю, что может быть многообещающим, чтобы использовать тот факт, что exp (log (x)) = x, и расширить дробь, чтобы вы могли переписать приведенную выше формулу как

, где C - постоянная по вашему выбору.Доказательство, которое следует той же идее, см. В приложении к этой статье .Если можно найти подходящее значение C , которое делает слагаемые в сумме управляемыми, проблема решается.Однако в моем случае p_l , q_l , p_m и q_m сильно отличаются по величине, поэтому независимо от того, какое значение C Я пытаюсь вычесть или добавить, я снова получу переполнение или переполнение.Таким образом, сейчас я действительно не знаю, как действовать дальше.Я ценю любые комментарии и советы.

Код и данные

Вы можете найти некоторые примеры данных (то есть векторы с вероятностями записи) здесь ,Код для итеративного процесса:

L <- 1000
M <- 5000
y <- numeric(length=100)
eval_downer <- numeric(length=L)
eval_upper <- numeric(length=M)
y[1] <- 0.5

for(t in 2:100){ 

  for(m in 1:M){
    up.m <- q_m[m]
    down.m <- (L * q_m[m]) + (M * p_m[m] / y[t-1])
    eval_downer[m]  <- up.m / down.m 
  }

  for(l in 1:L){
    up.l <- p_l[l]
    down.l <- (L * q_l[l]) + (M * p_l[l] / y[t-1])
    eval_upper[l]  <- up.l / down.l
  }

  upper <- mean(eval_upper)
  downer <- mean(eval_downer)

  y[t]  <-   upper / downer

  print(t)
}

Спасибо!

Answer 1

Лучше всего работать в логах, я думаю.Я буду работать с

log( (M * p_m[m] / y[t-1]) 

и т. Д.У вас есть суммы терминов в форме

t = exp( a)/(exp(b) + exp(c)) = 1 / (exp(b-a) + exp(c-a))

Лог l t

lt = (a-M) -  log1p( 1 + exp(m-M))

где

m = min(b,c), M = max(b,c)

Обратите внимание, что мМ не является положительным;если он очень отрицательный, exp (mM) может опуститься до 0, но это нормально, так как мы добавляем его к 1.

Тогда у нас есть суммы вида

s = Sum{ exp(l[i])}

Аналогично, его лог ls равен

ls = L + log( Sum{ 1 + exp(l[i]-L)})

, где L - максимум (a) l [].Теперь мы можем оценивать опыт, не опасаясь переполнения.