Creating needed NIMBLE functions
#Binomial coefficient calculator
nlchoose <- nimbleFunction(
run = function(n = double(0), k = double(0)){
return(sum(log((n-k+1):n)) - sum(log(1:k)))
}
)
#discrete uniform density function
dunifd <- nimbleFunction(
run = function(x = double(0), lower = double(0), upper = double(0), log = integer(0, default = 1)){
returnType(double(0))
logProb <- -log(upper-lower+1)
if(log) return(logProb)
else return(exp(logProb))
})
#discrete uniform RNG
runifd <- nimbleFunction(
run = function(n = integer(0, default = 1), lower = double(0), upper = double(0)){
returnType(double(0))
if(n!=1) print('rhyperg only allows n=1; using n = 1')
dev <- runif(1,0,1)
unitd <- 1/(upper-lower+1)
return(floor(lower + dev/unitd))
}
)
#hypergeometric density function
dhyperg <- nimbleFunction(
run = function(x = double(0), nx = double(0), N1 = double(0), N0 = double(0),
log = integer(0, default = 1)){
returnType(double(0))
logProb <- nlchoose(N1,x) + nlchoose(N0,nx-x) - nlchoose(N0+N1,nx)
if(log) return(logProb)
else return(exp(logProb))
})
#hypergeometric RNG
rhyperg <- nimbleFunction(
run = function(n = integer(0, default = 1), nx = double(0), N1 = double(0), N0 = double(0)){
returnType(double(0))
if(n!=1) print('rhyperg only allows n=1; using n = 1')
dev <- runif(1,0,1)
minx <- max(nx - N0, 0)
maxx <- min(nx, N1)
xs <- numeric(length = maxx-minx+1)
ps <- xs
xs[1] <- minx
ps[1] <- dhyperg(minx, nx, N1, N0, log = 0)
for(i in (minx+1):maxx){
xs[i-minx+1] <- i
ps[i-minx+1] <- dhyperg(i, nx, N1, N0, log = 0) + ps[i-minx]
}
return(max(xs[dev > ps]))
}
)
#Have to register Nimble functions (at least distributions) to R's global environment
assign('nlchoose', nlchoose, .GlobalEnv)
assign('dunifd', dunifd, .GlobalEnv)
assign('runifd', runifd, .GlobalEnv)
assign('dhyperg', dhyperg, .GlobalEnv)
assign('rhyperg', rhyperg, .GlobalEnv)Coding NIMBLE Model
A1_basic_nimb <-
nimbleCode({
for(v in 1:NSite)
for(t in 1:NYr){
S1_ct[v,t] ~ dbinom(Q[v,t], p_S1)
S3_ct[v,t] ~ dhyperg(n_S3[v,t], Q[v,t], P[v,t] - Q[v,t])
Q[v,t] ~ dunifd(0, P[v,t])
}
p_S1 ~ dunif(0,1)
})
#print(A1_basic_nimb)Defining Constants and Data
A1bConst <- list('NYr'=4, 'NSite'=6)
A1bData <- list('S1_ct' = matrix(c(2,5,3,1,2,0,1,10,4,0,0,2,3,2,6,0,0,1,1,1,1,0,0,0),nrow=6),
'S3_ct' = matrix(c(4,9,1,3,7,2,12,22,4,2,5,1,7,4,10,6,10,1,10,6,6,3,7,11),nrow=6),
'P' = matrix(c(2347,4465,2475,2708,4586,3442,1513,4412,2595,3129,3727,3449,1720,4401,3826,3328,4139,2673,1700,4432,1713,3022,2376,2552),nrow=6),
'n_S3' = matrix(c(1200,1266,1200,1205,1200,1202,1200,1200,1202,1200,1200,1213,1292,1200,1230,1218,1200,1297,1191,1176,1200,1209,1200,1277),nrow=6))
print(A1bData)## $S1_ct
## [,1] [,2] [,3] [,4]
## [1,] 2 1 3 1
## [2,] 5 10 2 1
## [3,] 3 4 6 1
## [4,] 1 0 0 0
## [5,] 2 0 0 0
## [6,] 0 2 1 0
##
## $S3_ct
## [,1] [,2] [,3] [,4]
## [1,] 4 12 7 10
## [2,] 9 22 4 6
## [3,] 1 4 10 6
## [4,] 3 2 6 3
## [5,] 7 5 10 7
## [6,] 2 1 1 11
##
## $P
## [,1] [,2] [,3] [,4]
## [1,] 2347 1513 1720 1700
## [2,] 4465 4412 4401 4432
## [3,] 2475 2595 3826 1713
## [4,] 2708 3129 3328 3022
## [5,] 4586 3727 4139 2376
## [6,] 3442 3449 2673 2552
##
## $n_S3
## [,1] [,2] [,3] [,4]
## [1,] 1200 1200 1292 1191
## [2,] 1266 1200 1200 1176
## [3,] 1200 1202 1230 1200
## [4,] 1205 1200 1218 1209
## [5,] 1200 1200 1200 1200
## [6,] 1202 1213 1297 1277A1bInit <- list(Q = matrix(rep(10,A1bConst$NYr*A1bConst$NSite),nrow=A1bConst$NSite),p_S1=0.5)
A1b <-try( nimbleModel(code=A1_basic_nimb, name='A1 basic version', constants = A1bConst, data=A1bData,
inits = A1bInit,debug=F))## defining model...## Registering the following user-provided distributions: dhyperg .
## NIMBLE has registered dhyperg as a distribution based on its use in BUGS code. Note that if you make changes to the nimbleFunctions for the distribution, you must call 'deregisterDistributions' before using the distribution in BUGS code for those changes to take effect.
## Registering the following user-provided distributions: dunifd .
## NIMBLE has registered dunifd as a distribution based on its use in BUGS code. Note that if you make changes to the nimbleFunctions for the distribution, you must call 'deregisterDistributions' before using the distribution in BUGS code for those changes to take effect.## building model...print(A1b)## [1] "Error in nndf_generateGetParamSwitchFunction(allParams[paramNamesToUse], : \n problem creating switch function for getParam from value=Q[getNodeFunctionIndexedInfo(INDEXEDNODEINFO_, 1), getNodeFunctionIndexedInfo(INDEXEDNODEINFO_, 2)],NA=NULL,NA=NULL\n"
## attr(,"class")
## [1] "try-error"
## attr(,"condition")
## <simpleError in nndf_generateGetParamSwitchFunction(allParams[paramNamesToUse], paramIDs[paramNamesToUse], type = "double", nDim = nDimSupported): problem creating switch function for getParam from value=Q[getNodeFunctionIndexedInfo(INDEXEDNODEINFO_, 1), getNodeFunctionIndexedInfo(INDEXEDNODEINFO_, 2)],NA=NULL,NA=NULL>