JAGS 和 R：获得特定 x 的后验预测分布

Question

我正在尝试从 Jags 中的简单线性回归中获得指定 x 值的后验预测分布。通过将此示例（来自 https://biometry.github.io/APES//LectureNotes/StatsCafe/Linear_models_jags.html）应用到我自己的数据，我可以使回归本身发挥作用。我在这里提供了一小部分数据，这样代码也可以在这里工作。

library(rjags)
library(R2jags)

#create data
dw=c(-15.2,-13.0,-10.0,-9.8,-8.5,-8.5,-7.7,-7.5,-7.2,-6.1,-6.1,-6.1,-5.5,-5.0,-5.0,-5.0,-4.5,-4.0,-2.0,-1.0,1.3)
phos=c(11.8,12.9,15.0,14.4,17.3,16.1,20.8,16.6,16.2,18.2,18.8,19.2,15.6,17.0,18.9,22.1,18.9,22.8,21.6,20.5,21.1)

#convert to list
jagsdwphos=list(dw=dw,phos=phos,N=length(phos))

#write model function for linear regression
lm1_jags <- function(){
    # Likelihood:
    for (i in 1:N){
    phos[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
    mu[i] <- intercept + slope * dw[i]
  }
    # Priors:
  intercept ~ dnorm(0, 0.01)
  slope ~ dnorm(0, 0.01)
  sigma ~ dunif(0, 100) # standard deviation
  tau <- 1 / (sigma * sigma)
}

#specifiy paramters of MCMC sampler, choose posteriors to be reported and      run the jags model
#set initial values for MCMC
init_values <- function(){
  list(intercept = rnorm(1), slope = rnorm(1), sigma = runif(1))
}
#choose paramters to report on
params <- c("intercept", "slope", "sigma")
#run model in jags
lm_dwphos <- jags(data = jagsdwphos, inits = init_values, parameters.to.save = params, model.file = lm1_jags,
          n.chains = 3, n.iter = 12000, n.burnin = 2000, n.thin = 10, DIC = F)

除了这个回归之外，我还想得到特定 phos 值的后验预测分布的输出，但我无法让它与我编写的这个简单示例一起使用。我在这里 https://doingbayesiandataanalysis.blogspot.com/2015/10/posterior-predicted-distribution-for.html 找到了一个教程，并尝试像这样实现它：

#create data
dw=c(-15.2,-13.0,-10.0,-9.8,-8.5,-8.5,-7.7,-7.5,-7.2,-6.1,-6.1,-6.1,-5.5,-5.0,-5.0,-5.0,-4.5,-4.0,-2.0,-1.0,1.3)
    phos=c(11.8,12.9,15.0,14.4,17.3,16.1,20.8,16.6,16.2,18.2,18.8,19.2,15.6,17.0,18.9,22.1,18.9,22.8,21.6,20.5,21.1)
#specifiy phos values to use for posterior predictive distribution
phosprobe=c(14,18,22)

#convert to list
jagsdwphos=list(dw=dw,phos=phos,N=length(phos),xP=phosprobe)

#write model function for linear regression
lm1_jags <- function(){
  # Likelihood:
  for (i in 1:N){
    phos[i] ~ dnorm(mu[i], tau) # tau is precision (1 / variance)
    mu[i] <- intercept + slope * dw[i]
  }
  # Priors:
  intercept ~ dnorm(0, 0.01) # intercept
  slope ~ dnorm(0, 0.01) # slope
  sigma ~ dunif(0, 100) # standard deviation
  tau <- 1 / (sigma * sigma) # sigma^2 doesn't work in JAGS
  nu <- nuMinusOne+1
  nuMinusOne ~ dexp(1/29.0)
  #prediction
  for(i in 1:3){
    yP ~ dt(intercept+slope*xP[i],tau,nu)
  }
}

#specifiy paramters of MCMC sampler, choose posteriors to be reported and run the jags model
#set initial values for MCMC
init_values <- function(){
  list(intercept = rnorm(1), slope = rnorm(1), sigma = runif(1))
}
#choose paramters to report on
params <- c("intercept", "slope", "sigma","xP","yP")
#run model in jags
lm_dwphos <- jags(data = jagsdwphos, inits = init_values, parameters.to.save = params, model.file = lm1_jags,
          n.chains = 3, n.iter = 12000, n.burnin = 2000, n.thin = 10, DIC = F)

但我收到以下错误消息：

错误 jags.model(model.file, data = data, inits = init.values, n.chains = n.chains, : RUNTIME ERROR: Compilation error on第 14 行。尝试重新定义节点 yP[1]

我承认我不太明白在我使用的那个例子中预测是如何实现的，也找不到关于 nu 到底是什么或这些数字从何而来的解释。所以我认为这是我在适应我的例子时犯了一些错误，但这是我能找到的 Jags 中唯一的教程，它给出了探测 x 的 y 值的整个分布，而不仅仅是平均值。

如有任何帮助或解释，我将不胜感激。

谢谢！

Answer 1

发生此错误是因为您没有索引 yP。你已经像这样写了这个循环：

  #prediction
  for(i in 1:3){
    yP ~ dt(intercept+slope*xP[i],tau,nu)
  }

随着 i 从 1 移动到 3，元素 yP 正在被覆盖。您需要像使用 xP.

一样对其进行索引

  #prediction
  for(i in 1:3){
    yP[i] ~ dt(intercept+slope*xP[i],tau,nu)
  }

JAGS 和 R：获得特定 x 的后验预测分布

JAGS and R: Obtain posterior predictive distribution for specific x

r

linear-regression

jags

r2jags