R 깁스 샘플러 for 베이지안 회귀

Question

R의 베이지안 회귀 모델 용 깁스 샘플러를 코딩하려고하는데 코드 실행에 문제가 있습니다. sigma.update 함수에서 베타 버전이 계속 진행되고있는 것 같습니다. 나는이 코드를 실행하면 나는 말한다 오류 얻을 "의 X %의 *의 %의 베타 오류 : 비 정합 인수"여기처럼 내 코드는 모습입니다 : 그것은 NaN을 반환

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가 정확하지 기능 . 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