5

我试图在 R 中计算这个后验分布。问题是分子太小，它是一堆 dbern(p_i, y_i) < 1 的乘积。（我的 n 约为 1500）。因此，R 吐出 0，所有 \theta 的后验值也为 0。

澄清一下，每个 y_i 都有自己的 p_i，这些 p_i 一起为 n y 组成一个包含 n 个元素的向量。每个 theta 都有自己的 n 元素向量 p_i。

在此处输入图像描述

一个可重复的例子（分子）

p <- sample(seq(0.001,0.999,by=0.01), 1500, replace=T)
y <- sample(c(0,1), 1500, replace=T)
dbern(y, p) # 1500-element vector, each element < 1
prod(dbern(y, p)) # produces 0
exp(sum(log(dbern(y, p)))) # produces 0

编辑（上下文）：我正在做贝叶斯变化点分析（jstor.org/stable/25791783 - Western and Kleykamp 2004）。与论文中的连续 y 不同，我的 y 是二进制的，所以我使用的是 Albert 和 Chib (1993) 中的数据增强方法。使用这种方法，y 的可能性是伯努利，p = cdf-normal(x'B)。

那么 p 如何取决于 theta 呢？这是因为θ是变化点。其中一个 x 是时间虚拟变量——例如，如果 theta=10，那么对于第 10 天之后的所有观测值，时间虚拟变量 = 1，对于第 10 天之前的所有观测值，时间虚拟变量 = 0。

因此，p 取决于 x，x 取决于 theta - 因此，p 取决于 theta。

我需要上述数量，因为它是吉布斯采样中 theta 的完整条件。

4

3 回答 3

5

解决此类精度问题的一种方法是在日志空间中工作。但是，这会引入分母的对数和积，这通常可能会很痛苦。

如果您出于优化的目的计算后验，请注意您可以完全删除分母：您无需进行归一化即可找到argmax.

于 2013-04-21T02:44:37.677 回答
3

好吧，我运行了您的示例，您（正如预期的那样）得到了，0因为有一个0

p <- sample(seq(0,1,by=0.01), 1500, replace=T)
y <- sample(c(0,1), replace=T)
x <- dbern(y, p)
any(x == 0)
## [1] TRUE
于 2013-04-21T06:37:07.903 回答
1

我也在 Cross Validated 上问了这个问题，glen_b在下面给了我这个（经过测试的）解决方案：

这是计算各种模型的可能性的常见问题；通常做的事情是在日志上工作，并使用一个通用的比例因子将值带入一个更合理的范围内。

在这种情况下，我建议：

第 1 步：选择一个相当“典型”的 θ，θ0。将通用项的分子和分母的公式除以 θ=θ0 的分子，以得到不太可能下溢的东西。

步骤2：在对数尺度上工作，这意味着分子是对数差异和的exp，分母是对数差异和的exp的总和。

注意：如果您的任何 p 是 0 或 1，请将它们单独拉出并且不要记录这些术语；它们很容易按原样评估！

分子中的常用项在大小上往往会更适中，因此在许多情况下，分子和分母都是相对合理的。

如果分母中有一系列尺寸，请先将较小的尺寸相加，然后再添加较大的尺寸。

如果一个或几个术语占主导地位，您应该将注意力集中在使那些相对准确的计算上。

于 2013-04-22T06:15:23.057 回答