Error: checkDistributionFunctions: random generation function...

[Keith Lau]

Hello,

I make up the 4-parameter beta distribution in NIMBLE. However, I always got a message:
checkDistributionFunctions: random generation function `rbeta4` is missing `returnType` or `returnType` does not match the type of the `x` argument to the corresponding density function.

The problem seems stemming from the random generation function. The code is provided below. I believe there is some problem in "rbeta4" below but I cannot spot it.

library("nimble")

dbeta4 <- nimbleFunction(
  run = function(x=double(0,default=0.5), a=double(0,default=1),b=double(0,default=1),lower=double(0,default=0),upper=double(0,default=1), log = logical(0, default = FALSE)) {
    returnType(double(0))
    de = (upper - lower)
    if (de==0){
      p = 0.0
    }else{
      p = dbeta((x - lower)/(upper - lower), a, b)/de
    }
   
    if (log) return(log(p))
    else return(p)
  })
assign('dbeta4', dbeta4, envir = .GlobalEnv)

rbeta4 <- nimbleFunction(
  run = function(n = integer(0), a=double(0,default=1),b=double(0,default=1),lower=double(0,default=0),upper=double(0,default=1)) {
    returnType(double(0))
    return((upper - lower) * rbeta(n, a, b) + lower)
  })
assign('rbeta4', rbeta4, envir = .GlobalEnv)



x = seq(from=0, to=3, by=0.01)

a = 2
b = 2
lower = 1
upper = 2
r = rbeta4(n=10000,a,b,lower,upper)


code <- nimbleCode({
 
  for (n in 1:N){
    r[n] ~ dbeta4(a,b,lower,upper)
  }
  a ~ dunif(0.01,5)
  b ~ dunif(0.01,5)
 
})

constants <- list(N=length(r),lower=lower,upper=upper)

data <- list(r=r)

inits <- list(a=0.01,b=0.01)

monitors = c("a","b")


MODEL <- nimbleModel(code = code, name = "MODEL", constants = constants,data = data, inits = inits)

cmodel <- compileNimble(MODEL)

conf <- configureMCMC(MODEL, monitors = monitors
                      ,useConjugacy=TRUE ,print = FALSE
                      ,enableWAIC=TRUE
                      ,onlySlice=FALSE
)

modelMCMC <- buildMCMC(conf)
cmodelMCMC <- compileNimble(modelMCMC, project = MODEL ,resetFunctions = TRUE)

burnin = 1000
keep = 10000
nchains = 1
rand.seed = 1

niter = burnin + keep

out <- system.time(samples <- runMCMC(cmodelMCMC, niter = niter, nburnin = burnin, nchains=nchains
                                              ,samplesAsCodaMCMC=TRUE
                                              ,samples= TRUE
                                              ,summary=TRUE
                                              ,WAIC = TRUE
                                              ,setSeed = rand.seed))

[Daniel Turek]

Keith, nice work with the rbeta4 distribution.  Overall it looks great, there's one small thing which is giving nimble trouble - which arguably should be ok, and is something I'll file an issue for us to look into.

Basically, when you define the dbeta4 density function, you specify the type of the random variable as the "x" argument, in this case a scalar value defined by "x=double(0)" in the list of arguments.  Here, when you define this, you cannot also give a default value for "x", as you had.  So the declaration needs to look like:

    run = function(x=double(0), a=double(0,default=1),b=double(0,default=1),lower=double(0,default=0),upper=double(0,default=1), log = logical(0, default = FALSE)) {

Is this a little heavy-handed of nimble's type checking?  Should default values for the "x" argument of the density function be allowed?  These are deeper questions.  But for now, changing that first line (specifically, to remove "default=0.5") should get you going.

Like I said, I'll file an issue for the development team to think more about this.

Cheers,

Daniel

--
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/4090e981-a45a-4e3b-a56b-e6e521af34fbn%40googlegroups.com.

[Keith Lau]

Thank you for the reply, Daniel. I have changed "x=double(0,default=0.5)" to "x=double(0)". It gave a message:

Number of dimensions 1 of the return() argument does not match number 0 given in the returnType() statement. This occurred for: return((ARG5_upper_ - ARG4_lower_) * rbeta(n=ARG1_n_,shape1=ARG2_a_,shape2=ARG3_b_) + ARG4_lower_) This was part of the call: { return((ARG5_upper_ - ARG4_lower_) * rbeta(n=ARG1_n_,shape1=ARG2_a_,shape2=ARG3_b_) + ARG4_lower_)

I struggled it for a long while. I found a solution is to change rbeta4:

rbeta4 <- nimbleFunction(
  run = function(n = integer(0), a=double(0),b=double(0),lower=double(0),upper=double(0)) {
    returnType(double(0))
    x <- (upper - lower) * rbeta(n, a, b) + lower
    return(x[1])
  })

where I must to specify "return(x[1]) " that only return a scalar to make NIMBLE work. I think this approach is awkward because rbeta4 is not vectorized. There may be a better workaround but I still am not sure how to keep rbeta4 vectorized.

Daniel Turek 在 2023年10月8日 星期日晚上7:38:06 [UTC+8] 的信中寫道：

[Daniel Turek]

Keith, can you change the rbeta(n, a, b) inside your rbeta4 function to instead be rbeta(1, a, b)?  I should also have suggested this in my first email.

Under the hood, nimble *always* uses n=1 for calculations in the model.  I know you used the "n" argument for generating many realizations of your distribution in your script, on the line:

r = rbeta4(n=10000,a,b,lower,upper)

But the nimble compiler (type-checking) will be happier if you use "rbeta(1, a, b)" in your function code, so that it knows (always, for certain) that the return value "x" is a scalar.  Then, you should not need to say "x[1]", as you found works.

In your script, then, you will also need to change 

r = rbeta4(n=10000,a,b,lower,upper)

to instead be something like:

r = replicate(10000, rbeta4(n=1,a,b,lower,upper))

Please let me know if these suggestions help.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5d7fa823-c802-413e-9df3-5c577aba3fe5n%40googlegroups.com.

[Keith Lau]

Thank you again, Daniel. Your suggestions all work well. It is a great lesson to know more how NIMBLE works under the hood!

Daniel Turek 在 2023年10月8日 星期日晚上8:21:59 [UTC+8] 的信中寫道：