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

Question

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

4

我在尝试用R编写贝叶斯回归模型的Gibbs采样器，但是运行代码时遇到了问题。sigma.update函数中的beta似乎出了点问题。当我运行代码时，会出现一个错误，提示“Error in x %*% beta：non-conformable arguments”。以下是我的代码：

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)

- user11235813

1

欢迎来到CV！我在您的代码中添加了一些空格，以使其更易读。此外，考虑使用r标签，因为这是与R相关的问题。 - Tim

您的错误提示表明您可能忘记转置某些内容了 - 很抱歉我现在没有时间来审查代码。 - Tim

1

beta.draws 应该是一个两列矩阵，beta.update 应该生成两个值，使用 rnorm(2,...)。这将解决错误，但你仍然应该检查方程和结果是否正确。
小贴士：对于形如 X'X 或 XX' 的矩阵乘积，使用函数 crossprod 或 tcrossprod。

- javlacalle

@javlacalle 我修复了你提到的问题，并且成功解决了我之前遇到的错误，但是现在在 sigma.draws 和 beta.draws 中除了第一个条目外，输出都是 NaN。有什么其他问题可能出错了吗？顺便再次感谢你。 - stochasticcrap

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- javlacalle · Accepted Answer

函数beta.update不正确，它返回NaN。您在参数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