bad parameters for distribution ... (No available re-parameterization found.) for a custom model

[Enrico Ryunosuke Crema]

Dear List,

I’m trying to extend some ideas introduced in this paper using a custom model in Nimble but I'm stuck with the error message “Error: bad parameters for distribution dAoristicExponentialGrowth(a = 4, b = 1, r = r). (No available re-parameterization found.)”

The code to reproduce my error is pasted below, but briefly, the idea is that I am trying to fit a statistical model of the time-frequency of my data. In the case of dAoristicExponentialGrowth, that will increase with an exponential rate r between points in time a and b (notice that a>b as time is expressed in numbers of years from present). Time is represented as a discrete distribution, so effectively dAoristicExponentialGrowth is just a wrapper for a multinomial distribution. The problem is that my observations are characterised by some uncertainty, and this is expressed as vector of probabilities for each dscrete time unit. It is worth noting that: 1) the vector has the same length as the interval between a and b in dAoristicExponentialGrowth, and 2) the uncertainty cannot be expressed by using some statistical distribution. The way I approached this was to create the function dAoristicExponentialGrowth, which would expect as an observation a vector of probabilities. The function seems to work fine, i.e. it does generate the correct log-likelihood. However when I use nimbleModel I get the error above. Not sure if this is something that can be done in nimble at all or if I am doing something wrong entirely. Any help would be greatly appreciated!

All Best,

Enrico

#### Code ###

library(nimble)

# Define custom model

dAoristicExponentialGrowth=nimbleFunction(

  run = function(x = double(1),a=double(0),b=double(0),r=double(0), log = integer(0)) {

    returnType(double(0))

    t = 1:(abs(b-a)+1)

    n = numeric(abs(b-a))

    tfinal = abs(b-a)+1

    for (i in 1:tfinal)

    {

      n[i] = (1+r)^t[i]

    }

    p = n/sum(n)

    pg = x * p

    logProb = log(sum(pg))

    if(log) {

      return(logProb)

    } else {

      return(exp(logProb))

    }

  })   

# register model 

registerDistributions(list(

  dAoristicExponentialGrowth = list(

    BUGSdist = "dAoristicExponentialGrowth(x,a,b,r)",

    Rdist = "dAoristicExponentialGrowth(x,a,b,r)",

    pqAvail = FALSE

  )))

# create some mock data

d  <- list()

d$theta  <- rbind(c(0,0,0.5,0.5),c(0,1,0,0))

constants <- list()

constants$N  <- 2

model <- nimbleCode({

for (i in 1:N)

{

theta[i,] ~ dAoristicExponentialGrowth(a=4,b=1,r=r)

}

r ~ dexp(1)

})

model <- nimbleModel(model, constants = constants, data = d, inits = list(r=0.001))

###

[Sally paganin]

Hi Enrico, 

I haven't fully thought about the details of the distribution you are implementing, but I have a few comments/suggestions on the code.

Since you are just defining a density, you can skip the call to registerDistributions, and let nimble automatically register your distribution when your call nimbleModel. I tested your code, and it runs.

There is more about user distributions and cases where to use registerDistributions in Ch 12.2 of the user manual, but a couple of things I noticed in your use of registerDistributions:

• for the BUGSdist argument, one needs to provide the name of the distribution followed by only the parameters names in parenthesis; so it would be "dAoristicExponentialGrowth(a,b,r)" omitting "x"
• the argument Rdist is for providing alternative parameterization of the distribution and is optional, so you would not need that in this case. 

Hope this helps!

Best,

Sally

--
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/e3c97e30-5bc2-4421-83cb-9a8f1d9ab12cn%40googlegroups.com.

[Enrico R Crema]

Hi Sally,

Many thanks for the swift reply. Not registering the distribution solved the issue with nimbleModel(), but then running the following:

model <- nimbleModel(model, constants = constants, data = d, inits = list(r=0.001))

cModel <- compileNimble(model)
conf <- configureMCMC(model)
MCMC <- buildMCMC(conf)
cMCMC <- compileNimble(MCMC, project = cModel)
samples <- runMCMC(cMCMC, niter = 5000, thin=1,nburnin = 1000,samplesAsCodaMCMC = T)

Crashed the R session with the following message:

Run-time size error: expected dim(ARG_x_)[1] == dim(p)[1]
Run-time size error: expected dim(ARG_x_)[1] == dim(p)[1]
munmap_chunk(): invalid pointer

I tried then registering the distribution as follows:

registerDistributions(list(
  dAoristicExponentialGrowth = list(
    BUGSdist = "dAoristicExponentialGrowth(a,b,r)",
    pqAvail = FALSE,
    types = c('value = double(1)', 'a = double(0)','b = double(0)','r=double(0)')
  )))

But I still get a similar error message:

Run-time size error: expected dim(ARG_x_)[1] == dim(p)[1]
Run-time size error: expected dim(ARG_x_)[1] == dim(p)[1]
malloc(): invalid size (unsorted)

Any thoughts?

All Best,
Enrico

[Perry de Valpine]

Hi Enrico,

I think the problem is fixed by:

n = numeric(abs(b-a)+1)

If I follow, you want n to be the same length as t.  On my machine, it compiles and runs and seems to mix with that change.

The error message is not very clear (it is difficult to generate great error-messages from within generated C++ code), but it is basically trying to say that some sizes are not commensurate in some operation. In this case I believe it is `pg = x * p` where it would determine that x and p have different lengths.

In R, `x[i] <- value` will automatically extend the length of x if necessary, but that will not happen in compiled nimbleFunction code.

FWIW, the approach I used to find the problem was to run the MCMC uncompiled, with debugging in your dist:

debug(dAoristicExponentialGrowth)

MCMC$run(niter = 5)

Then I could step through and see the issue by inspection.

Thanks for posting and HTH!

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/541639e9-61f1-6d2c-d1af-96ff9d54b43f%40gmail.com.

[Enrico R Crema]

Hi Perry,

Silly me! That was a stupid mistake!!! Thank you so much and thanks for the tip for future reference!

All Best,
Enrico