r - 带有 R 的 JAGS：模拟数据并一次拟合多个模型

Question

我有三个非常相似的模型，我想通过 LOO 和 WAIC 比较哪个模型具有更高的预测准确性。我想多次模拟数据并取它们的对数似然的平均值。手动执行此操作将花费很长时间，因为我需要为每个模拟执行以下步骤： 1. 通过一个模型模拟数据；2. 保存第一个模型的数据和可能性；3.将数据拟合到第二个模型并保存第二个模型的可能性；3. 将数据拟合到第三个模型并保存对数似然。然后运行第二个模拟并重复所有四个步骤。

有什么方法可以同时模拟数据，例如多次模拟数据，每次都将数据拟合到所有三个模型并保存它们的对数似然？或者，至少，对于每个模拟，将数据拟合到所有三个模型一次？（比如将 1 个数据块和 2 个模型块放在一个 txt 文件中？）

以下是我的三个模型中的两个（两个模型共享相同的数据列表）。

loglik是计算对数似然。

注意：我使用 R2jags

模型1：例如，使用这个模型来模拟数据

# simulated data
data{
for (i in 1:Nsubjects)
{
   for (m in 1:Ntrials)
{

    expectancy [i,m] ~ dnorm (b_true [i] * v_true[i,m] + (10 - b_true [i]),tau_true[i]) T(0,10)
}
      for (j in 2: Ntrials)
{
 
    v_true [i,j] <- v_true [i,j-1] + alpha_true [i] * pe_true [i,j-1]

 
    pe_true [i,j-1] <- shock [i,j-1] - v_true[i,j-1]
}

    sigma_true [i] ~ dunif (0,1)


    tau_true [i] <- 1/pow(sigma_true [i],2)  


    alpha_true [i] ~ dbeta(1,1)


    b_true [i] ~ dunif (0,10)
}
}

# model
model{
  # data
for (i in 1:Nsubjects)
{
   for (m in 1:Ntrials)
{

      expectancy [i,m] ~ dnorm (b [i] * v[i,m] + (10 - b[i]),tau[i])


      predk [i,m] ~ dnorm (b [i] * v[i,m] + (10 - b[i]),tau[i])

      # log-likelihood
      loglik[i,m] <- logdensity.norm(expectancy [i,m], b [i] * v[i,m] + (10 - b[i]),tau[i])

}
    for (j in 2:Ntrials)
    {
 
      v [i,j] <- v [i,j-1] + alpha [i]  * pe [i,j-1]


      pe [i,j-1] <- shock [i,j-1] - v [i,j-1]


  # priors
    }
    sigma[i] ~ dunif (0,5)
    tau [i] <- 1/pow(sigma [i],2) 

    b [i] ~ dnorm (b_mu,b_tau) T(0,10)
    alpha [i] ~ dnorm (alpha_mu,alpha_tau) T(0,1) 
}

   # hyperpriors
     alpha_mu ~ dbeta (1,1)
     alpha_tau ~ dscaled.gamma (2,1)
     alpha_sd <- 1/sqrt (alpha_tau)
     b_mu ~ dunif (0,10)
     b_tau ~ dscaled.gamma (2,1)
     b_sd <- 1/sqrt (b_tau)
}

模型 2

model{
for (i in 1:Nsubjects)
{
  for (m in 1:Ntrials)
{
 
      expectancy [i,m] ~ dnorm (b [i] * v [i,m] + (10 - b [i]),tau[i])

      # posteiror predictive
      predk [i,m] ~ dnorm (b [i] * v [i,m] + (10 - b [i]),tau[i])

      # log-likelihood
      loglik[i,m] <- logdensity.norm (expectancy [i,m],b [i] * v [i,m] + (10 - b [i]),tau[i])
}
    for (j in 2:Ntrials)
    {
  
      v [i,j] <- v [i,j-1] + alpha [i,j-1]  * pe [i,j-1]


      alpha [i,j] <- gamma [i] * abs (pe [i,j-1]) + (1-gamma [i]) * alpha [i,j-1]   


      pe [i,j-1] <- shock [i,j-1] - v [i,j-1]

    }
  # priors
    sigma[i] ~ dunif (0,5)
    tau [i] <- 1/pow(sigma [i],2) 

    b [i] ~ dnorm (b_mu,b_tau) T(0,10)
    gamma [i] ~ dnorm (gamma_mu,gamma_tau) T(0,1)
    alpha [i,1] ~ dnorm (alpha_mu,alpha_tau) T(0,1) 
}

   # hyperpriors
     alpha_mu ~ dbeta (1,1)
     alpha_tau ~ dscaled.gamma (2,1)
     alpha_sd <- 1/sqrt (alpha_tau)
     b_mu ~ dunif (0,10)
     b_tau ~ dscaled.gamma (2,1)
     b_sd <- 1/sqrt (b_tau)
     gamma_mu ~ dbeta(1,1)
     gamma_tau ~ dscaled.gamma (2,1)
     gamma_sd <- 1/sqrt (gamma_tau)
}

标签： rjagsr2jags