0

I have set up a Metropolis-Hastings algorithm, and now I am trying to run the algorithm using parallel computing. I have set up a single-chain function

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)
}

, which returns an p x iter matrix, with p parameters and iter iterations. 

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

UPDATE: When I tried to simplify the problem the code suddenly seemed to work. Even though I purposely tried to make errors, the code ran perfectly and the results were as wanted. So for others with similar problems unfortunately I cannot give an answer.

A follow-up question: My initial purpose was to built up an entire function, such that

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

but the foreach-loop does not seem to read all the required functions like:

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

Can anybody answer how to implement custom functions inside a foreach loop? Am I to input it as MCmulti <- function(FUN,...) FUN() and call MCmulti(MCsingle,...) ?

svick
  • 236,525
  • 50
  • 385
  • 514
Troels
  • 1
  • 1
  • You probably need a `.combine` argument. – tchakravarty May 12 '15 at 15:52
  • I have tried that. First of all I just need an output. – Troels May 12 '15 at 18:24
  • Show us `MCSingle()` or reproduce with simple `MCSingle()`. – tchakravarty May 12 '15 at 18:25
  • The MCsingle() function is included now. – Troels May 12 '15 at 18:56
  • Still not reproducible: `getinits()` is missing. – tchakravarty May 12 '15 at 19:31
  • Code should be reproducible now. – Troels May 12 '15 at 20:20
  • Your code is not reproducible at all. Package loading statements are missing, `inits` not defined, `r` not defined. Please test your reproducible example in a new R environment to understand what is missing. – tchakravarty May 13 '15 at 11:28
  • I appologize for the lack of reproduction, since it is my first strike at `foreach` I thought I lacked some kind of general knowledge, such that the problem could be solved from a quick view at it. In the effort of reproducing the code, however, I seemed to solve the problems, as I was unable to reproduce the error, so I build it up to work in my own setting too. – Troels May 14 '15 at 14:04
  • That is usually the case -- keeping a clean workspace is usually a good idea. Good to hear that you solved the problem. – tchakravarty May 15 '15 at 05:49

0 Answers0