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 個のパラメーターと反復反復を含むp x反復行列を返します。
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-loop は、次のような必要なすべての関数を読み取っていないようです:
Error in FUN() : task 1 failed - "could not find function "geninits""
foreach
ループ内にカスタム関数を実装する方法を誰かが答えることができますか? と入力してMCmulti <- function(FUN,...) FUN()
呼び出すのMCmulti(MCsingle,...)
ですか?