HMC with user-defined distribution

[Miguel Ángel Martínez Beneito]

Hi all,

I am having some problems when trying to use HMC with a distribution defined by myself. The distribution is the Leroux et al. CAR distribution which I have programmed (https://onlinelibrary.wiley.com/doi/10.1002/sim.10166) as:

  dcar_leroux <- nimbleFunction(
    name = 'dcar_leroux',
    run = function(x = double(1),        # Spatial random effect (vector)
                 rho = double(0),      # Amount of spatial dependence (scalar)
                 sd.theta = double(0), # Standard deviation (scalar)
                 Lambda = double(1),   # Eigenvalues of matrix D - W
                 from.to = double(2),  # Matrix of distinct pairs of neighbors from.to[, 1] < from.to[, 2]
                 log = integer(0, default = 0)) {
      #returnType(double(0))
      # Number of small areas
      NMuni <- dim(x)[1]
      # Number of distinct pairs of neighbors
      NDist <- dim(from.to)[1]
      # Required vectors
      x.from <- nimNumeric(NDist)
      x.to <- nimNumeric(NDist)
      for (Dist in 1:NDist) {
        x.from[Dist] <- x[from.to[Dist, 1]]
        x.to[Dist] <- x[from.to[Dist, 2]]
      }
   
      logDens <- sum(dnorm(x[1:NMuni], mean = 0, sd = sd.theta * pow(1 - rho, -1/2), log = TRUE)) -
        NMuni/2 * log(1 - rho) +  sum(log(rho * (Lambda[1:NMuni] - 1) + 1))/2 -
        pow(sd.theta, -2) * rho * sum(pow(x.from[1:NDist] - x.to[1:NDist], 2))/2
      if(log) return(logDens)
      else return(exp(logDens))
      returnType(double())
    },
    buildDerivs= list(run = list(ignore = 'Dist'))
  )

  rcar_leroux <- nimbleFunction(
    name = 'rcar_leroux',
    run = function(n = integer(0),
                   rho = double(0),
                   sd.theta = double(0),
                   Lambda = double(1),
                   from.to = double(2)) {
      returnType(double(1))
      nimStop("user-defined distribution dcar_leroux provided without random generation function.")
      x <- nimNumeric(length(Lambda))
      return(x)
    }
  )

  assign('dcar_leroux', dcar_leroux, envir = .GlobalEnv)
  assign('rcar_leroux', rcar_leroux, envir = .GlobalEnv)

The problem comes when I try to build the corresponding model with

> buildMCMC(configuredModel)

as this returns the error:

> Error: parameterTransform doesn't handle 'dcar_leroux' distributions

It seems as if the parameterTransform function, which is called by buildMCMC, does not consider the case of having this, or any, user-defined distribution. Can this be solved in some way?

I do not attach any reproducible case study since the model that I am working on is quite intrincated and hardly illustrative of the problem mentioned. However, if necessary, I could build a simpler version of the model for illustrative purposes.

Thanks in advance.

[Perry de Valpine]

Dear Miguel,

Thanks for the very clear report and example. Indeed, the parameter transformation system handles user-defined scalar distributions (by making use of any provided bounds), but it does not handle user-defined multivariate distributions. That is an important need that should be filled. Thanks to Daniel we can provide a branch of nimble with a very small modification that should work for your case. Can you please try installing from

remotes::install_github("nimble-dev/nimble", ref = "transform_user_dist", subdir = "packages/nimble")

and let us know if that works.

Note that it still doesn't handle parameter transformation for user-provided multivariate distributions in any general way. All it does is default to using no transformation (i.e. identity transformation) rather than throwing an error. It will also give a friendly warning about the situation.

HTH

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 visit https://groups.google.com/d/msgid/nimble-users/ec460af0-a86e-4bb1-a1a8-0d2ed77e7e38n%40googlegroups.com.

[Miguel Ángel Martínez Beneito]

Perry (and Daniel),

First of all, thank you very much for your assistance. 

I have installed the suggested version of Nimble. As a result, now the buildMCMC function returns, as expected, several repeated warnings ([Warning] `parameterTransform` system detected multivariate user-defined distribution dcar_leroux. No transformation will be applied to any dimension of 'x' values of dcar_leroux. If some values of 'x' are not valid in dcar_leroux, you may encounter errors.  This warning can be disabled using nimbleOptions(parameterTransformWarnUserDists = FALSE)). However, an additional error makes buildMCMC to crash now. The error reports the following:

One or more errors occurred looking for variables in a model (first 10 shown below).
These messages may be cryptic, but generally the variable or expression somewhere in each message was not valid in a model:
Error in 1 %.s% 73 : could not find function "%.s%"
Error in 1 %.s% 73 : could not find function "%.s%"
Error in 1 %.s% 73 : could not find function "%.s%"

...

By the way, 73 is the length of the vector following the Leroux distribution in my model.

Hope this helps ...

Best wishes (and thanks again).

Miguel.

[Perry de Valpine]

Hi Miguel,

This is likely a completely unrelated error. The "%.s%" is used in some internal processing when setting up the DAG for a model. If it ever appears as in the error message you're seeing, it means something is wrong. Unfortunately, there is a bug where particular ways of using subsets of non-scalar stochastic nodes can cause this error. It is pretty unusual, but then again if it's arising from a model you need, it might not seem unusual at all. There was a thread on the topic in November. For now, the hope is that once you see the kind of situation that can create the problem, you might be able to change how you write your model to work around the problem. To diagnose further or make specific suggestions for you, I think we'd need to see your model code with at least the constants so that we can build the model (i.e. call nimbleModel) and inspect it. I could be wrong about the cause, but it does sound like the same issue as the November thread I linked.

HTH

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/6edcb9d7-fd09-4815-be47-56a48d2a6b2fn%40googlegroups.com.

[Miguel Ángel Martínez Beneito]

Hi Pierre,

Sorry for the late reply but I have several problems over my desk right now.

I just wanted to let you know that the github development version that you suggested seems to fix the issue with HMC and our Leroux distribution, therefore I would appreciate so much if this was included in future versions of Nimble. At the end, I have tried to fit a simple model with Leroux distributed random effects with HMC and it works. Thus, as you suggested, the "%.s%" problem would be related to additional particular issues of my more complex version of the model (curiously this is a problem just for HMC but not for regular MCMC). In fact, my complex model has an array structure with inner (Leroux) multivariate nodes as the thread on the topic in November that you pointed out so it definitely seems a related issue. If I figure out something else about this issue I will let you know.

Finally, just a minor comment, since the installation of the provisional github version of nimble, Rstudio complains ("Error in fetch(key): lazy-load database '/home/.../nimble.rdb' is corrupt") every time that I try to access Nimble help either directly or indirectly (Rstudio tries to do it while writting code). This is not important since everything else seems to work fine, but at the end it results a bit annoying. I am just mentioning this just in case you should do something in this regard in future versions of Nimble that include the Leroux distribution fix just implemented.

Best wishes.

[Miguel Ángel Martínez Beneito]

Hi,

I am retaking this issue. 

A time ago I defined a multivariate distribution to be used in my models, with HMC if possible. A problem arised (Error: parameterTransform doesn't handle 'dcar_leroux' distribution) since, as mentioned, the parameterTransform function is not intended to work with multivariate distributions. As a (succesful) solution you uploaded a nimble patch to Github that skipped the parameterTransform function for my distribuion (dcar_Leroux at that time). With this patch I made it work and everything went smoothly. Now I am having a similar problem with a new distribution that we are developin, "dSPDE" in this case. I wonder, would it be possible to develop a similar patch for HMC to correctly deal with this distribution?

All the best.

[Perry de Valpine]

Hi Miguel,

There is a change that will be part of an upcoming release (should be fairly soon) that makes the default behavior in this case be to apply an identity transformation and issue a warning, which can be turned off by changing a nimbleOption. I think this change was made in recognition of your previous issue!

Please try installing from our development branch, "devel", and see if that works.

remotes::install_github("nimble-dev/nimble", ref = "devel", subdir = "packages/nimble")

Thanks

Perry

To view this discussion visit https://groups.google.com/d/msgid/nimble-users/ec860be6-ba86-41f7-ad6f-d29c872051fcn%40googlegroups.com.