not able to set different deterministic starting values for each chain

[Rebecca Taylor]

Hi all.

I am new to NIMBLE, and I am trying unsuccessfully to set different deterministic initial values for each chain. The help file for runMCMC states that the inits argument can be "a list of length nchains, each element being a named list of initial values which be used for one MCMC chain." However, NIMBLE appears to only use the first list of initial values I give it. In the following example, I ask for two chains, in each of two different runs. In the first run, I ask for two different sets of initial values:

inits = list(inits1=I1, inits2=I2).

In the second run I ask for  the same set of initial values for each chain:

inits = list(inits1=I1, inits2=I1).

I get identical chains in both runs. 

I will reprint the code below from the run where I ask for different initial values, and I will reprint the output I get regardless of the run.

set.seed(10120)
myData <- rgamma(1000, shape = 0.4, rate = 0.8)
myCode <- nimbleCode({
  a ~ dunif(0, 100)
  b ~ dnorm(0, 100)    
  for (i in 1:length_y) {
    y[i] ~ dgamma(shape = a, rate = b)
  }
})
I1 <- list(a=0.3,b=0.9)
I2 <- list(a=0.5, b=0.7)
myModel <- nimbleModel(code = myCode,
                       data = list(y = myData),
                       constants = list(length_y = 1000),
                       inits = list(inits1=I1, inits2=I2))
CmyModel <- compileNimble(myModel)  
myMCMC <- buildMCMC(CmyModel)
CmyMCMC <- compileNimble(myMCMC)
results <- runMCMC(CmyMCMC, niter = 10, setSeed = c(1,1), nchains=2)
 
cbind(as.numeric(results[[1]][,"a"]),as.numeric(results[[2]][,"a"]),as.numeric(results[[1]][,"b"]),as.numeric(results[[2]][,"b"]))

#       a1   a2      b1       b2
# [1,]  0.3  0.3 0.9000000 0.6052796
# [2,]  0.3  0.3 0.9000000 0.6052796
# [3,]  0.3  0.3 0.6052796 0.6052796
# [4,]  0.3  0.3 0.6052796 0.6052796
# [5,]  0.3  0.3 0.6052796 0.6052796
# [6,]  0.3  0.3 0.6052796 0.6052796
# [7,]  0.3  0.3 0.6052796 0.6052796
# [8,]  0.3  0.3 0.6052796 0.6052796
# [9,]  0.3  0.3 0.6052796 0.6052796
# [10,] 0.3  0.3 0.6052796 0.6052796

If anyone can offer me some advice, I'd very much appreciate it. Once I figure out how to do this with sequential chains, I hope to do it with parallel chains. (I can get everything else to work in parallel.)

Thank you!!

Rebecca

[Perry de Valpine]

Hi Rebecca,

Thanks for posting.

One thing I see is that the documentation about an inits list is for runMCMC, but you are providing an inits list to nimbleModel and not to runMCMC.

Can you try providing the inits list to runMCMC?

Fingers crossed.

Perry

--
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/55a137f7-4f77-4fb5-90c5-70510eb08644n%40googlegroups.com.

[Taylor, Rebecca L]

Perry, 

Thank you!!! (I should have seen that.) When I supply inits to runMCMC, it works as it is supposed to 🙂.

This does still leave me with the question about how to specify different initial values when I parallelize the chains. When I supply inits to runMCMC via a parLapply function, it still only wants one set of initial values. It's like I need to be able to write a function that says, "if you are worker #1, use inits1; if you are worker#2, use inits2..." However, given my rudimentary knowledge of both NIMBLE and package "parallel", I don't know how to proceed from here.

The following code works in every other respect (and yes, it is nearly identical to an example one of you wrote).

library(parallel)

this_cluster <- makeCluster(2)

set.seed(10120) 

myData <- rgamma(1000, shape = 0.4, rate = 0.8)

I1 <- list(a=0.3,b=0.9)

I2 <- list(a=0.5, b=0.7)

run_MCMC_allcode <- function(seed, data, inits) {

  library(nimble) 

  myCode <- nimbleCode({

    a ~ dunif(0, 100)

    b ~ dnorm(0, 100)   

    for (i in 1:length_y) {

      y[i] ~ dgamma(shape = a, rate = b)

    }

  }) 

  myModel <- nimbleModel(code = myCode,

                         data = list(y = data),

                         constants = list(length_y = 1000)) 

  CmyModel <- compileNimble(myModel) 

  myMCMC <- buildMCMC(CmyModel)

  CmyMCMC <- compileNimble(myMCMC)

 

  results <- runMCMC(CmyMCMC, niter = 10000, setSeed = seed, inits = inits)

  return(results)

}

chain_output <- parLapply(cl = this_cluster, X = 1:2, fun = run_MCMC_allcode, data = myData, inits = I1)

Rebecca

Rebecca Taylor

Research Statistician

Alaska Science Center

U. S. Geological Survey

(907) 786-7004

Pronouns: She/her

Please note I sometimes send and answer email outside normal work hours. 

This helps me with my schedule but is not intended to pressure others to do the same. 

---

From: nimble...@googlegroups.com <nimble...@googlegroups.com> on behalf of Perry de Valpine <pdeva...@berkeley.edu>
Sent: Saturday, July 9, 2022 4:12 PM
To: Taylor, Rebecca L <rebecc...@usgs.gov>
Cc: nimble-users <nimble...@googlegroups.com>
Subject: [EXTERNAL] Re: not able to set different deterministic starting values for each chain

 

 

This email has been received from outside of DOI - Use caution before clicking on links, opening attachments, or responding.

You received this message because you are subscribed to a topic in the Google Groups "nimble-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/nimble-users/h098rf08Cu4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CABeTmqjgskOqCOACsgZXMH4t4AdY7-tk6-sdbFaseZbupNF1Eg%40mail.gmail.com.

[Perry de Valpine]

Hi Rebecca,

There's an example at r-nimble.org/examples on parallelizing nimble.  If you haven't been through that yet, please see if it sets you on course and let us know otherwise.  The main idea is that each thread needs its own model object and subsequent building steps.

Best wishes,

Perry

[Taylor, Rebecca L]

Hi Perry,

Thanks--yes--I've been through that extensively, and I have been able to modify that simple example in all other respects for my more complicated model. But I'm stuck on the initial values part of it.

Rebecca

Rebecca Taylor

Research Statistician

Alaska Science Center

U. S. Geological Survey

(907) 786-7004

Pronouns: She/her

Please note I sometimes send and answer email outside normal work hours. 

This helps me with my schedule but is not intended to pressure others to do the same. 

---

From: Perry de Valpine <pdeva...@berkeley.edu>
Sent: Saturday, July 9, 2022 5:12 PM

To: Taylor, Rebecca L <rebecc...@usgs.gov>
Cc: nimble-users <nimble...@googlegroups.com>

Subject: Re: [EXTERNAL] Re: not able to set different deterministic starting values for each chain

 

[Perry de Valpine]

Hi Rebecca,

My apologies -- I was hasty in wanting to get you a quick answer but didn't quite see that you were clearly working from our example.  Oops.

I think there are two ways you could do this.  The first is that you could make the X argument in parLapply be a list of lists, and then each inner list will be passed to the function.  Here is an example of how that looks (not using nimble).

# Modified from the examples in help(parLapply)

cl <- makeCluster(getOption("cl.cores", 2))

parLapply(cl, X = list(list(a = 1, b = 2), list(a = 2, b = 3)), function(ab) ab$a+ ab$b )

With this trick, each inner list could contain arbitrarily many objects needed inside run_MCMC_allcode, including the inits object.  This could then be different for each thread.

Another option would be to use clusterMap, which like mapply iterates over multiple inputs together.

I'm not an ace with parallelization, so if these don't give a path forward, someone else can jump in and help. 

I hope it works!

Perry

[Taylor, Rebecca L]

Thank you, Perry--no apologies needed! I will see what I can do with those ideas, and then I will either return to the list for more help or post my success to the list so others who may be interested can see the resolution. It will take me a few days to return to this (sigh).

Rebecca

Rebecca Taylor

Research Statistician

Alaska Science Center

U. S. Geological Survey

(907) 786-7004

Pronouns: She/her

Please note I sometimes send and answer email outside normal work hours. 

This helps me with my schedule but is not intended to pressure others to do the same. 

---

From: Perry de Valpine <pdeva...@berkeley.edu>

Sent: Saturday, July 9, 2022 5:35 PM

[Taylor, Rebecca L]

Hi Perry,

I'm still stuck. I've tried your suggestions and I've read some tutorials on parallelization in R, and I'm just not getting it. Might Daniel or Chris be willing to jump in?

Thanks!

Rebecca

Rebecca Taylor

Research Statistician

Alaska Science Center

U. S. Geological Survey

(907) 786-7004

Pronouns: She/her

Please note I sometimes send and answer email outside normal work hours. 

This helps me with my schedule but is not intended to pressure others to do the same. 

---

From: Perry de Valpine <pdeva...@berkeley.edu>

Sent: Saturday, July 9, 2022 5:35 PM

[Perry de Valpine]

Hi Rebecca,

Here are a couple of more worked out options for you.  In both options, you would set up a list of lists of inputs you want, such as inits and a seed for random number generation.  Then you could either export that to cluster nodes, which makes it available to each node, and access values based on the simple integer value of "X" passed via parLapply.  Or, you would provide the list as the X argument to parLapply, so an element of the list would be the X argument in a node, and unpack it from there.  The options are very similar and I hope just help in seeing how things work.

Note that I have also modified the output of "run_MCMC_allcode" to return a list with the results and also the seed value, for the purpose of inspection to see that each seed value was used in one cluster node.

-Perry

library(parallel)
this_cluster <- makeCluster(2) # Use only 2 for simplicity

set.seed(10120)
# Simulate some data

myData <- rgamma(1000, shape = 0.4, rate = 0.8)

# Option 1: make global lists and use the X input from parLapply
# to look up elements.
run_MCMC_allcode <- function(X, data) {

  library(nimble)

  myCode <- nimbleCode({
    a ~ dunif(0, 100)
    b ~ dnorm(0, 100)

    for (i in 1:length_y) {
      y[i] ~ dgamma(shape = a, rate = b)
    }
  })

  inputs <- inputsList[[X]]
  inits <- inputs$inits
  seed <- inputs$seed

  myModel <- nimbleModel(code = myCode,
                          data = list(y = data),

                          constants = list(length_y = 1000),
                          inits = inits)

  CmyModel <- compileNimble(myModel)

  myMCMC <- buildMCMC(CmyModel)
  CmyMCMC <- compileNimble(myMCMC)

  results <- runMCMC(CmyMCMC, niter = 10000, setSeed = seed)

  return(list(results = results, seed = seed))
}

inputsList <- list(
  list(seed = 123, # The seeds are arbitrary numbers
       inits = list(a = 0.5, b = 0.5)),
  list(seed = 234,
       inits = list(a = 0.3, b = 0.3)))
clusterExport(this_cluster, c("inputsList"))

chain_output <- parLapply(cl = this_cluster, X = 1:2,
                          fun = run_MCMC_allcode,

                          data = myData)
plot(chain_output[[1]]$results[, 'a'])

# Option 2: Make the X input to parLapply be a list
# including elements such as seed and inits
#
run_MCMC_allcode <- function(X, data) {

  library(nimble)

  myCode <- nimbleCode({
    a ~ dunif(0, 100)
    b ~ dnorm(0, 100)

    for (i in 1:length_y) {
      y[i] ~ dgamma(shape = a, rate = b)
    }
  })

  inits <- X$inits

  myModel <- nimbleModel(code = myCode,
                          data = list(y = data),

                          constants = list(length_y = 1000),
                          inits = inits)

  CmyModel <- compileNimble(myModel)

  myMCMC <- buildMCMC(CmyModel)
  CmyMCMC <- compileNimble(myMCMC)

  seed <- X$seed
  results <- runMCMC(CmyMCMC, niter = 10000, setSeed = seed)

  return(list(results = results, seed = seed))
}

# Xlist is the same as inputsList above.
# I'm  just naming it Xlist to make clear it will be
# the value passed as the X argument.
Xlist <- list(
  list(seed = 123,
       inits = list(a = 0.5, b = 0.5)),
  list(seed = 234,
       inits = list(a = 0.3, b = 0.3)))
chain_output <- parLapply(cl = this_cluster, X = Xlist,
                          fun = run_MCMC_allcode,
                          data = myData)

plot(chain_output[[1]]$results[, 'a'])

stopCluster(this_cluster)

[Taylor, Rebecca L]

Thank you Perry!!  I'll dig into this.

Rebecca

Rebecca Taylor

Research Statistician

Alaska Science Center

U. S. Geological Survey

(907) 786-7004

Pronouns: She/her

Please note I sometimes send and answer email outside normal work hours. 

This helps me with my schedule but is not intended to pressure others to do the same. 

---

From: Perry de Valpine <pdeva...@berkeley.edu>

Sent: Thursday, July 14, 2022 10:16 AM

[Taylor, Rebecca L]

Perry,

Thank you SOOO much! I am running at full speed now.

Rebecca

Rebecca Taylor

Research Statistician

Alaska Science Center

U. S. Geological Survey

(907) 786-7004

Pronouns: She/her

Please note I sometimes send and answer email outside normal work hours. 

This helps me with my schedule but is not intended to pressure others to do the same. 

---

From: Perry de Valpine <pdeva...@berkeley.edu>

Sent: Thursday, July 14, 2022 10:16 AM