对于小 (1e-4) 参数值，无法用 rethinking::quap 拟合贝叶斯模型

Question

我正在尝试使用 rethinking 软件包版本 2.13 计算 p 的后验概率：

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)
)

它抛出一个错误：

Error in rethinking::quap(alist(n ~ dbinom(N, p), p ~ dnorm(1e-04, 5e-05)), : non-finite finite-difference value [1] Start values for parameters may be too far from MAP. Try better priors or use explicit start values. If you sampled random start values, just trying again may work. Start values used in this attempt: p = 1e-04

start其实就是MAP
它适用于更大的 p 值。
如何解决问题？

Answer 1

tl;dr 您正在运行解决数值逼近问题；添加control=list(ndeps=1e-6)调整有限差分容差是解决问题的最快方法。

或者，您可以更改模型以适应对数刻度 p：

alist(
    n ~ dbinom2(N, exp(logp)),
    logp ~ dnorm(-4, 2)
)

更一般地说，如果您正在处理 p 可能大于 0.5 的问题，您可能应该适合对数几率或对数比例而不是对数比例，使用 plogis(logit_p) 而不是 exp(log_p)。（我使用 log/exp 而不是 logit (qlogis)/plogis 因为我首先想到了它，并且因为它对大多数人来说更熟悉。）

我尝试的两种方法没有工作（我认为应该有；我需要进一步挖掘以找出原因）：(1) 指定预期的比例通过 control=list(parscale=1e-4) 的参数； (2) 使用method="L-BFGS-B", lower=0设置概率参数的下限。

可悲的现实是，用于可能性、后验等的数值评估的工具非常脆弱，以至于您很快就会达到需要了解幕后情况的地步。

我们可以通过编写 dbinom2() 函数来诊断发生了什么，该函数只是 dbinom:

的“嘈杂”包装器

dbinom2 <- function(x, size, prob, log) {
    r <- dbinom(x, size, prob, log)
    cat(x,size,prob,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，然后 prob=-9e-04（负概率值），这是运行有麻烦了。

该错误消息提供了更多有关正在发生的事情的线索（如果您已经对事情的幕后工作方式了解很多）； quap() 默认使用“BFGS”方法（它直接将其传递给 optim()，这意味着它正在尝试使用 衍生物 或 gradients 找到 MAP 估计的函数。
因为我们没有明确的梯度表达式，R 做了次好的事情（这很糟糕，但对于简单的问题工作正常），即尝试通过 [=63 来估计梯度=]finite differences — 基本上将导数的定义复制为 lim(dx→0) (f(x+dx)-f(x))/dx，但使用较小的（“有限”）dx 值而不是尝试取极限。
问题出现在决定 dx 时。一路进入optim()(?optim())的帮助页面，你可以看到：

‘ndeps’ A vector of step sizes for the finite-difference approximation to the gradient, on ‘par/parscale’ scale. Defaults to ‘1e-3’.

这意味着 optim() 正在使用 dx=0.001，这就是您的问题所在。它会尝试 f(x-dx) 以及 f(x+dx)，并且由于 x-dx 小于 0，这会中断计算。使用 ndeps=1e-6 指定 dx 应该更小（足够小以至于 optim() 不会尝试评估负概率的对数似然）。

对于小 (1e-4) 参数值，无法用 rethinking::quap 拟合贝叶斯模型

Cannot fit bayesian model with rethinking::quap for small (1e-4) parameter values

r

bayesian

numerical-methods