在R中使用Arm函数进行Gibbs采样

问题描述

我是gibbs采样的新手，我编写了以下代码，但由于某些原因，输出仍然为空。你能告诉我我的错误在哪里，它是基于y =（y1，y2，y3，y4）的，其中y.mis = y4是缺失的，给定（y1，y2，y3），我们应该生成y4。

library(armspp)

log.pdf.yi <- function(y.mis,y.obs,Trt = c(0,1,1),Time = c(0,3,6,9),r = c(1,0),beta = c(5,-0.2,-0.1),sigmasqu = 0.701,sigmasqe = 2.429,phi = c(-1,0.5,0.25,-0.25)) {

  yi <- c(y.obs,y.mis)
  Xi <- cbind(x0 = 1,Trt,Time)
  ni <- nrow(Xi)
  D <- matrix(sigmasqu,ni,ni)
  Vi <- diag(rep(sigmasqe,ni)) + D
  Vi.inv <- solve(Vi)
  ei <- yi - c(Xi %*% beta)
  log.fyi <- -0.5 * c(t(ei) %*% Vi.inv %*% ei)

  z1 <- c(1,yi[1],Trt[1])
  z2 <- c(1,yi[2],Trt[2])
  z3 <- c(1,yi[3],Trt[3])
  z4 <- c(1,yi[4],Trt[4])

  zi <- rbind(z1,z2,z3,z4)
  etai <- zi %*% phi
  pi <- exp(etai)/(1 + exp(etai))
  log.fr <- sum(log(dbinom(r,size = 1,prob = pi)))

  log.fyi + log.fr
}

samples <- function(n_sam = 5,k = 10,n = 4,beta=c(5,sigmasqu=0.701,sigmasqe=2.429,-0.25)) {
  id <- rep(c(1:k),rep(n,k))
  u.id <- unique(id)
  trt  <- rep(c(0,n*k/2)
  time <- rep(c(0,k)
  u1 <- rnorm(k,mean = 0,sd = sqrt(sigmasqu))
  u  <- rep(u1,k))
  e <- rnorm(k*n,sd = sqrt(sigmasqe))
  x_beta <- cbind(1,trt,time)
  Y = as.vector(x_beta %*% beta + u + e)

  g_mat <- matrix(0,ncol = k,nrow = n_sam)
  for (i in 1:k) {
    Yi <- Y[id == u.id[i]]
    ni <- length(Yi)
    yi.obs <- as.vector(Yi[-ni])
    print(yi.obs)
    y0 <- arms(n_samples = n_sam,log.pdf.yi,-100,100,arguments = list(y.obs = yi.obs))
    g_mat[,i] <- y0
  }
  g_mat
}


samples()

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）