1

我已阅读 JAGS 手册，但没有找到将相同的先验分布分配给 JAGS / R2JAGS 模型中的多个参数的方法。

例如，目前我必须重复很多这样的代码：

reg.model <- function() {
  # Model structure
  for(i in 1:N){
    Y[i] ~ dnorm(mu[i], phi)
    mu[i] <- beta0 + beta1*x1[i] + beta2*x2[i]
  }
  sigma2 <- 1 / phi

  # Priors (Lots of re-typing here for beta0, beta1, beta2)
  phi   ~ dgamma(1,1)
  beta0 ~ dnorm(0, 0.01*phi)
  beta1 ~ dnorm(0, 0.01*phi)
  beta2 ~ dnorm(0, 0.01*phi)
}

如何干燥此代码？

4

1 回答 1

2

如果将自变量作为矩阵传递，将系数放入向量中，然后使用矩阵乘法（使用inprod或%*%）来计算线性预测变量，则可以执行此操作。

例如：

M <- function() {
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], sd^-2)
    mu[i] <- X[i, ] %*% beta
  }

  sd ~ dunif(0, 100)
  for (j in 1:J) {
    beta[j] ~ dnorm(0, 0.0001)
  }
}

我们生成一些虚拟数据，如下：

set.seed(1)
n <- 1000
J <- 3
X <- cbind(1, matrix(runif(n * J), ncol=J))

head(X)
#      [,1]      [,2]       [,3]      [,4]
# [1,]    1 0.2655087 0.53080879 0.8718050
# [2,]    1 0.3721239 0.68486090 0.9671970
# [3,]    1 0.5728534 0.38328339 0.8669163
# [4,]    1 0.9082078 0.95498800 0.4377153
# [5,]    1 0.2016819 0.11835658 0.1919378
# [6,]    1 0.8983897 0.03910006 0.0822944

b <- rnorm(J + 1, 0, 10)
sigma <- runif(1)
y <- c(X %*% b) + rnorm(n, 0, sigma)

并拟合模型：

library(R2jags)
out <- jags(list(y=y, X=X, n=n, J=ncol(X)), NULL, c('beta', 'sd'), M)

现在将估计的系数和标准差与其真实值进行比较：

out$BUGSoutput$summary[, '50%']
#     beta[1]      beta[2]      beta[3]      beta[4]     deviance           sd 
#   8.5783315   -9.3849277    8.9632598   -9.4242272 2196.1535645    0.7264754 

b # True betas
# [1]  8.500435 -9.253130  8.935812 -9.410097

sigma # True sd
# [1] 0.70504
于 2014-11-17T08:26:18.300 回答