How to create local variable in nimble code

[Lu Zhang]

Hello,

     I am working on the implementation of a hierarchical spatial model in nimble.  To get the posterior simulations of the spatial residual, I need to update covariance matrix in a multivariate normal distribution. The function for generating covariance matrix is

CovMat <- nimbleFunction(

  run = function(phi, dist = double(2), N, sigmasq) {

    

    # function for getting covariance matrix: sigmasq*exp(-phi*dist)

    

    CovM <- matrix(0, ncol = N, nrow = N)

    

    for (j in 1:N){

      for (k in 1:N){

        CovM[j, k] <- sigmasq * exp(-phi * dist[j, k]) 

      }

    }

    

    returnType(double(2))

    return(CovM[1:N, 1:N])

  }

)

When I called this function directly in my nimble code:

FullGPCode <- nimbleCode({

  

  w[1:N] ~ dmnorm( meanv[1:N], cov = CovMat(phi, D[1:N, 1:N], N, 1 / invsigmasq))

  for (i in 1:N){

    y[i] ~ dnorm(inprod(x[i, 1:2], beta[1:2]) + w[i], sd = sqrt(1 / invtausq))

  }

  

  # Hyperparameters:

  phi ~ dgamma(ap, rate = bp)                    # 1/range

  invsigmasq ~ dgamma(as, rate = bs)             # partial sill

  invtausq ~ dgamma(at, rate = bt)               # nugguts

  beta[1:2] ~ dmnorm(mu_b[1:2], V_b[1:2 , 1:2])  # regression coefficients

}

)

I got this error:

Error in checkMultivarExpr() : 

  Error with parameter 'cov' of distribution 'dmnorm': multivariate parameters cannot be expressions; please define the expression as a separate deterministic variable and use that variable as the parameter.

Then I tried to use a local variable "Covm" to store the updated covariance matrix and then call in my multivariate normal distribution:

FullGPCode <- nimbleCode({

  

  Covm[1:N, 1:N] <-  CovMat(phi, D[1:N, 1:N], N, 1 / invsigmasq)

  w[1:N] ~ dmnorm( meanv[1:N], cov = Covm[1:N, 1:N])

  for (i in 1:N){

    y[i] ~ dnorm(inprod(x[i, 1:2], beta[1:2]) + w[i], sd = sqrt(1 / invtausq))

  }

  

  # Hyperparameters:

  phi ~ dgamma(ap, rate = bp)                    # 1/range

  invsigmasq ~ dgamma(as, rate = bs)             # partial sill

  invtausq ~ dgamma(at, rate = bt)               # nugguts

  beta[1:2] ~ dmnorm(mu_b[1:2], V_b[1:2 , 1:2])  # regression coefficients

}

)

I got this error when I tried to compile the code

compiling... this may take a minute. On some systems there may be some compiler warnings that can be safely ignored.

Error in x[[1]] : object of type 'symbol' is not subsettable

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

 {

    model$Covm[1:200, 1:200] <<- CovMat(phi = model$phi[1], dist = model$D[1:200, 1:200], N = 200, sigmasq = 1/model$invsigmasq[1])

}

So my questions are: 

   1. Is there a way to store a matrix output from a user-defined function in a local variable in nimble code? 

   2. How to define a separate deterministic variable and use that variable as the parameter (from the first error)? 

Thank you!

Lu

[Perry de Valpine]

Hi Lu,

Thanks for the query. I see two small issues with your code, one of which you fixed but then the second one tripped you up.

Your fix to do this:

Covm[1:N, 1:N] <-  CovMat(phi, D[1:N, 1:N], N, 1 / invsigmasq)

  w[1:N] ~ dmnorm( meanv[1:N], cov = Covm[1:N, 1:N])

is the correct idea.

But I think the next problem arose because you need to declare a type for all arguments to the run function in CovMat, even if they are scalars.  I think you’ll want phi = double(), N = integer(), sigmasq = double().  Note that double() and double(0) are equivalent.  Your declaration for dist looks good.

For those interested in why the call the CovMat can’t be written directly as the cov argument: it is because nimble’s processing of BUGS code does not do as much size processing as one would like, but we must be sure that what is provided for cov is valid and one way to do that is to limit it to only a variable name with index ranges.  We’ll be able to do better when we find some time in the future, but even so there will be tricky cases because for BUGS code we must determine the sizes of all objects in order to build the model correctly.

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 post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/ef98c0cc-4b9d-47ae-a528-7fa109dd70fc%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Lu Zhang]

Hi Perry,

    Cool! Thank you very much for your timely reply. I will try in my code.

Best,

Lu

[Lu Zhang]

Update: 

 Thank you for Perry's help. Now I can compile my code. However, the simulations I got behave very weird. This is the plot of the samples for w[1]. It did not update. And the estimators of phi and sigmasq are almost 0. 

The simulations of the intercept and slope are good. So I think my code may fail to pass the correct covariance matrix when generating the MCMC chain. So would anyone help me finding the bug in my code? 

The following is my code: 

library("nimble")

CovMat <- nimbleFunction(

  run = function(phi = double(), dist = double(2), N = integer(), 

                 sigmasq = double()) {

    

    # function for getting covariance matrix: sigmasq*exp(-phi*dist)

    

    CovM <- matrix(NA, ncol = N, nrow = N)

    

    for (j in 1:N){

      for (k in 1:N){

        CovM[j, k] <- sigmasq * exp(-phi * dist[j, k]) 

      }

    }

    

    returnType(double(2))

    return(CovM[1:N, 1:N])

  }

)

FullGPCode <- nimbleCode({

  
    # This is the part I think fails to work properly:

    Covm[1:N, 1:N] <- CovMat(phi, D[1:N, 1:N], N, (1 / invsigmasq))

    w[1:N] ~ dmnorm(meanv[1:N], cov = Covm[1:N, 1:N])

    for (i in 1:N){

      y[i] ~ dnorm(inprod(x[i, 1:2], beta[1:2]) + w[i], sd = sqrt(1 / invtausq))

    }

  

    # Hyperparameters: parameters in blue do not get correct simulation

    phi ~ dgamma(ap, rate = bp)                    # 1/range

    invsigmasq ~ dgamma(as, rate = bs)             # partial sill

    invtausq ~ dgamma(at, rate = bt)               # nugguts

    beta[1:2] ~ dmnorm(mu_b[1:2], V_b[1:2, 1:2])   # regression coefficients

  }

)

Constants <- list(N = n, P = 2, ap = 12, bp = 1, as = 2 , bs = 1, 

                  at = 2, bt = 0.1, mu_b = c( 1 , 5 ), 

                  V_b = matrix (c(100, 0, 0, 10000), nrow = 2 , ncol = 2),

                  D = D[1:n, 1:n], meanv = rep(0, n))

## list specifying model data

Data <- list(y = y, x = X)

## list specifying initial values

Inits <- list(invsigmasq = 2, invtausq = 5, 

              beta = c(5.0, 1.0), phi = 5, w = rep(0.01, n))

## build the R model object

FullGP_build <- nimbleModel(code = FullGPCode, name = 'FullGP',

                          constants = Constants, data = Data, inits = Inits)

CFullGP <- compileNimble(FullGP_build)

###########################################

##### MCMC configuration and building #####

###########################################

## generate the default MCMC configuration

FullGP_conf <- configureMCMC(FullGP_build, print = T)

FullGP_conf$addMonitors(c('invsigmasq', 'invtausq', 'beta','phi', 'w'))

## check the samplers assigned by default MCMC configuration

FullGP_conf$printSamplers()

## double-check our monitors, and thinning interval

FullGP_conf$getMonitors()

## build the executable R MCMC function

build_mcmc_FullGP <- buildMCMC(FullGP_conf)

###################################

##### compile to C++, and run #####

###################################

## compile the the model

CFullGPMCMC <- compileNimble(build_mcmc_FullGP, project = FullGP_build)

## run the default MCMC function

niter <- 5000

set.seed(0)

t<-proc.time()

CFullGPMCMC$run(niter)

proc.time() - t

Fsamples <- as.matrix(CFullGPMCMC$mvSamples)

Thank you!

Lu

 

[Perry de Valpine]

Hi Lu,

Can you send an example that is fully reproducible, with a value of n and generation of data if possible?

If I follow, w[1:N] is a vector random effect. By default it will be sampled by a block sampler (I don’t think conjugacy with the y[i] would be determined in this case), and adaptation of a block sampler can sometimes not work well.  A couple of suggestions would be:

You could add univariate samplers for scalar elements of w.  You could do that like this:

for(node in FullGP_build$expandNodeNames(‘w’, returnScalarComponents = TRUE)) 

  FullGP_conf$addSampler(target = node, type = ‘RW’)

Alternatively, there is an argument to configureMCMC called multivariateNodesAsScalars.  If you set that to TRUE, the configuration will not use the block sampler and will use univariate samplers instead.

Another question is whether you are providing reasonable initial values.  That could matter.

Also I noticed that you shouldn't really need the CovMat function.  It should work to do

Covm[1:N, 1:N] <- sigmsq * exp(-phi * dist[1:N, 1:N])

directly in the BUGS code.

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 post to this group, send email to nimble...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/cbab2cb9-05f5-4065-9519-f7f9dd182182%40googlegroups.com.

[Lu Zhang]

Hi Perry,

    Thank you so much for your help. My code works pretty well after setting the multivariateNodesAsSaclers to be true in configureMCMC. I attached the zip file for all my code. There are three subfolders in my folder: data, projects, and results. You can generate and get the data in subfolder "data". In the folder "project", there are three files, the one I posted here is "FullGP.R". I also have questions from "nngp_gibbs.R". I may open a new topic about that. The third one (nngp.R) is still underdeveloped. I plan to write a user-defined density function in that one, but I am still confused by how to do that. And I put all simulation results in subfolder "results".

Thank you very much!

Best,

Lu

[nk...@berkeley.edu]

Hello,

We are having a similar problem to the one described by Lu Zhang with the implementation of a spatial hierarchical model using NIMBLE. A minimal reproducible example including data to run the model is attached as the file "minimal_example_KCNK.R".

We have written a two-level hierarchical model with four hyperparameters: alpha, sigma, tau and lambda as follows:

# define hierarchical model with BUGS code

siteHierarchyCode <- nimble::nimbleCode({

  # ---

  # Define priors for hyperparameters

  # ---

  # mean of means

  alpha ~ dnorm(mean = 0,sd = 1000)

  # sd of means : approximates Jeffrey's prior (inverse prior)

  tau ~ dgamma(shape = 0.0001, rate = 0.0001)

  # sd at each site : approximates Jeffrey's prior (inverse prior)

  sigma ~ dgamma(shape = 0.0001, rate = 0.0001)

  # correlation length at each site : approximates Jeffrey's prior (inverse prior)

  lambda ~ dgamma(shape = 0.0001, rate = 0.0001)

  # ---

  # Define hierarchical model at each site

  # ---

  for (i in 1:I){ # for each site

    # the mean is drawn from a univariate gaussian distribution

    mu[i] ~ dnorm(mean = alpha,sd = tau) # distribution of mean at site i

    mat_mu_i[i,1:J_i[i]] <- mu[i] * mat_ones[i,1:J_i[i]]

    # the covariance matrix is defined by an exponential variogram function

    Sigma[1:J_i[i],1:J_i[i],i] <-  sigma^2 * (1 - exp( -(d_ij[1:J_i[i],

                                                              1:J_i[i],

                                                              i]/lambda) ))

    # observations are multivariate normal

    y[i,1:J_i[i]] ~ dmnorm(mat_mu_i[i,1:J_i[i]],

                           cov = Sigma[1:J_i[i],1:J_i[i],i])

  }

})

We run the model using the following code:

#'

#' # Create NIMBLE model and run #

#'

# actually create the model

siteModel <-

  nimble::nimbleModel(code = siteHierarchyCode,

                      constants = siteConst,

                      data = siteData,

                      inits = siteInit,

                      name = 'siteModel')

# compile the model

CsiteModel <- nimble::compileNimble(siteModel)

# create compile and run MCMC

siteModelConf <- nimble::configureMCMC(model = siteModel,

                                       print=T,

                                       multivariateNodesAsScalars=T)

siteModelConf$addMonitors(c('alpha','tau','sigma','lambda'))

siteModelMCMC <-

  nimble::buildMCMC(siteModelConf)

CsiteModelMCMC <-

  nimble::compileNimble(siteModelMCMC, project = siteModel)

CsiteModelMCMC$run(10^5)

While the model complies and runs, the hyperparameters lambda and sigma remain at their initial values, while hyperparameters alpha and tau appear to have been sampled by the MCMC. We have tried modifying the initial values and have so far had no success.

We would like to know, is this problem due to the way the model is defined, or, is there a specific method for running this model such that lambda and sigma are computed?

Thank you,

Karina Cucchi and Nura Kawa

Department of Civil and Environmental Engineering, U.C. Berkeley

[nimble project]

Hi Nura,

Thanks for the inquiry.

I built your model and checked on which parts are validly initialized from your siteInit and resulting calculations or not.  To do this I ran CsiteModel$calculate() and then CsiteModel$getLogProb("y[1, 1:3]") and similar for other nodes.  And I inspected some of the values, e.g. CsiteModel$Sigma.  (I could just as well have used the uncompiled version).  It looks like your Sigma has zeros on the diagonal.  That can't be right is making the log probability values for the y nodes be NaN.  Were you maybe exploring different model formulations and really want sigma^2 * exp(-d_ij[....] / lambda), which would put sigma^2 on the diagonal?

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

To post to this group, send email to nimble...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/98524c1f-bed7-4ee1-9b03-d5b008f6977d%40googlegroups.com.

[nk...@berkeley.edu]

Hello,

Thank you very much for your reply! As you pointed out I had defined sigma incorrectly, and, after changing it, the code now works!

Nura

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

To post to this group, send email to nimble...@googlegroups.com.