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))
###--
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.
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]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:
Any thoughts?
All Best,To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/541639e9-61f1-6d2c-d1af-96ff9d54b43f%40gmail.com.
Hi Perry,
Silly me! That was a stupid mistake!!! Thank you so much and
thanks for the tip for future reference!