问题描述
我正在尝试使用 p 包版本 2.13 计算 rethinking 的后验：
rethinking::quap(
alist(
n ~ dbinom(N,p),p ~ dnorm(1e-4,0.5e-4)
),data = list(n = 10,N = 10e4),start = list(p = 10 / 10e4)
)
重新思考时出错::quap(alist(n ~ dbinom(N,p ~ dnorm(1e-04,5e-05)),: 非有限有限差分值 [1] 参数的起始值可能与 MAP 相差太远。 尝试更好的先验或使用明确的起始值。 如果您对随机起始值进行采样，只需再试一次即可。 本次尝试中使用的起始值： p = 1e-04
start 实际上是 MAP。
它适用于更大的 p 值。
如何解决问题？
解决方法
tl;dr 您遇到了数值近似问题；添加control=list(ndeps=1e-6)来调整有限差分容差是解决问题的最快方法。
或者，您可以更改模型以在对数刻度上拟合 p：
alist(
n ~ dbinom2(N,exp(logp)),logp ~ dnorm(-4,2)
)
更一般地，如果您正在处理 p 可能大于 0.5 的问题，您可能应该使用对数赔率或 logit 标度而不是对数标度，使用 {{1 }} 而不是 plogis(logit_p)。 （我用 log/exp 而不是 logit (qlogis)/plogis 是因为我最先想到的，也因为它对大多数人来说更熟悉。）
我尝试了两种没有的方法（我认为应该有；我需要进一步挖掘以找出原因）：（1）通过{指定参数的预期比例{1}}; (2) 使用exp(log_p)设置概率参数的下限。
可悲的现实是，对可能性、后验等进行数值评估的工具非常脆弱，您很快就会达到需要了解幕后情况的程度。
我们可以通过编写一个 control=list(parscale=1e-4) 函数来诊断发生了什么，该函数只是一个 method="L-BFGS-B",lower=0 的“嘈杂”包装器：
dbinom2()然后用 dbinom 代替 dbinom2 <- function(x,size,prob,log) {
r <- dbinom(x,log)
cat(x,r,"\n")
return(r)
}
重新运行估计，我们得到：
dbinom2(...)在这种情况下，dbinom(...) 和 10 1e+05 1e-04 -2.078512
10 1e+05 0.0011 -78.1496
10 1e+05 -9e-04 NaN
值保持不变；搜索算法首先尝试x，然后是size，然后是prob=1e-4（负概率值），这就是它遇到麻烦的地方。
- 错误消息提供了更多关于正在发生的事情的线索（如果您已经对引擎盖下的工作原理有相当多的了解）；
prob=0.0011默认使用“BFGS”方法（它直接将其传递给prob=-9e-04，这意味着它正在尝试使用 derivatives 或 gradients em> 的函数来找到 MAP 估计。 - 因为我们没有一个明确的梯度表达式，R 做了下一个最好的事情（这很糟糕，但对于简单的问题可以正常工作），即尝试通过有限差分来估计梯度 — 基本上将导数的定义复制为
quap()，但使用小的（“有限”）optim()值而不是尝试取极限。 - 在决定
lim(dx→0) (f(x+dx)-f(x))/dx时会出现问题。进入dx(dx) 的帮助页面方式，您可以看到：
‘ndeps’ 有限差分的步长向量 在“par/parscale”尺度上近似梯度。 默认为“1e-3”。
这意味着 optim() 正在使用 ?optim()，这就是您的问题所在。它尝试 optim() 和 dx=0.001，并且由于 f(x-dx) 小于 0，这会破坏计算。使用 f(x+dx) 指定 x-dx 应该更小（足够小，以便 ndeps=1e-6 不会尝试评估负概率的对数似然）。