if (FALSE) {
## generate random trap locations
nTraps <- 200
traps_xmin <- 0
traps_ymin <- 0
traps_xmax <- 100
traps_ymax <- 200
set.seed(0)
traps_xCoords <- round(runif(nTraps, traps_xmin, traps_xmax))
traps_yCoords <- round(runif(nTraps, traps_ymin, traps_ymax))
trap_coords <- cbind(traps_xCoords, traps_yCoords)
## buffer distance surrounding sides of rectangular discretization grid
## which overlays trap locations
buffer <- 10
## resolution of rectangular discretization grid
resolution <- 10
## creates grid and makeID function,
## for grid overlaying trap locations,
## and to lookup nearest grid cell to any AC
makeGridReturn <- makeGrid(xmin = traps_xmin, xmax = traps_xmax,
ymin = traps_ymin, ymax = traps_ymax,
buffer = buffer,
resolution = resolution)
grid <- makeGridReturn$grid
makeID <- makeGridReturn$makeID
## maximum radis within an individual AC to perform trap calculations,
dmax <- 30
## n = localTraps[i,1] gives the number of local traps
## localTraps[i, 2:(n+1)] gives the indices of the local traps
localTraps <- findLocalTraps(grid, trap_coords, dmax)
plotTraps <- function(i, grid, trap_coords, localTraps) {
plot(grid[,1], grid[,2], pch = '.', cex=2)
points(trap_coords[,1], trap_coords[,2], pch=20, col='forestgreen', cex=1)
if(!missing(i)) {
i <- max(i %% dim(grid)[1], 1)
n <- localTraps[i,1]
trapInd <- numeric(0)
if(n > 0) trapInd <- localTraps[i,2:(n+1)]
theseTraps <- trap_coords[trapInd,, drop = FALSE]
points(theseTraps[,1], theseTraps[,2], pch = 20, col = 'red', cex=1.5)
points(grid[i,1], grid[i,2], pch = 'x', col = 'blue', cex=3)
}
}
## visualise some local traps
plotTraps(10, grid, trap_coords, localTraps)
plotTraps(200, grid, trap_coords, localTraps)
plotTraps(380, grid, trap_coords, localTraps)
## example model code
## using local trap calculations
code <- nimbleCode({
sigma ~ dunif(0, 100)
p0 ~ dunif(0, 1)
for(i in 1:N) {
S[i,1] ~ dunif(0, xmax)
S[i,2] ~ dunif(0, ymax)
Sdiscrete[i,1] <- round(S[i,1]/res) * res
Sdiscrete[i,2] <- round(S[i,2]/res) * res
id[i] <- makeID( Sdiscrete[i,1:2] )
nLocalTraps[i] <- getNumLocalTraps(id[i], localTraps[1:LTD1,1], LTD1)
localTrapIndices[i,1:maxTraps] <-
getLocalTrapIndices(maxTraps, localTraps[1:LTD1,1:LTD2], nLocalTraps[i], id[i])
d[i, 1:maxTraps] <- calcLocalTrapDists(
maxTraps, nLocalTraps[i], localTrapIndices[i,1:maxTraps],
S[i,1:2], trap_coords[1:nTraps,1:2])
g[i, 1:nTraps] <- calcLocalTrapExposure(
nTraps, nLocalTraps[i], d[i,1:maxTraps], localTrapIndices[i,1:maxTraps], sigma, p0)
y[i, 1:nTraps] ~ dbinom_vector(prob = g[i,1:nTraps], size = trials[1:nTraps])
}
})
## generate random detection data; completely random
N <- 100
set.seed(0)
y <- array(rbinom(N*nTraps, size=1, prob=0.8), c(N, nTraps))
## generate AC location initial values
Sinit <- cbind(runif(N, traps_xmin, traps_xmax),
runif(N, traps_ymin, traps_ymax))
constants <- list(N = N,
nTraps = nTraps,
trap_coords = trap_coords,
xmax = traps_xmax,
ymax = traps_ymax,
res = resolution,
localTraps = localTraps,
LTD1 = dim(localTraps)[1],
LTD2 = dim(localTraps)[2],
maxTraps = dim(localTraps)[2] - 1)
data <- list(y = y, trials = rep(1,nTraps))
inits <- list(sigma = 1,
p0 = 0.5,
S = Sinit)
## create NIMBLE model object
Rmodel <- nimbleModel(code, constants, data, inits,
calculate = FALSE, check = FALSE)
## use model object for MCMC, etc.
}
Run the code above in your browser using DataLab