R Gibbs Sampler for Bayesian Regression

Question

Я пытаюсь закодировать пробоотборник Gibbs для байесовской модели регрессии в R, и у меня возникают проблемы с запуском моего кода. Кажется, что что-то происходит с бета-версией в функции sigma.update. Когда я запускаю код я получаю сообщение об ошибке «Ошибка при й% *% бете: без облегать аргументы» Вот то, что мой код выглядит следующим образом:

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 , Вы определяете матрицу в аргументе sd, который передается в rnorm, в этом аргументе ожидается вектор. Я думаю, что вы пытаетесь сделать, может быть сделано таким образом:

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