可能性太小的乘积 - R只给出0
我试图在R中计算这个后验分布.问题是分子是一堆dbern(p_i,y_i)<1的乘积,它太小了.(我的约是1500).因此,R吐出0,所有\ theta的后验值也为0.
为了澄清,每个y_i都有自己的p_i,这些p_i一起构成n个元素的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和Kleykamp 2004).与论文中的连续y不同,我的y是二元的,所以我在Albert和Chib(1993)中使用数据增强方法.使用该方法,y的可能性是伯努利,p = cdf-normal(x'B).
那么p如何依赖于theta？这是因为θ是变化点.其中一个x是时间虚拟 - 例如,如果θ= 10,那么对于第10天之后的所有观察,时间dummy = 1,并且对于第10天之前的所有观察,= 0.
因此,p取决于x,x取决于θ - 因此,p取决于θ.
我需要上面的数量,因为这是吉布斯采样中theta的完全条件.
解决这类精度问题的一种方法是在日志空间中工作.然而,这引入了分母的对数和积,这通常可能是痛苦的.
如果您为了优化而计算后验,请注意您可以完全删除分母:您不需要标准化来查找argmax.
我也在 Cross Validated 上问了这个问题，glen_b给了我下面这个（经过测试的）解决方案：
\n\n\n\n
这是计算各种模型的可能性的常见问题；通常要做的事情是处理日志，并使用通用的缩放因子将值带入更合理的范围。
\n\n在这种情况下，我建议：
\n\n步骤1：选择一个相当“典型”的\xce\xb8、\xce\xb80。将通项的分子和分母公式除以 \xce\xb8=\xce\xb80 的分子，以获得下溢可能性较小的值。
\n\n步骤2：使用对数尺度，这意味着分子是对数差值之和的exp之和，分母是对数差值之和的exp之和。
\n\n注意：如果您的任何 p 为 0 或 1，请将它们单独取出，并且不要记录这些项；它们很容易按原样进行评估！
\n\n分子中的常用项的大小往往比较适中，因此在许多情况下分子和分母都相对合理。
\n\n如果分母有多个大小范围，请先将较小的值相加，然后再将较大的值相加。
\n\n如果一个或几个术语占主导地位，您应该集中精力使这些术语的计算相对准确。
\n