0

我已经建立了 Metropolis-Hastings 算法，现在我正在尝试使用并行计算运行该算法。我已经设置了单链功能

library(parallel)
library(foreach)
library(mvtnorm)
library(doParallel)

n<-100
mX <- 1:n
vY <- rnorm(n)
chains <- 4
iter <- n
p <- 2

#Loglikelihood
post <- function(y, theta)  dmvnorm(t(y), rep(0,length(y)), theta[1]*exp(-         abs(matrix(rep(mX,n),n) - matrix(rep(mX,each=n),n))/theta[2]),log=TRUE)

geninits <- function()  list(theta = runif(p, 0, 1))

dist <- 0.01
jump <- function(x, dist)  exp(log(x) + rmvnorm(1,rep(0,p),diag(rep(dist,p))))

MCsingle <- function(){ # This is part of a larger function, so no input are needed
 inits <- geninits()
 theta.post <- matrix(NA,nrow=p,ncol=iter)
 for (i in 1:p) theta.post[i,1] <- inits$theta[i]
 for (t in 2:iter){
  theta_star <- c(jump(theta.post[, t-1],dist))
  pstar <- post(vY, theta = theta_star) # post is the loglikelihood using dmvnorm.
  pprev <- post(vY, theta = theta.post[,t-1])
  r <- min(exp(pstar - pprev) , 1)
  accept <- rbinom(1, 1, prob = r)
  if (accept == 1){
   theta.post[, t] <- theta_star
   } else {
    theta.post[, t] <- theta.post[, t-1]
   }
 }
return(theta.post)
}

，它返回一个p x iter矩阵，带有 p 个参数和iter迭代。

cl<-makeCluster(4)
registerDoParallel(cl)
posterior <- foreach(c = 1:chains) %dopar% {
MCsingle()  }

更新：当我试图简化问题时，代码突然似乎起作用了。尽管我故意尝试出错，但代码运行完美，结果符合预期。因此，不幸的是，对于其他有类似问题的人，我无法给出答案。

一个后续问题：我最初的目的是建立一个完整的功能，这样

MCmulti <- function(mX,vY,iter,chains){
    posterior <- foreach(c = 1:chains) %dopar% {
    MCsingle()  }
    return(posterior)
   }

但是 foreach 循环似乎没有读取所有必需的功能，例如：

Error in FUN() : task 1 failed - "could not find function "geninits""

谁能回答如何在foreach循环中实现自定义函数？我要输入它MCmulti <- function(FUN,...) FUN()并调用MCmulti(MCsingle,...)吗？

4

0 回答 0