checkDistributionFunctions: density function for is not available for parallel in R

[Huong Nguyen]

Hi all,

I created the new distribution for piMOM (product of inverse moments) and run the code successfully individually. However, when I put the code into the big function so that I can use parallel (foreach) to run in multiple cores, it sent out errors that: "checkDistributionFunctions: density function for dpiMOM is not available. It must be a nimbleFunction (with no setup code)". Does anyone have any suggestions of what's wrong? Thank you so much in advance.

final <- c()
output <- c()
output2 <- c()

library(nimble)
library(rjags)
library(knitr)


myFunction <- function (k){

  # piMOM density
  dpiMOM <- nimbleFunction(
    run = function(x = double(0), t = double(0), r = double(0),
                   log = integer(0, default = 0)) {
      returnType(double(0))
      logProb <- (r/2)*log(t)-lgamma((r/2))+log(abs(x)^(-r-1))-t/(x^2)
      if(log) return(logProb)
      else return(exp(logProb))
    })
  

# rpiMOM for simulation
  rpiMOM <- nimbleFunction(
    run = function(n=integer(0), t = double(0), r = double(0)){
      returnType(double(0))
      if (n!=1) print("rpi_MOM only allows n = 1; using n=1.")
      dev <- (t^(r/2))/gamma(r/2)*exp(1)^(-t)
      return(dev)
    })
 
  set.seed(k)
  n <- 60 # number of observations
  p <- 1000 # number of covariates
  X <- matrix(rnorm(n*p),n,p)
  sigma <- 1
  beta <- rep(0,p)
  beta[1:10] <- (-1)^(1:10)
  X<- scale(X)
  Y <- rnorm(n, X %*% beta, sigma)

  m_piMOM <- nimbleCode({
    for(i in 1:n){
      Y[i]   ~ dnorm(alpha+inprod(X[i,1:p],beta[1:p]),sigma_e)
    }
    for(j in 1:p){
      beta[j] <- gamma[j]*delta[j]
      gamma[j] ~ dbern(q)
      delta[j] ~ dpiMOM(t,r)
    }
   
    alpha ~ dnorm(0,0.01)
    sigma_e ~ dgamma(0.1,0.1)
   
    q ~ dbeta(1,1)
    t <-5
    r <-2
  })
 
  data   <- list(Y=Y,X=X)
 
  samples <- nimbleMCMC(m_piMOM,data = data, constants = list(p=p, n=n),
                        monitors = c("beta", "gamma"), thin = 5,
                        niter = 5000, nburnin = 1000, nchains = 1)
 
  library(coda)
  coda.samples <- as.mcmc(samples)
 
  burn <- 1000
  # log-score for non-local prior
  p_piMOM = 0
  logloss_piMOM = 0
  gamma_set_piMOM <- coda.samples[,1001:2000]
 
  for (j in 1:p){
    total=0
    total = total + sum(gamma_set_piMOM[,j])
    p_piMOM[j] = total/burn
    logloss_piMOM[j] = -(sum(gamma_set_piMOM[,j])*log(p_piMOM[j])+sum(1-gamma_set_piMOM[,j])*log(1-p_piMOM[j]))
    logloss_piMOM[j] = logloss_piMOM[j]/p
  }
 
  # MSE for non-local betas
  beta_set_piMOM <- coda.samples[,1:1000]
  mse_piMOM = 0
  for (j in 1:p){
    total = 0
    total = total + (beta_set_piMOM[,j]-beta)^2
    mse_piMOM[j] = sum(total)/(burn)
  }
 
  output <- c(mean(logloss_piMOM, na.rm=TRUE), mean(mse_piMOM, na.rm=TRUE))
  output2 <- c(logloss_piMOM, mse_piMOM)
  final <- list(output,output2)
  return(final)
}

library(doFuture)
library(doRNG)
registerDoFuture()
plan(multisession)

id <- sample(1000,20)
result0419pi <- foreach(i = id, .combine='rbind', .packages="nimble") %dopar% {
  myFunction(i)
}

[Chris Paciorek]

hi Huong,

Please see the bottom of this nimble example that shows how to deal with scoping so that the user-defined distribution can be found when used in parallel.

-chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/785ab043-8b73-4e4f-b971-0e6fc963b706n%40googlegroups.com.