1

我尝试使用迭代过程来估计模型选择的贝叶斯设置中的边际似然。如果您对具体细节感兴趣，请参阅this或this paper。我与数字问题斗争了很长一段时间，所以我想我不妨在这里试一试，试试运气。我试图使下面的内容尽可能简短、可重复和通用。

那个设定

目标是运行表单的迭代过程

在此处输入图像描述

获得对我的模型的边际似然y的估计（这个迭代过程将很快收敛到y的某个值）。您在公式中看到的所有其他变量都是标量或已计算的参数可能性，因此它们是固定值。每次迭代中唯一变化的项是y的值。所以基本上我有长度为 L 的向量p_l和q_l以及长度为 M 的向量p_m和q_m。到目前为止一切都很好。

问题

不幸的是，上面公式中使用的似然值不是对数似然（我计算并存储在向量中），而是实际的似然。这是我遇到的数字问题的根源。您可能知道对数量级很大的负对数似然求幂的老问题。我的对数似然性非常负，以至于 exp(log似然性) 几乎总是为 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)
}

谢谢！

4

1 回答 1

0

我认为最好在日志中工作。我会一起工作

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

等你有表格的条款总和

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

t 的 log lt 为

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

在哪里

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

注意 mM 不是正数；如果它是非常负的 exp(mM) 可能会下溢到 0，但这没关系，因为我们将它添加到 1。

然后我们有形式的总和

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

同样，它的日志 ls 是

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

其中 L 是 l[] 的 (a) 最大值。我们现在可以评估 exp 而不必担心溢出。

于 2018-06-12T08:58:14.257 回答