bayesian - 在 JAGS 中运行逻辑模型 - 你可以矢量化而不是循环个别案例吗？

Question

我对 JAGS 还很陌生，所以这可能是个愚蠢的问题。我正在尝试在 JAGS 中运行一个模型，该模型预测一维随机游走过程在穿过边界 B 之前穿过边界 A 的概率。该模型可以通过以下逻辑模型分析求解:

Pr(A,B) = 1/(1 + exp(-2 * (d/sigma) * theta))

其中“d”是平均漂移率(正值表示向边界 A 漂移)，“sigma”是该漂移率的标准差，“theta”是起点和边界之间的距离(假设两个边界都相等)。

我的数据集由 50 名参与者组成，每名参与者提供 1800 个观察结果。我的模型假设 d 由观察到的环境变量(我将其称为“x”)和将 x 与 d 相关的加权系数(我将其称为“beta”)的特定组合确定。因此，存在三个参数:beta、sigma 和 theta。我想为每个参与者估计一组参数。我的意图是最终运行一个分层模型，其中组级别参数影响个人级别参数。然而，为简单起见，在这里我将只考虑一个模型，在该模型中我为一个参与者估计一组参数(因此该模型不是分层的)。

我在 rjags 中的模型如下:

model{
for ( i in 1:Ntotal ) {   
 d[i] <- x[i] * beta
 probA[i] <- 1/(1+exp(-2 * (d[i]/sigma) * theta )   )
 y[i] ~ dbern(probA[i])    
 }

 beta ~ dunif(-10,10)
 sigma ~ dunif(0,10)
 theta ~ dunif(0,10)
}

此模型运行良好，但需要很长时间才能运行。我不确定 JAGS 如何执行代码，但如果此代码在 R 中运行，效率会相当低，因为它必须遍历案例，为每个案例单独运行模型。因此，随着样本量的增加，运行分析所需的时间会迅速增加。我有一个相当大的样本，所以这是一个问题。

有没有办法将这段代码向量化，以便它可以一次计算所有数据点的可能性？例如，如果我将其作为一个简单的最大似然模型来运行。我会对模型进行矢量化，并计算参与者提供的所有 1800 个案例的给定特定参数值的数据概率(因此不需要 for 循环)。然后我会记录这些可能性并将它们加在一起，为参与者给出的所有观察结果提供一个单一的对数可能性。这种方法可以节省大量时间。在 JAGS 中有没有办法做到这一点？

---

编辑

感谢您的回复，并感谢您指出我展示的模型中的参数可能无法识别。我应该指出该模型是一个简化版本。完整模型如下:

model{
  for ( i in 1:Ntotal ) {
    aExpectancy[i] <- 1/(1+exp(-gamma*(aTimeRemaining[i] - aDiscrepancy[i]*aExpectedLag[i]) ) )
    bExpectancy[i] <- 1/(1+exp(-gamma*(bTimeRemaining[i] - bDiscrepancy[i]*bExpectedLag[i]) ) )
    aUtility[i] <- aValence[i]*aExpectancy[i]/(1 + discount * (aTimeRemaining[i]))
    bUtility[i] <- bValence[i]*bExpectancy[i]/(1 + discount * (bTimeRemaining[i]))
    aMotivationalValueMean[i] <- aUtility[i]*aQualityMean[i]
    bMotivationalValueMean[i] <- bUtility[i]*bQualityMean[i]
    aMotivationalValueVariance[i] <- (aUtility[i]*aQualitySD[i])^2 + (bUtility[i]*bQualitySD[i])^2
    bMotivationalValueVariance[i] <- (aUtility[i]*aQualitySD[i])^2 + (bUtility[i]*bQualitySD[i])^2
    mvDiffVariance[i] <- aMotivationalValueVariance[i] + bMotivationalValueVariance[i]
    meanDrift[i] <- (aMotivationalValueMean[i] - bMotivationalValueMean[i])
    probA[i] <- 1/(1+exp(-2*(meanDrift[i]/sqrt(mvDiffVariance[i])) *theta ) )
    y[i] ~ dbern(probA[i])
  }

在这个模型中，估计的参数是theta、discount和gamma，这些参数是可以恢复的。当我根据单个参与者 (Ntotal = 1800) 的观察结果运行模型时，模型运行大约需要 5 分钟，这完全没问题。然而，当我在整个样本上运行模型时(45 名参与者 x 每个 1800 个案例 = 78,900 个观察值)，我已经运行了 24 小时并且它还不到 50%。这看起来很奇怪，因为我预计它只需要 45 倍的时间，所以最多 4 或 5 个小时。我错过了什么吗？

Answer 1

我希望我没有误读这种情况(如果我误读了，我之前会道歉)，但你的问题似乎来自对 JAGS 工作原理(或 WinBUGS 或 OpenBUGS 就此而言)的概念性误解。

你的程序实际上没有运行，因为你写的不是用编程语言写的​​。所以向量化没有帮助。

您只写了模型的描述，因为 JAGS 的语言是描述性语言。

JAGS 读取您的模型后，它会组装一个转换矩阵来运行 MCMC，其平稳分布是给定(观察到的)数据的参数的后验分布。 JAGS 不会对您的程序执行任何其他操作。

您一直在等待程序运行的所有时间实际上是在等待(并希望)达到您的 MCMC 的放松时间。

因此，导致您的程序运行时间过长的原因是生成的转换矩阵必须具有不良的松弛属性或类似的东西。

这就是为什么对只读取和运行一次的程序进行矢量化几乎没有帮助的原因。

那么，您的问题出在其他地方。

希望对您有所帮助，如果没有，请见谅。

祝一切顺利