2

我试图从 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/postterior-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 line 14. 尝试重新定义节点 yP[1]

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

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

谢谢！

4

1 回答 1

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)
  }
于 2018-07-16T13:27:39.703 回答