用于贝叶斯回归的 R Gibbs 采样器

Question

我正在尝试为 R 中的贝叶斯回归模型编写 Gibbs 采样器，但在运行我的代码时遇到了问题。 sigma.update 函数中的 beta 似乎发生了一些事情。当我运行代码时，我收到一条错误消息，上面写着“x %*% beta 中的错误：不符合要求的参数”这是我的代码的样子：

x0 <- rep(1, 1000)
x1 <- rnorm(1000, 5, 7)
x <- cbind(x0, x1)
true_error <- rnorm(1000, 0, 2)
true_beta <- c(1.1, -8.2)
y <- x %*% true_beta + true_error

beta0 <- c(1, 1)
sigma0 <- 1  
a <- b <- 1
burnin <- 0
thin <- 1
n <- 100

gibbs <- function(n.sims, beta.start, a, b,
                  y, x, burnin, thin) {
   beta.draws <- matrix(NA, nrow=n.sims, ncol=1)
   sigma.draws<- c()
   beta.cur <- beta.start
   sigma.update <- function(a,b, beta, y, x) {
        1 / rgamma(1, a + ((length(x)) / 2),
                   b + (1 / 2) %*% (t(y - x %*% beta) %*% (y - x %*% beta)))
     }
   beta.update <- function(x, y, sigma) {
        rnorm(1, (solve(t(x) %*% x) %*% t(x) %*% y),
              sigma^2 * (solve(t(x) %*%x)))
     }
   for (i in 1:n.sims) {
     sigma.cur <- sigma.update(a, b, beta.cur, y, x)
     beta.cur <- beta.update(x, y, sigma.cur)
     if (i > burnin & (i - burnin) %% thin == 0) {
       sigma.draws[(i - burnin) / thin ] <- sigma.cur
       beta.draws[(i - burnin) / thin,] <- beta.cur
       }
     }
   return (list(sigma.draws, beta.draws) )
   }

gibbs(n, beta0, a, b, y, x, burnin, thin)

Answer 1

函数beta.update不正确，它返回NaN。您在传递给rnorm 的参数sd 中定义一个矩阵，该参数中需要一个向量。我认为您正在尝试做的事情可以通过这种方式完成：

beta.update <- function(x, y, sigma) {
  rn <- rnorm(n=2, mean=0, sd=sigma)
  xtxinv <- solve(crossprod(x))
  as.vector(xtxinv %*% crossprod(x, y)) + xtxinv %*% rn
}

请注意，您正在计算一些在所有迭代中都固定的元素。例如，您可以定义一次t(x) %*% x，然后将此元素作为参数传递给其他函数。通过这种方式，您可以避免在每次迭代时都执行这些操作，从而节省一些计算，可能还节省一些时间。

编辑

根据您的代码，这就是我所做的：

x0 <- rep(1, 1000)
x1 <- rnorm(1000, 5, 7)
x <- cbind(x0, x1)
true_error <- rnorm(1000, 0, 2)
true_beta <- c(1.1, -8.2)
y <- x %*% true_beta + true_error

beta0 <- c(1, 1)
sigma0 <- 1  
a <- b <- 1
burnin <- 0
thin <- 1
n <- 100

gibbs <- function(n.sims, beta.start, a, b, y, x, burnin, thin) 
{
  beta.draws <- matrix(NA, nrow=n.sims, ncol=2)
  sigma.draws<- c()
  beta.cur <- beta.start
  sigma.update <- function(a,b, beta, y, x) {
    1 / rgamma(1, a + ((length(x)) / 2),
    b + (1 / 2) %*% (t(y - x %*% beta) %*% (y - x %*% beta)))
  }
  beta.update <- function(x, y, sigma) {
    rn <- rnorm(n=2, mean=0, sd=sigma)
    xtxinv <- solve(crossprod(x))
    as.vector(xtxinv %*% crossprod(x, y)) + xtxinv %*% rn
  }
  for (i in 1:n.sims) {
    sigma.cur <- sigma.update(a, b, beta.cur, y, x)
     beta.cur <- beta.update(x, y, sigma.cur)
     if (i > burnin & (i - burnin) %% thin == 0) {
       sigma.draws[(i - burnin) / thin ] <- sigma.cur
       beta.draws[(i - burnin) / thin,] <- beta.cur
     }
  }
  return (list(sigma.draws, beta.draws) )
}

这就是我得到的：

set.seed(123)
res <- gibbs(n, beta0, a, b, y, x, burnin, thin)
head(res[[1]])
# [1] 3015.256257   13.632748    1.950697    1.861225    1.928381    1.884090
tail(res[[1]])
# [1] 1.887497 1.915900 1.984031 2.010798 1.888575 1.994850
head(res[[2]])
#          [,1]      [,2]
# [1,] 7.135294 -8.697288
# [2,] 1.040720 -8.193057
# [3,] 1.047058 -8.193531
# [4,] 1.043769 -8.193183
# [5,] 1.043766 -8.193279
# [6,] 1.045247 -8.193356
tail(res[[2]])
#            [,1]      [,2]
# [95,]  1.048501 -8.193550
# [96,]  1.037859 -8.192848
# [97,]  1.045809 -8.193377
# [98,]  1.045611 -8.193374
# [99,]  1.038800 -8.192880
# [100,] 1.047063 -8.193479