RJMCMC model won't compile in latest version of nimble 0.12.1

[Daniel Eacker]

Dear nimble users/developers,

I received some help getting the sampler_RJ_indicator function to work for a problem where I wanted to select among both categorical and continuous variables and have the algorithm sample the categorical variable as one variable (although it is represented by multiple indicator variables in the model). 

This seemed to work fine in the last version of nimble that I have on my PC, but when I updated the version and tried to run it on my Mac it fails to compile and reports the error that is related to a modified version of the sampler_RJ_indicator that I created with some help (thanks Sally!). 

I get this error related to the custom sampler (my_sampler_RJ_indicator):

mcmc5 = compileNimble(mcmc, project = RJexampleModel, resetFunctions = TRUE, showCompilerOutput = TRUE)

Compiling

  [Note] This may take a minute.

  [Note] On some systems there may be some compiler warnings that can be safely ignored.

Error : $ operator is invalid for atomic vectors

Error: There is some problem at the Eigen processing step for this code:

I'm not sure why the code would run fine in a previous version and now it is throwing this error. The inner workings of nimble are a bit over my head, so I was hoping this would be more obvious to someone involved with the development.

Here is the example code that I was trying to run the function with:

## example MCMC rj

# my_sampler_RJ_indicator #######################################

my_sampler_RJ_indicator <- nimbleFunction(

  name = 'my_sampler_RJ_indicator',

  contains = sampler_BASE,

  setup = function(model, mvSaved, target, control) {

    ## note: target is the indicator variable,

    ## control$targetNode is the variable conditionally in the model

    ## control list extraction

    coefNode      <- model$expandNodeNames(control$targetNode, returnScalarComponents = TRUE)

    coefNodes     <- control$targetNode

    proposalScale <- control$scale

    proposalMean  <- control$mean

    ## node list generation

    calcNodes <- model$getDependencies(c(coefNode, target))

    calcNodesReduced <- model$getDependencies(target)

  },

  run = function() {

    currentIndicator <- model[[target]]

    if(currentIndicator == 0) {   ## propose addition of coefNode

      currentLogProb <- model$getLogProb(calcNodesReduced)

      proposalCoef <- numeric(length(coefNode))

      logProbForwardProposal <- 0

      for(l in 1:length(coefNode)){

        proposalCoef[l] <- rnorm(1, proposalMean, proposalScale)

        logProbForwardProposal <- logProbForwardProposal + dnorm(proposalCoef[l], proposalMean, proposalScale, log = TRUE)

      }

      model[[coefNodes]] <<- proposalCoef

      model[[target]] <<- 1

      proposalLogProb <- model$calculate(calcNodes)

      logAcceptanceProb <- proposalLogProb - currentLogProb - logProbForwardProposal

    } else {                      ## propose removal of coefNode

      currentLogProb <- model$getLogProb(calcNodes)

      currentCoef <-  model[[coefNodes]]

      logProbReverseProposal<- 0

      for(l in 1:length(coefNode)){

        logProbReverseProposal <- logProbReverseProposal + dnorm(currentCoef[l], proposalMean, proposalScale, log = TRUE)

      }

      model[[coefNodes]] <<- 0

      model[[target]] <<- 0

      model$calculate(calcNodes)

      logAcceptanceProb <- model$getLogProb(calcNodesReduced) - currentLogProb + logProbReverseProposal

    }

    accept <- decide(logAcceptanceProb)

    if(accept) { copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)

    } else     { copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }

  },

  methods = list(

    reset = function() { }

  )

)

code <- nimbleCode({

  beta0 ~ dnorm(0, sd = 100)

  beta[1] ~ dnorm(0, sd = 100)

  beta[2] ~ dnorm(0, sd = 100)

  beta[3] ~ dnorm(0, sd = 100)

  sigma ~ dunif(0, 100) 

  for(i in 1:2){

    z[i] ~ dbern(psi)

  }

  #z1 ~ dbern(psi)  ## indicator variable associated with beta1

  #z2 ~ dbern(psi)  ## indicator variable associated with beta2

  psi ~ dbeta(1, 1) ## hyperprior on inclusion probability

  for(i in 1:N) {

    Ypred[i] <- beta0 + z[1]*beta[1] * x1[i] + z[1] *beta[2] * x2[i] + beta[3] * x3[i] * z[2]

    Y[i] ~ dnorm(Ypred[i], sd = sigma)

  }

})

## simulate some data

set.seed(1)

N <- 100

x1 <- rbinom(N, 1, 0.5)

x2 <- rbinom(N, 1, 0.5) #

x3 = std2(runif(N, -10, 10))

Y <- rnorm(N, 1.5 + 2 * x1 + 0.4 * x2 + 0.8 * x3, sd = 1)

## build the model and configure default MCMC

RJexampleModel <- nimbleModel(code, constants = list(N = N),

                              data = list(Y = Y, x1 = x1, x2 = x2, x3=x3), 

                              inits = list(beta0 = 0, beta=c(0,0,0),

                                           sigma = sd(Y), z = rbinom(2,1,0.5), psi = 0.5))

cRJexampleModel <- compileNimble(RJexampleModel)

RJexampleConf <- configureMCMC(RJexampleModel)

configureRJ(conf = RJexampleConf,     

            targetNodes = c("beta[3]"),

            indicatorNodes = c('z[2]'),

            control = list(mean = c(0), scale = 2))

nodeControl  = list(mean = 0, scale = 1)

nodeControl$targetNode <- c("beta[1:2]")

RJexampleConf$removeSamplers("z[1]")

RJexampleConf$addSampler(type = my_sampler_RJ_indicator,

                target = "z[1]",

                control = nodeControl)

currentConf <- RJexampleConf$getSamplers("beta[1]")

RJexampleConf$removeSampler("beta[1]")

RJexampleConf$addSampler(type = sampler_RJ_toggled,

                         target = "beta[1]",

                         control = list(samplerType = currentConf[[1]]))

currentConf <- RJexampleConf$getSamplers("beta[2]")

RJexampleConf$removeSampler("beta[2]")

RJexampleConf$addSampler(type = sampler_RJ_toggled,

                         target = "beta[2]",

                         control = list(samplerType = currentConf[[1]]))

mpleConf$printSamplers()

RJexampleConf$addMonitors(c('z'))

mcmc = buildMCMC(RJexampleConf)

cmcmc5 = compileNimble(mcmc, project = RJexampleModel, resetFunctions = TRUE, showCompilerOutput = TRUE)

#tic()

results_samplerRJind <- runMCMC(cmcmc4, niter = 20000, nburnin=5000,thin=1,nchains=1, setSeed = 500)

Thanks!

Sincerely,

Daniel

[Daniel Eacker]

Here is the function I included to standardize by 2 SDs and center the continuous variable (Gelman 2007):

# simple function to standardize variables
    std2=function(x){(x - mean(x,na.rm=TRUE))/(2*sd(x,na.rm=TRUE))}

Thanks,

Daniel

[Brian Hatfield]

Hi Daniel,

I am attempting to use the my_sampler_RJ_indicator in your post above for a similar situation (using a single indicator for quadratic terms of a glm).  I know it has been a while, but I am wondering if you have found a solution to this problem?  When I run your my_sampler_RJ_indicator  code above I get the same error that you listed on compilation, (latest version of nimble, v 0.12.2):

in my case:

Compiling
  [Note] This may take a minute.

  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Error : $ operator is invalid for atomic vectors
Error: There is some problem at the Eigen processing step for this code:

 {
    currentIndicator <- model[[target]]
    if (currentIndicator == 0) {
        currentLogProb <- model$getLogProb(nodes = calcNodesReduced)
        proposalCoef <- nimNumeric(length = length(coefNode), value = 0, init = TRUE, fillZeros = TRUE, recycle = TRUE)
        logProbForwardProposal <- 0
        for (l in 1:length(coefNode)) {
            proposalCoef[l] <- rnorm(n = 1, mean = proposalMean, sd = proposalScale)
            logProbForwardProposal <- logProbForwardProposal + dnorm(x = proposalCoef[l], mean = proposalMean, sd = proposalScale, log = TRUE)

        }
        model[[coefNodes]] <<- proposalCoef
        model[[target]] <<- 1

        proposalLogProb <- model$calculate(nodes = calcNodes)

        logAcceptanceProb <- proposalLogProb - currentLogProb - logProbForwardProposal
    }
    else {

        currentLogProb <- model$getLogProb(nodes = calcNodes)

        currentCoef <- model[[coefNodes]]

Thank you,

Brian

[Perry de Valpine]

Hi Brian,

Thanks for the post.  Yikes, it looks like we (nimble development team) dropped the ball on this last fall.  We will try to look at it soon.

-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/ffc2fce8-b293-42ec-8c3e-970782543d84n%40googlegroups.com.

[Daniel Eacker]

Hi Brian,

I just reverted back an earlier version of nimble for now to get this to run.

Perry, If you get time to work on this problem that would be awesome! I'm sure there are a lot of people that will encounter the same problem when trying to use RJMCMC with factors involved.

Thanks,

Dan

[Brian Hatfield]

Hi Perry, Dan,

Thanks for getting back to me. A little more info that may be helpful?

Using Dan's code from the first post I tested just compiling the  my_sampler_RJ_indicator rather than the entire model:

my_sampler_RJ_indicator_test<-my_sampler_RJ_indicator(model = RJexampleModel,
                                                      mvSaved = RJexampleModel,

                                                      target = 'z[1]',
                                                      control = nodeControl
                                                      )

##this runs successfully :

my_sampler_RJ_indicator_test$run()

##compiling throws the same error:
my_sampler_RJ_indicator_test_c <- compileNimble(my_sampler_RJ_indicator_test)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Error : $ operator is invalid for atomic vectors
Error: There is some problem at the Eigen processing step for this code:
 {
    currentIndicator <- model[[target]]
    if (currentIndicator == 0) {
        currentLogProb <- model$getLogProb(nodes = calcNodesReduced)
        proposalCoef <- nimNumeric(length = length(coefNode), value = 0, init = TRUE, fillZeros = TRUE, recycle = TRUE)
        logProbForwardProposal <- 0
        for (l in 1:length(coefNode)) {
            proposalCoef[l] <- rnorm(n = 1, mean = proposalMean, sd = proposalScale)
            logProbForwardProposal <- logProbForwardProposal + dnorm(x = proposalCoef[l], mean = proposalMean, sd = proposalScale, log = TRUE)
        }
        model[[coefNodes]] <<- proposalCoef
        model[[target]] <<- 1
        proposalLogProb <- model$calculate(nodes = calcNodes)
        logAcceptanceProb <- proposalLogProb - currentLogProb - logProbForwardProposal
    }
    else {
        currentLogProb <- model$getLogProb(nodes = calcNodes)
        currentCoef <- model[[coefNodes]]

the one thing I tried was to move all the instances of coefNodes (in the my_sampler_RJ_indicator code) into the l loops, so I could remove the coefNodes and only use the coefNode..

my_sampler_RJ_indicator_2 <- nimbleFunction(

  name = 'my_sampler_RJ_indicator',
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    ## note: target is the indicator variable,
    ## control$targetNode is the variable conditionally in the model
    ## control list extraction
    coefNode      <- model$expandNodeNames(control$targetNode, returnScalarComponents = TRUE)

    #coefNodes     <- control$targetNode

    proposalScale <- control$scale
    proposalMean  <- control$mean
    ## node list generation
    calcNodes <- model$getDependencies(c(coefNode, target))
    calcNodesReduced <- model$getDependencies(target)
  },
  run = function() {
    currentIndicator <- model[[target]]
    if(currentIndicator == 0) {   ## propose addition of coefNode
      currentLogProb <- model$getLogProb(calcNodesReduced)
      proposalCoef <- numeric(length(coefNode))
      logProbForwardProposal <- 0

      for(l in 1:length(coefNode)){

        proposalCoef[l] <- rnorm(1, proposalMean, proposalScale)

        logProbForwardProposal <- logProbForwardProposal + dnorm(proposalCoef[l], proposalMean, proposalScale, log = TRUE)

        model[[coefNode[l]]] <<- proposalCoef[l]

      }
      #model[[coefNodes]] <<- proposalCoef

      model[[target]] <<- 1
      proposalLogProb <- model$calculate(calcNodes)
      logAcceptanceProb <- proposalLogProb - currentLogProb - logProbForwardProposal
    } else {                      ## propose removal of coefNode
      currentLogProb <- model$getLogProb(calcNodes)

      #currentCoef <-  model[[coefNodes]]
      logProbReverseProposal<- 0

      for(l in 1:length(coefNode)){

        currentCoef[l] <-  model[[coefNode[l]]]

        logProbReverseProposal <- logProbReverseProposal + dnorm(currentCoef[l], proposalMean, proposalScale, log = TRUE)

        model[[coefNode[l]]] <<- 0

        }
      #model[[coefNodes]] <<- 0

      model[[target]] <<- 0
      model$calculate(calcNodes)
      logAcceptanceProb <- model$getLogProb(calcNodesReduced) - currentLogProb + logProbReverseProposal
    }
    accept <- decide(logAcceptanceProb)
    if(accept) { copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
    } else     { copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }
  },
  methods = list(
    reset = function() { }
  )
)

my_sampler_RJ_indicator_2_test<-my_sampler_RJ_indicator_2(model = RJexampleModel,
                                                      mvSaved = RJexampleModel,

                                                      target = 'z[1]',
                                                      control = nodeControl
)

## again this runs without an error:

my_sampler_RJ_indicator_2_test$run()

## compiling throws a new error:

my_sampler_RJ_indicator_2_test_c <- compileNimble(my_sampler_RJ_indicator_2_test)

Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.

Error in exists(as.character(varOrNodeExpr), x, inherits = FALSE) :
  first argument has length > 1

I'm at a dead end with this error at the moment.

Brian

[Perry de Valpine]

Hi Brian and Dan,

Brian, thanks for the additional detective work.

Dan, thanks for your patience since we dropped this ball feel free to last fall.

All readers, we (nimble development team) really try to respond to all messages and not drop anyone's inquiry.  We do also sometimes get really busy, so please feel free to ping the list again or email one of us directly if we've dropped something. 

Below is a version of the code that works.  In the situation that coefNodes could be one or multiple nodes, it is better to use values(model, coefNodes) than model[[coefNodes]].  There is limited support for model[[coefNodes]] if coefNodes is something like (here) "beta[1:2]", which is really a node-name-like vector of two nodes.  However a difficulty would arise if sometimes (in some cases of the added sampler) coefNodes is a single scalar node and other times multiple nodes.  That's why the values(model, coefNodes) will be more reliable here.   In addition the coefNodes would be expanded internally, so I've used the expanded version you already have on hand, i.e. values(model, coefNode).  What evidently broke in a nimble update last September is assignment like "model[[nodes]] <- scalar" if nodes is not a single scalar node.  (The description of ways to access model nodes in nimbleFunction run code is in User Manual chapter 15.). We will need to pin that down and fix it, but, in the meantime, the code below should let you move ahead and not be stuck in an older version of nimble.

(Also, it is a bit better to get the length of coefNode is setup code and use that in run code.  If you use "length(coefNode)" in run code, then coefNode will become a member data object of the resulting C++ class and it will be repeatedly queried for its length.  In your code it is trivially short, so it doesn't really matter either way.)

Perry

my_sampler_RJ_indicator <- nimbleFunction(
  name = 'my_sampler_RJ_indicator',
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    ## note: target is the indicator variable,
    ## control$targetNode is the variable conditionally in the model
    ## control list extraction
    coefNode      <- model$expandNodeNames(control$targetNode, returnScalarComponents = TRUE)

   # coefNodes     <- control$targetNode

    proposalScale <- control$scale
    proposalMean  <- control$mean

    len_coefNode <- length(coefNode) # It is better to do this in setup code and use below

    ## node list generation
    calcNodes <- model$getDependencies(c(coefNode, target))
    calcNodesReduced <- model$getDependencies(target)
  },
  run = function() {
    currentIndicator <- model[[target]]
    if(currentIndicator == 0) {   ## propose addition of coefNode
      currentLogProb <- model$getLogProb(calcNodesReduced)

      proposalCoef <- numeric(len_coefNode)
      logProbForwardProposal <- 0
      for(l in 1:len_coefNode) {

        proposalCoef[l] <- rnorm(1, proposalMean, proposalScale)
        logProbForwardProposal <- logProbForwardProposal + dnorm(proposalCoef[l], proposalMean, proposalScale, log = TRUE)
      }

      values(model, coefNode) <<- proposalCoef

      model[[target]] <<- 1
      proposalLogProb <- model$calculate(calcNodes)
      logAcceptanceProb <- proposalLogProb - currentLogProb - logProbForwardProposal
    } else {                      ## propose removal of coefNode
      currentLogProb <- model$getLogProb(calcNodes)

      currentCoef <-  values(model, coefNode)
      logProbReverseProposal<- 0
      for(l in 1:len_coefNode) {

        logProbReverseProposal <- logProbReverseProposal + dnorm(currentCoef[l], proposalMean, proposalScale, log = TRUE)
      }

      values(model, coefNode) <<- rep(0, len_coefNode)

      model[[target]] <<- 0
      model$calculate(calcNodes)
      logAcceptanceProb <- model$getLogProb(calcNodesReduced) - currentLogProb + logProbReverseProposal
    }
    accept <- decide(logAcceptanceProb)
    if(accept) { copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
    } else     { copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) }
  },
  methods = list(
    reset = function() { }
  )
)

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/a609a51f-ca7f-45c1-b470-10aa50b315afn%40googlegroups.com.

[Daniel Eacker]

Thank you for diagnosing this Perry! I was able to move on with the analysis at the time so I wasn't worried about it, and it seems that the requests for help with nimble keep growing. That's just a reflection of how useful the nimble package is for researchers.

I'll check out your detailed solution and test out the code.

Sincerely,

Dan Eacker

[Brian Hatfield]

I just ran this on some data and it worked!

Thanks Perry for the solution, pretty sure I would not have figured this out on my own.  Ill try and note values() for future use.  Thanks Dan for the initial coding etc. 

Brian