Speed-up a model with Dirichlet processes involved

[Evgeny Lobzaev]

Hi All!

First of all, thank you for such an amazing tool as NIMBLE!

I am working with gene expression data and trying to fit a model that involves 2 Dirichlet processes. 

Essentially I am assuming that gene expression is distributed according to Negative Binomial distribution, particularly, this parametric from is used (often used in genetics): https://mc-stan.org/docs/2_20/functions-reference/nbalt.html. I assume that both parameters (on log-scale) are dependent on some covariates (X1,X2...). So it is essentially like a GLM but with covariates used for both parameters. I am trying to cluster the coefficient of this linear combinations separately for each distribution parameter and this is why I am using 2 separate Dirichlet processes.

The model code is the following:

modelCode <- nimbleCode({

  #### DP

  # for mean

  for(i in 1:(H-1)){

    v[i] ~ dbeta(1, alpha)

  }

  w[1:H] <- stick_breaking(v[1:(H-1)])

  

  #for phi

  for(i in 1:(HH-1)){

    vv[i] ~ dbeta(1, alpha2)

  }

  ww[1:HH] <- stick_breaking(vv[1:(HH-1)])

  

  ### sample from G0

  # mean

  for(h in 1:H){

    a0[h] ~ dnorm(0.0,1.5)

    a1[h] ~ ddexp(0.0,1.5)

    a2[h] ~ ddexp(0.0,1.5)

    a3[h] ~ ddexp(0.0,1.5)

    a4[h] ~ ddexp(0.0,1.5)

  }

  # phi

  for(hh in 1:HH){

    b0[hh] ~ dnorm(0.0,1.5)

    b1[hh] ~ ddexp(0.0,1.5)

    b2[hh] ~ ddexp(0.0,1.5)

    b3[hh] ~ ddexp(0.0,1.5)

    b4[hh] ~ ddexp(0.0,1.5)

  }

  

  

  # likelihood

  for(g in 1:G){

    # sample category: zeta[g] in [1......H], and zeta2[g] in [1...HH]

    zeta[g] ~ dcat(w[1:H])  

    zeta2[g] ~ dcat(ww[1:HH])

    

    for(t in 1:T){

      # sample gene expression

      Y[g,t] ~ dnegbin(p[g,t],phi[g,t])

      

      # construct value p[g,t] in terms of mu[g,t] and phi[g,t]

      p[g,t]   <- phi[g,t]/(scale_factors[t]*exp(logmu[g,t])+phi[g,t])

      phi[g,t] <- exp(logphi[g,t])

      

      # we need to define logmu[g,t]

      logmu[g,t] <- a0[zeta[g]] + a1[zeta[g]] * X1[t]+ a2[zeta[g]] * X2[t] + a3[zeta[g]]* X3[t]+ a4[zeta[g]]* X4[t]

      

      # define phi[g,t]

      logphi[g,t] <- b0[zeta2[g]] + b1[zeta2[g]] * X1[t] + b2[zeta2[g]] * X2[t] + b3[zeta2[g]] * X3[t] + b4[zeta2[g]] * X4[t]

    }

  }

})

The data passing and model compilation happens here:

constants  <- list(G      = G,

                   H      = n_clusters_mean,

                   HH     = n_clusters_phi,

                   T      = T,

                   alpha  = concentration_mean,

                   alpha2 = concentration_phi

                   )

data       <- list(Y             = Y, 

                   X1            = X1, 

                   X2            = X2,

                   X3            = X3,

                   X4            = X4,

                   scale_factors = scale_factors

                   )

inits      <- list(a0     = rnorm(constants$H, 0, 1), 

                   a1     = rnorm(constants$H, 0, 1),

                   a2     = rnorm(constants$H, 0, 1),

                   a3     = rnorm(constants$H, 0, 1),

                   a4     = rnorm(constants$H, 0, 1), 

                   

                   b0     = rnorm(constants$HH, 0, 1), 

                   b1     = rnorm(constants$HH, 0, 1),

                   b2     = rnorm(constants$HH, 0, 1), 

                   b3     = rnorm(constants$HH, 0, 1), 

                   b4     = rnorm(constants$HH, 0, 1), 

                   

                   v      = rbeta(constants$H, 1, 1),

                   vv     = rbeta(constants$HH, 1, 1)

                   )

# this is R object, not C++ object

Rmodel     <- nimbleModel(modelCode,constants,data,inits)

# lets monitor just cluster assignments z, we can change samplers here

monitors   <- c("zeta","zeta2")

RmodelConf <- configureMCMC(Rmodel, monitors = monitors,thin=lag)

# create mcmc object

Rmcmc      <- buildMCMC(RmodelConf) 

print("R objects created")

# compile everything in C++

Cmodel     <- compileNimble(Rmodel,showCompilerOutput = TRUE)

Cmcmc      <- compileNimble(Rmcmc,project = Cmodel,showCompilerOutput = TRUE)

print("C++ objects created")

runMCMC_samples <- runMCMC(Cmcmc, nburnin = burnin, niter = n_steps)

The model worked fine with small dataset (200 genes and 300 cells; ran for 5.5 hrs) but there were memory issues (and, as the result, it crashed on HPC; probably not enough memory was allocated) with 1000 genes and 250 cells. 

In fact the true input should be around 6K-by-600, so I am wondering what can I potentially do (if anything) to speed up the process? I read another thread here, where you suggested using calculate = False for nimbleModel(). I can also use vectorized form to compute logmu and logphi. I haven't checked it yet but I am about to. 

Is there anything else I can do to speed up the model and make it less memory-consuming?

Thank you for your help.

Regards,

Evgenii

[Perry de Valpine]

Hi Evgeny,

Thanks for the message and I'm very interested to see this application.

Let me point you to my message of June 26 on this thread:

https://groups.google.com/forum/#!topic/nimble-users/8vaTHQoJz2M

I commented on why nimble encounters memory problems (which we hope to improve).  Then I gave a strategy for keeping large objects out of the model.  That strategy might work well for your problem.

Looking at your model, I wonder if the following would work: Vectorize logmu and logphi over t (since any change in a parameter will affect all t).  And use a vector of dummy data values to push calculation of phi[g, 1:T], p[g, 1:T] and log probabilities of y[g, 1:T] to R (or C++) using the trick I showed in the June 26 message.  Your model code would have something like:

dummy[g] ~ dCustom(logmu[g, 1:T],  logphi[g, 1:T], g ) # do the phi, p and y steps all together

All the actual values of y would live in dCustom.

This could also be done for logmu and/or logphi at a couple of levels.  You could have

logmu[g, 1:T] <- calc_logmu(a0[zeta[g]], a1[zeta[g]], etc.)

All of the values of X would live in calc_logmu.

Going further, you could store logmu itself outside the model, using another dummy variable to serve as a low-dimensional node to establish the graph relationships in the model.  This would look something like

logmu_dummy[g] <- calc_logmu(g, a0[zeta[g]], a1[zeta[g]], etc.)

where calc_logmu would have to return a value to logmu_dummy[g] that wouldn't really matter.  The actual values of logmu would be stored in R (or C++) in the same scope from where they'd be used by dCustom.

This is all something of a workaround to the scaling issues, but it is fairly straightforward to set up once you see how the pieces fit together and allows nimble to manage the model graph and MCMC system.   I'm happy to help sort it out if it's not clear how it all can fit together.

A couple of final thoughts on your model: Sampling for dcat can be inefficient because on every iteration every possible value is tried.  It would be easy to set up a sampler that for example tries one alternative or a random subset of alternatives.  Computational cost vs. mixing is hard to predict and would depend on the model.  The dcat case is a tradeoff I've been curious about. Also, you may be good with what you've written, but nimble does have Bayesian non-parametric distributions (see blog posts linked from our documentation page at r-nimble.org) and they will be improved in an upcoming release.

HTH and let me know how it goes.

-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/03185087-8187-4a8a-ac75-d99eb60383efo%40googlegroups.com.

[Evgeny Lobzaev]

Hi Perry!

Thank you for such a quick reply!

That 26th of June post was exactly the one I was referring to me in my initial post.

I think before I ask for additional help, I will struggle for a bit myself (and then ask for help if I do not manage to make it work!) as this will allow me to better understand how NIMBLE works. 

One thing that I'd like to ask right now is about BNP models in NIMBLE. You are saying that NIMBLE has BNP models. But this is essentially what I am using right now: I just took a stick-breaking DP representation from chapter 10.3.2 in the documentation. In that example categorical distribution is used. Are you suggesting to replace the stick-breaking representation with the CRP (to avoid using categorical distribution)? Or did you mean something else?

Evgenii

On Wednesday, July 15, 2020 at 3:41:48 PM UTC+1, pdevalpine wrote:

Hi Evgeny,

Thanks for the message and I'm very interested to see this application.

Let me point you to my message of June 26 on this thread:

https://groups.google.com/forum/#!topic/nimble-users/8vaTHQoJz2M

I commented on why nimble encounters memory problems (which we hope to improve).  Then I gave a strategy for keeping large objects out of the model.  That strategy might work well for your problem.

Looking at your model, I wonder if the following would work: Vectorize logmu and logphi over t (since any change in a parameter will affect all t).  And use a vector of dummy data values to push calculation of phi[g, 1:T], p[g, 1:T] and log probabilities of y[g, 1:T] to R (or C++) using the trick I showed in the June 26 message.  Your model code would have something like:

dummy[g] ~ dCustom(logmu[g, 1:T],  logphi[g, 1:T], g ) # do the phi, p and y steps all together

All the actual values of y would live in dCustom.

This could also be done for logmu and/or logphi at a couple of levels.  You could have

logmu[g, 1:T] <- calc_logmu(a0[zeta[g]], a1[zeta[g]], etc.)

All of the values of X would live in calc_logmu.

Going further, you could store logmu itself outside the model, using another dummy variable to serve as a low-dimensional node to establish the graph relationships in the model.  This would look something like

logmu_dummy[g] <- calc_logmu(g, a0[zeta[g]], a1[zeta[g]], etc.)

where calc_logmu would have to return a value to logmu_dummy[g] that wouldn't really matter.  The actual values of logmu would be stored in R (or C++) in the same scope from where they'd be used by dCustom.

This is all something of a workaround to the scaling issues, but it is fairly straightforward to set up once you see how the pieces fit together and allows nimble to manage the model graph and MCMC system.   I'm happy to help sort it out if it's not clear how it all can fit together.

A couple of final thoughts on your model: Sampling for dcat can be inefficient because on every iteration every possible value is tried.  It would be easy to set up a sampler that for example tries one alternative or a random subset of alternatives.  Computational cost vs. mixing is hard to predict and would depend on the model.  The dcat case is a tradeoff I've been curious about. Also, you may be good with what you've written, but nimble does have Bayesian non-parametric distributions (see blog posts linked from our documentation page at r-nimble.org) and they will be improved in an upcoming release.

HTH and let me know how it goes.

-Perry

To unsubscribe from this group and stop receiving emails from it, send an email to nimble...@googlegroups.com.

[Daniel Turek]

I'm not an expert on the BNP facilities offered in NIMBLE, but I think I can comment: the stick-breaking representation (as you're using), and the CRP distribution, are modelling the same Dirichlet process.  They differ in the MCMC sampling algorithms which are applied to each, automatically, by NIMBLE's MCMC.  Translating your code between these two (mathematically identical) representations is a relatively small and straight-forward exercise.  The relative performance of the different sampling algorithms, as applied to each representation, could differ, would depend on the model and the data, and is generally difficult to predict.

Changing between the stick-breaking representation and the CRP representation, the categorical distribution would persist in your model.  And Perry's comments are correct, that (in particular for a large number of categories, e.g. large values of H or HH in your code) the categorical sampler (which is applied to the categorical distribution) would become arbitrarily inefficient, in that it evaluates the posterior density for every category value (H or HH times), then samples directly from that posterior distribution.  And that you could hopefully reduce this inefficiency by writing a customized sampling strategy (a custom-written MCMC sampler), to update your categorical distributions, in a more efficient (less wasteful) manner.

Hope this helps,

Daniel

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/f2afb677-391e-4ff4-b798-5285d6c5ff69o%40googlegroups.com.

[Evgeny Lobzaev]

I see. 

Well, to be honest, I have no idea how to write such a customized sampler? Like where do I start if I have never done it before?

Evgenii

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/f2afb677-391e-4ff4-b798-5285d6c5ff69o%40googlegroups.com.

[Daniel Turek]

You may or may not need to go down the road of writing a customized sampler for your model.  It really depends on your model, whether this is worthwhile.  I would see how performance (speed and mixing) of the MCMC is sampling is, before starting with that.

But, if you do decide to look into it, then I'd suggest first reading section15.5 of the NIMBLE user's manual ("Example: writing user-defined samplers to extend NIMBLE’s MCMC engine")

https://r-nimble.org/html_manual/cha-progr-with-models.html#sec:user-samplers

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/55f6e41a-43c5-4ae7-846a-efdcc2883754o%40googlegroups.com.

[Evgeny Lobzaev]

Ok, need to take a look at that chapter! Thank you, Daniel.

Perry,

I am not sure I managed to incorporate all your ideas in the code (as I did not fully understand all of the points). However, I did the following and it seemed to improve the speed. I am still wondering if additional improvement can be achieved (particularly, I am still sending count data matrix Y to NIMBLE, I did not understand how to keep it in R as you were trying to suggest).

Here is the code:

makeRfuns <- function(keepInR){

  

  # get samples

  X1 <- keepInR$X1

  X2 <- keepInR$X2

  X3 <- keepInR$X3

  

  # I can use X1,X2,X3 in all functions here

  # generic function, can be used for both logmu and logphi

  calc_log <- function(c0,c1,c2,c3){

    return(c0 + c1 * X1 + c2 * X2 + c3 * X3)

  }

  return(list(calc_log = calc_log))

}

####### multivariate negative binomial distribution

# density

dNBVector <- nimbleFunction(

  run = function(x      = double(1),

                 logmu  = double(1),

                 logphi = double(1),

                 log    = integer(0, default = 0)) { ## default should be zero

    returnType(double(0))

    # compute p,r out of logmu and logphi

    r       = exp(logphi)

    p       = r/(exp(logmu)+r)

    logProb = sum(dnbinom(x,r,p,log=TRUE))

    if(log) return(logProb) else return(exp(logProb))

  }

)

# random number generation

rNBVector<-nimbleFunction(

  run = function(n      = integer(0),

                 logmu  = double(1),

                 logphi = double(1)) { ## default should be zero

    returnType(double(1))

    # compute p,r out of logmu and logphi

    r  = exp(logphi)

    p  = r/(exp(logmu)+r)

    return(rnbinom(n,r,p))

  }

)

##### obtain Y,X1,X2,X3 somehow

# create a function that computes logmu/logphi

keepInR <- list(X1=X1,X2=X2,X3=X3)

# obtain pointer to calc_log function

Rfuns             <- makeRfuns(keepInR)

calc_log_internal <- Rfuns$calc_log

# now need to register with NIMBLE

calc_log <- nimbleRcall(function(c0=double(0),c1=double(0),c2=double(0),c3=double(0)){},

                        returnType = double(1),

                        Rfun = "calc_log_internal"

                       )

#### create nimble model

modelCode <- nimbleCode({

  #### DP

  # for mean

  for(i in 1:(H-1)){

    v[i] ~ dbeta(1, alpha)

  }

  w[1:H] <- stick_breaking(v[1:(H-1)])

  

  #for phi

  for(i in 1:(HH-1)){

    vv[i] ~ dbeta(1, alpha2)

  }

  ww[1:HH] <- stick_breaking(vv[1:(HH-1)])

  

  ### sample from G0

  # mean

  for(h in 1:H){

    a0[h] ~ dnorm(0.0,1.5)

    a1[h] ~ dnorm(0.0,1.5)

    a2[h] ~ dnorm(0.0,1.5)

    a3[h] ~ dnorm(0.0,1.5)

  }

  # phi

  for(hh in 1:HH){

    b0[hh] ~ dnorm(0.0,1.5)

    b1[hh] ~ dnorm(0.0,1.5)

    b2[hh] ~ dnorm(0.0,1.5)

    b3[hh] ~ dnorm(0.0,1.5)

  }

  

  for(g in 1:G){

    # sample category: zeta[g] in [1......H], and zeta2[g] in [1...HH]

    zeta[g] ~ dcat(w[1:H])  

    zeta2[g] ~ dcat(ww[1:HH])

    

    # try to use our function

    logmu[g,1:T]  <- calc_log(a0[zeta[g]],a1[zeta[g]],a2[zeta[g]],a3[zeta[g]])

    logphi[g,1:T] <- calc_log(b0[zeta2[g]],b1[zeta2[g]],b2[zeta2[g]],b3[zeta2[g]])

    

    # likelihood

    Y[g,1:T] ~ dNBVector(logmu[g,1:T],logphi[g,1:T])

    

  }

})

Basically, I did not understand how to return dummy values (either deterministic or stochastic) and keep all the data (Y,X1,X2...) inside a particular scope (I think this is what you meant, but I can be wrong).

Evgenii

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/55f6e41a-43c5-4ae7-846a-efdcc2883754o%40googlegroups.com.

[Claudia Wehrhahn Cortes]

Hi Evgeny,

I am really interested to see this application too!

Right, your model is a nonparametric model with a stick-breaking representation for the random measures. Using the CRP representation might be more efficient (well this depends on what values you are using for H and HH and how many cluster are determine by the CRP distribution). I think is worth trying it (see code below). 

The functionality for using the CRP in a model such as yours is not yet in our released version (but should be in a month or so) but  you can install the nimble version from devel to use the CRP for your model doing this:

Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS="true")

remotes::install_github('nimble-dev/nimble','devel', 'packages/nimble')

If you use the CRP representation you can define H1 < G and HH1 < G cluster parameters (a0, a1, a2, a3, a4, and b0, b1, b2, b3, b4) to reduce compilation time. This will create a warning in the “buildMCMC” step and if ever the sampler needs more cluster parameters than the ones you have available, the sampler will warn you and you should increase the values of H1 and HH1. Also, the CRP has the advantage that only cluster parameters that are occupied are updated so this reduces the run time of the MCMC (again, how much it reduces depends on how many clusters are created).

Under this representation the random measures are integrated out so you don’t get posterior samples of them directly from the MCMC’s output but you can use the getSamplesDPmeasure function which will provide posterior samples of the random measures.

One more thing, if you use the CRP representation and consider gamma prior distributions for the concentration parameters NIMBLE will assign specialized samplers for them. Recall that this parameter affects the number of clusters that are created.

 

I hope this comments together with Perry’s and Daniel’s inputs are helpful.

Please, let me know if you have any questions or something is not clear enough.

----------------------------------------------------------

Here is how your code would look like:

modelCodeCRP <- nimbleCode({

  # if you want to place prior distributions for the concentration parameter:

  alpha ~ dgamma(c1, d1) 

  alpha2 ~ dgamma(c2, d2) 

  

  for(h in 1:H1){ # H1 can be smaller or equal to G

    a0[h] ~ dnorm(0.0,1.5)

    a1[h] ~ ddexp(0.0,1.5)

    a2[h] ~ ddexp(0.0,1.5)

    a3[h] ~ ddexp(0.0,1.5)

    a4[h] ~ ddexp(0.0,1.5)

  }

  for(hh in 1:HH1){ # HH1 can be smaller or equal to G

    b0[hh] ~ dnorm(0.0,1.5)

    b1[hh] ~ ddexp(0.0,1.5)

    b2[hh] ~ ddexp(0.0,1.5)

    b3[hh] ~ ddexp(0.0,1.5)

    b4[hh] ~ ddexp(0.0,1.5)

  }

  

  zeta[1:G] ~ dCRP(conc = alpha, size = G)  

  zeta2[1:G] ~ dCRP(conc = alpha2, size = G) 

  

  # likelihood

  for(g in 1:G){

    

    for(t in 1:T0){

      Y[g,t] ~ dnegbin(p[g,t],phi[g,t])

      

      p[g,t]   <- phi[g,t]/(scale_factors[t]*exp(logmu[g,t])+phi[g,t])

      phi[g,t] <- exp(logphi[g,t])

      

      logmu[g,t] <- a0[zeta[g]] + a1[zeta[g]] * X1[t]+ a2[zeta[g]] * X2[t] + a3[zeta[g]]* X3[t]+ a4[zeta[g]]* X4[t]

      

      logphi[g,t] <- b0[zeta2[g]] + b1[zeta2[g]] * X1[t] + b2[zeta2[g]] * X2[t] + b3[zeta2[g]] * X3[t] + b4[zeta2[g]] * X4[t]

    }

  }

})

Best,

Claudia

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/222a5526-84a6-4570-bf49-953e692a2c17o%40googlegroups.com.

[Evgeny Lobzaev]

Hi Claudia,

Thank you for your response. I may be doing something wrong but I am receiving the following error:

• sampler_CRP: At the moment, NIMBLE can only sample when there is one variable being clustered for each cluster membershipID. When I use non-vectorized code (i.e. as the one in your snippet).
  
• sampler_CRP: CRP_nonconjugate sampler not yet set up to handle clustering of more than three variables. When I use vectorized form (i.e. as the one which I posted in one of the messages above).

I have installed the dev version of NIMBLE (it shows as 0.10).

Evgenii

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/222a5526-84a6-4570-bf49-953e692a2c17o%40googlegroups.com.

[Chris Paciorek]

Hi Evgeny,

I think something went wrong with your installation from the devel branch, but I'm not sure what if it does show that it is version 0.10.0.

That error message is the old one from version 0.9.1, and the message is not in the current code in our devel branch. You might try uninstalling your current nimble and

installing again using the instructions Claudia provided.

Alternatively, if you are familiar with Git you might clone our repo from github.com/nimble-dev/nimble, then change to the directory "nimble/packages" and then do:

R CMD build nimble

R CMD INSTALL nimble_0.10.0.tar.gz 

to install nimble.

-chris

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/96622458-5a19-4425-9735-13e37e101adao%40googlegroups.com.

[Evgeny Lobzaev]

Hi Chris,

Ok, that may have helped with previous issue. But the new error shows up in configureMCMC() function: Error in unrolledIndices[, loopIndex[varIdx]] : 

  invalid subscript type 'list'. Any suggestion what it might be? I am using the exact model that Claudia presented except I don't put priors on the CRP concentration parameter (i.e. nothing is vectorized for now) but that is clearly irrelevant to the error.

Evgenii

On Thursday, July 16, 2020 at 4:16:28 PM UTC+1, Chris Paciorek wrote:

Hi Evgeny,

I think something went wrong with your installation from the devel branch, but I'm not sure what if it does show that it is version 0.10.0.

That error message is the old one from version 0.9.1, and the message is not in the current code in our devel branch. You might try uninstalling your current nimble and

installing again using the instructions Claudia provided.

Alternatively, if you are familiar with Git you might clone our repo from github.com/nimble-dev/nimble, then change to the directory "nimble/packages" and then do:

R CMD build nimble

R CMD INSTALL nimble_0.10.0.tar.gz 

to install nimble.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/96622458-5a19-4425-9735-13e37e101adao%40googlegroups.com.

[Chris Paciorek]

That's strange -- I just tried Claudia's code (also without priors on alpha and alpha2), and I don't get that error. Can you provide the actual code you are using, with the constants and data you use?

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/44ba2783-31cc-4ee7-b8ce-10f945acbcefo%40googlegroups.com.

[Evgeny Lobzaev]

Hi Chris!

Absolutely.

Evgenii

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/44ba2783-31cc-4ee7-b8ce-10f945acbcefo%40googlegroups.com.

[Perry de Valpine]

Hi Evgenii,

Well, this has been a very interesting case on multiple fronts.  Let me respond on the issues of vectorizing and keeping large objects out of the model and in R instead.  It looks like you made a lot of progress on both.  Great to see that you got those pieces working.

Let me explain more about what happens with a stochastic declaration and what I meant about keeping the large Y matrix out of the model.

A stochastic declaration like

Y[g,1:T] ~ dNBVector(logmu[g,1:T],logphi[g,1:T])    #(I am just using the most recent code in front of me.)

tells nimble that Y[g, 1:T] is a (vector, in this case) stochastic node.  This means two things: Y[g, 1:T] will be built into the graph for determining dependencies among nodes.  And the way to calculate the log probability is:

logProb_Y[g, 1] <- dNBVector(Y[g, 1:T], logmu[g,1:T],logphi[g,1:T])

So you could do the following:

proxy_Y[g] ~ dNBVector(logmu[g,1:T],logphi[g,1:T], g)

provide proxy_Y[g] as data to the model.  Its values won't matter.

Now the two things above will be: proxy_Y[g] will be built into the graph, doing the job of Y[g, 1:T.  For example, proxy_Y[g] will be a dependency of logmu[g, 1:T].  And the way to calculate the log probability will be:

logProb_proxy_Y[g] <- dNBVector(proxy_Y[g], logmu[g,1:T], logphi[g,1:T], g)

Then you can create dNBVector as follows:

Add it to the list of functions returned from makeRfuns.

Include Y in the keepInR list so dNBVector can access its values by scoping just as calc_log accesses X1, X2, X3.

Use the same steps as you use for calc_log, i.e. dNBVector_internal and then a nimbleRcall called dNBVector.

Notice that if dNBVector takes g as an argument (as shown above), then it can pick out the correct row of Y to use.

dNBVector must take proxy_Y[g] as its first argument, but it can just ignore it.

Have dNBVector return the log probability of all of Y[g, 1:T].  This will be stored as the log probability of proxy_Y[g].

I hope it makes enough sense to follow!  If I had more time at the moment I would draft it out more fully for you, but let me know if it doesn't work smoothly.

-Perry

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/222a5526-84a6-4570-bf49-953e692a2c17o%40googlegroups.com.

[Evgeny Lobzaev]

Hi Perry!

Thank you for detailed explanation. I think I understood what needs to be done. I created a code and it works.

I have few concerns/points to mention.

1.  When I declared a proxy_Y node, I had to make it vector-valued(i.e. proxy_Y[g,1:T] instead of suggested proxy_Y[g]). Otherwise, the c++ compilation fails. 
2. I tested 2 instances of the code: one where I keep matrix Y in R and one where I supply matrix Y to nimble. The second case surprisingly runs a bit faster. Is it because of C++ <-> R  data passing overhead or what can be the reason?
3. In both versions I had to provide not just a density function for multivariate Negative Binomial distribution (dNBVector) but also a random number generating function (rNBVector). In 2 different situations (provide/not provide Y to nimble) the dNBVector function is coded slightly differently. I wonder if that can contribute to performance issues as well?

In fact, I will share both versions (that is whether or not Y is supplied to nimble) and, if you have some time/desire, you can take a look.

I have to say, this is very kind of you and other team members of NIMBLE to help me. You are really helping me with my project and I appreciate that a lot!

Evgenii

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/222a5526-84a6-4570-bf49-953e692a2c17o%40googlegroups.com.

[Chris Paciorek]

Hi Evgeny,

On the question of the error when you use the CRP approach, there was
a difference between NIMBLE's behavior on R 3.6 and R 4.0 with regard
to
setting up the MCMC sampler for the CRP node, and I suspect you are using R 4.0.
I've fixed that discrepancy, so if you re-install from the devel
branch, I think it should be fine for you now.

-chris

> 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/1b309a05-501f-4836-a5e8-a738bef7b834o%40googlegroups.com.

[Evgeny Lobzaev]

Hi Chris,

Indeed it was R 4.0.2 (just didn't think it would matter).

It all works now!

Evgenii

On Friday, July 17, 2020 at 3:45:48 PM UTC+1, Chris Paciorek wrote:

Hi Evgeny,

On the question of the error when you use the CRP approach, there was
a difference between NIMBLE's behavior on R 3.6 and R 4.0 with regard
to
setting up the MCMC sampler for the CRP node, and I suspect you are using R 4.0.
I've fixed that discrepancy, so if you re-install from the devel
branch, I think it should be fine for you now.

-chris

> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1b309a05-501f-4836-a5e8-a738bef7b834o%40googlegroups.com.

[Perry de Valpine]

Thanks for reporting on some experimental results.

On Fri, Jul 17, 2020 at 3:29 AM Evgeny Lobzaev <evgeny....@gmail.com> wrote:

Hi Perry!

Thank you for detailed explanation. I think I understood what needs to be done. I created a code and it works.

I have few concerns/points to mention.

1.  When I declared a proxy_Y node, I had to make it vector-valued(i.e. proxy_Y[g,1:T] instead of suggested proxy_Y[g]). Otherwise, the c++ compilation fails. 

This doesn't make sense to me.  In the nimbleRcall prototype are you changing x=double(1) to x=double(0) to tell it the data value will be a scalar?

 

1. I tested 2 instances of the code: one where I keep matrix Y in R and one where I supply matrix Y to nimble. The second case surprisingly runs a bit faster. Is it because of C++ <-> R  data passing overhead or what can be the reason?

It could be the copy overhead.  It could also be the cost of the R evaluator.  In theory things should generally run faster in compiled nimble than in R.  The point of going to R was to avoid scaling problems in nimble.  So, yes and yes, a mix of tradeoffs involved. 

1. In both versions I had to provide not just a density function for multivariate Negative Binomial distribution (dNBVector) but also a random number generating function (rNBVector). In 2 different situations (provide/not provide Y to nimble) the dNBVector function is coded slightly differently. I wonder if that can contribute to performance issues as well?

Did you mean to say the rNBVector was coded slightly differently?  The "r" function typically isn't used much, unless there are no inits or there are posterior predictive nodes following the distribution.  If the dNBVector is coded differently, yes that could make a difference -- but I think that is self-evident.

 

In fact, I will share both versions (that is whether or not Y is supplied to nimble) and, if you have some time/desire, you can take a look.

I have to say, this is very kind of you and other team members of NIMBLE to help me. You are really helping me with my project and I appreciate that a lot!

Great! Glad to help.  Limited time to dig further into the performance issues right now.  Glad you're up and running.

-Perry

 

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/1b309a05-501f-4836-a5e8-a738bef7b834o%40googlegroups.com.

[Claudia Wehrhahn Cortes]

Hi Evgeny,

I'm glad everyhting is working now, please let us know if you have further questions.

Best,

Claudia 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CABeTmqherSWXeKzYxMD4nx7B-Xk%2B-CFi1n62AfLf9sYUmssZXw%40mail.gmail.com.

[Evgeny Lobzaev]

Thank you so much guys!

You really helped a lot! Hope I don't need to bother you again!

Evgenii

On Friday, July 17, 2020 at 8:25:34 PM UTC+3, Claudia Wehrhahn Cortes wrote:

Hi Evgeny,

I'm glad everyhting is working now, please let us know if you have further questions.

Best,

Claudia 

On Fri, Jul 17, 2020 at 9:45 AM Perry de Valpine <pdeva...@berkeley.edu> wrote:

Thanks for reporting on some experimental results.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1b309a05-501f-4836-a5e8-a738bef7b834o%40googlegroups.com.

--
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...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CABeTmqherSWXeKzYxMD4nx7B-Xk%2B-CFi1n62AfLf9sYUmssZXw%40mail.gmail.com.

[Claudia Wehrhahn Cortes]

Hi Evgeny,

you are very welcome. And please, keep bothering if required! I am really interested to know what NIMBLE users will do with the BNP functionality!

Best,

Claudia

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/a4069848-ad2b-482c-a8b5-7ccc4915237co%40googlegroups.com.