3

我对 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 个小时。我错过了什么吗？

4

3 回答 3

5

我希望我没有误读这种情况（如果我错了，我之前道歉），但您的问题似乎来自对 JAGS 工作原理（或 WinBUGS 或 OpenBUGS 的工作原理）的概念性误解。

你的程序实际上并没有运行，因为你写的不是用编程语言写的。所以矢量化将无济于事。

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

一旦 JAGS 读取了您的模型，它就会组装一个转换矩阵来运行一个 MCMC，其平稳分布是给定（观察到的）数据的参数的后验分布。JAGS 不会对您的程序做任何其他事情。

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

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

这就是为什么向量化一个只被读取和运行一次的程序将没有什么帮助。

所以，你的问题出在其他地方。

我希望它有帮助，如果没有，对不起。

一切顺利。

于 2015-09-10T13:50:59.393 回答
2

您不能以与 R 中相同的方式矢量化，但如果您可以使用相同的概率表达式（即常见的 d[i]）对观察进行分组，那么您可以使用二项分布而不是伯努利分布，这将有很大帮助。如果每个观察值都有一个唯一的 d[i] ，那么恐怕你会被卡住。

另一种选择是查看 Stan，它对于像您这样的大型数据集通常更快。

马特

于 2015-09-07T07:32:17.543 回答
-1

感谢您的回复。是的，您提出了一个很好的观点，即我展示的模型中的参数可能无法识别。

我应该指出该模型是简化版本。完整模型如下：

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 ~ dunif(0,10)  
  discount ~ dunif(0,10)
  gamma ~ dunif(0,1)
}

在这个模型中，估计的参数是theta、discount和gamma，并且这些参数可以被恢复。

当我对单个参与者 ( Ntotal = 1800) 的观察结果运行模型时，模型运行大约需要 5 分钟，这完全没问题。

但是，当我在整个样本上运行模型时（45 名参与者 X 1800 个案例 = 78,900 个观察值），我已经运行了 24 小时，而且还不到 50%。

这看起来很奇怪，因为我预计它只需要 45 倍的时间，所以最多 4 或 5 个小时。我错过了什么吗？

于 2015-09-06T11:30:45.943 回答