如何使用来自 rjags / JAGS 的估计值来预测值

Posted 2023-02-19

【中文标题】如何使用来自 rjags / JAGS 的估计值来预测值

【英文标题】：How to predict values using estimates from rjags / JAGS

【发布时间】：2016-07-03 06:46:49

【问题描述】：

在建立模型并使用 Gibbs Sampling 对其进行训练后，我得到了所有隐藏值预测的结果：

jags <- jags.model('example.bug',
               data = data,
               n.chains = 4,
               n.adapt = 100)

update(jags, 1000)

samples <- jags.samples(jags,
         c('r','alpha','alpha_i','alpha_u','u','i'),
         1000)

r 是一个评级列表，其中一些被保留用于模型的预测。假设我可以使用r[test] 获得它们，其中test 是一个整数列表，表示保留评级的索引。但是当我试图让他们使用这种方式时：

summary(samples$r, mean)[test]

我刚刚得到这个：

$drop.dims
iteration     chain 
 1000         4 

你能告诉我如何得到期望值吗？提前谢谢！

【问题讨论】：

为什么不mean(r[test])？

@effel 不，r 应该是评级列表，其中一些是 NA。它在模型中使用。在这种情况下，JAGS 生成的预测值应该是samples。

您可能会在这里找到有用的答案：***.com/questions/33662987/…

【参考方案1】：

画samples

如果没有您的数据或模型，我将使用简单示例 here 进行演示，经过修改后，jags 可以监控预测结果。

library(rjags)

# simulate some data    
N <- 1000
x <- 1:N
epsilon <- rnorm(N, 0, 1)
y <- x + epsilon

# define a jags model
writeLines("
  model 
    for (i in 1:N)
      y[i] ~ dnorm(y.hat[i], tau)
      y.hat[i] <- a + b * x[i]
    
    a ~ dnorm(0, .0001)
    b ~ dnorm(0, .0001)
    tau <- pow(sigma, -2)
    sigma ~ dunif(0, 100)
  
", con = "example2_mod.jags")

# create a jags model object
jags <- jags.model("example2_mod.jags",
                   data = list('x' = x,
                               'y' = y,
                               'N' = N),
                   n.chains = 4,
                   n.adapt = 100)

# burn-in
update(jags, 1000)

# drawing samples gives mcarrays
samples <- jags.samples(jags, c('a', 'b'), 1000)
str(samples)
# List of 2
#  $ a: mcarray [1, 1:1000, 1:4] -0.0616 -0.0927 -0.0528 -0.0844 -0.06 ...
#   ..- attr(*, "varname")= chr "a"
#  $ b: mcarray [1, 1:1000, 1:4] 1 1 1 1 1 ...
#   ..- attr(*, "varname")= chr "b"
# NULL

提取预测

我们的结果samples 是mcarray 对象的列表，其维度为 1 x 迭代 x 链。此时您真的很想运行诊断，但我们将跳转到从后验中总结样本以进行预测。一种方法是对链和迭代取平均值。

# extract posterior means from the mcarray object by marginalizing over 
# chains and iterations (alternative: posterior modes)
posterior_means <- apply(samples$y.hat, 1, mean)
head(posterior_means)
# [1] 0.9201342 1.9202996 2.9204649 3.9206302 4.9207956 5.9209609

# reasonable?
head(predict(lm(y ~ x)))
#         1         2         3         4         5         6 
# 0.9242663 1.9244255 2.9245847 3.9247439 4.9249031 5.9250622 

样本外预测

或者，您可以通过以下方式进行样本外预测。我将重用我们现有的特征向量x，但这可能是测试数据。

# extract posterior means from the mcarray object by marginalizing over chains and iterations (alternative: posterior modes)
posterior_means <- lapply(samples, apply, 1, "mean")
str(posterior_means)
# List of 3
#  $ a    : num -0.08
#  $ b    : num 1
#  $ y.hat: num [1:1000] 0.92 1.92 2.92 3.92 4.92 ...
# NULL


# create a model matrix from x
X <- cbind(1, x)
head(X)
#        x
# [1,] 1 1
# [2,] 1 2
# [3,] 1 3
# [4,] 1 4
# [5,] 1 5
# [6,] 1 6

# take our posterior means 
B <- as.matrix(unlist(posterior_means[c("a", "b")]))
#          [,1]
# a -0.07530888
# b  1.00015874

根据模型，预测结果是 a + b * x[i]，正如我们在 jags 中所写的那样。

# predicted outcomes are the product of our model matrix and estimates
y_hat <- X %*% B
head(y_hat)
#           [,1]
# [1,] 0.9248499
# [2,] 1.9250086
# [3,] 2.9251673
# [4,] 3.9253261
# [5,] 4.9254848
# [6,] 5.9256436

【讨论】：

感谢您的回答？但是如何extract e.g. posterior means from from your mcarrayobject 和X'B 是什么意思？

你能提供data吗？

更新演示。