run_discrete | R Documentation |
A wrapper function to conduct Markov Chain Monte Carlo (MCMC) sampling using 'nimble'.
run_discrete( model, data, constants, inits, params, dimensions = NULL, nimFun = NULL, niter = 1000, nburnin = 100, thin = 1, nchains = 1, parallel = FALSE, RNGseed = NULL )
model |
The |
data |
A list of data inputs needed to run SCR models in |
constants |
A list of constants needed to run SCR models in |
inits |
Starting values for stochastic parameters to begin MCMC sampling. |
params |
A vector of character strings that define the parameters to trace in the MCMC simulation. |
dimensions |
A named list of dimensions for variables. Only needed for variables used with empty indices in model code that are not provided in constants or data. |
nimFun |
A list of nimbleFunction(s) to be used in the model when using
|
niter |
An integer value of the total number of MCMC iterations to run per chain. |
nburnin |
An integer value of the number of MCMC iterations to discard as 'burnin'. |
thin |
An integer value of the amount of thinning of the chain. For
example, |
nchains |
An integer value for the number of MCMC chains to run |
parallel |
A logical value indicating whether MCMC chains shoud be run
in parallel processing. Default is |
RNGseed |
An integer value specifying the random number generating seed
using in parallel processing. Also used when |
This function provides a wrapper to easily run Bayesian SCR models
using nimble
.
A list of MCMC samples for each parameter traced with length equal to the number of chains run.
Daniel Eacker
## Not run: # simulate a single trap array with random positional noise x <- seq(-800, 800, length.out = 5) y <- seq(-800, 800, length.out = 5) traps <- as.matrix(expand.grid(x = x, y = y)) set.seed(200) traps <- traps + runif(prod(dim(traps)),-20,20) mysigma = c(220, 300) # simulate sex-specific scaling parameter mycrs = 32608 # EPSG for WGS 84 / UTM zone 8N pixelWidth = 100 # store pixelWidth or grid resolution # create state-space Grid = grid_classic(X = traps, crs_ = mycrs, buff = 3*mysigma, res = pixelWidth) # create polygon for mask library(sf) poly = st_sfc(st_polygon(x=list(matrix(c(-1765,-1765,1730,-1650,1600,1650,0,1350, -800,1700,-1850,1000,-1765,-1765),ncol=2, byrow=TRUE))), crs = mycrs) # create habitat mask hab_mask = mask_polygon(poly = poly, grid = Grid$grid, crs_ = mycrs, prev_mask = NULL) # simulate data for uniform state-space and habitat mask data3d = sim_classic(X = traps, ext = Grid$ext, crs_ = mycrs, sigma_ = mysigma, prop_sex = 0.6, N = 200, K = 4, base_encounter = 0.15, enc_dist = "binomial", hab_mask = hab_mask, setSeed = 100) # total augmented population size M = 400 # get initial activity center starting values s.st = initialize_classic(y=data3d$y, M=M, X=traps, ext=Grid$ext, hab_mask = hab_mask) # convert traps and starting locations to discrete format d_list = discretize_classic(X=traps, grid = Grid$grid, s.st = s.st, crs_ = mycrs, hab_mask = hab_mask) # inspect discrete data list str(d_list) # prepare data data = list(y=data3d$y) data$y = data$y[which(apply(data$y, 1, sum)!=0),,] # remove augmented records # covert to 2d by summing over individuals and traps data$y = apply(data$y, c(1,2), sum) # add discretized traps data$X = d_list$X100 # convert units to km # add grid data$grid = d_list$grid/1000 # convert units to km # prepare constants (note get density in activity center/100 m2 rather # than activity centers/m2) constants = list(M = M,n0 = nrow(data$y),J=dim(data$y)[2], K=dim(data3d$y)[3], nPix=d_list$nPix,pixArea = 0.1^2, sigma_upper = 1, A = (sum(hab_mask)*(pixelWidth/100)^2)) # augment sex data$sex = c(data3d$sex,rep(NA,constants$M-length(data3d$sex))) # add z and zeros vector data for latent inclusion indicator data$z = c(rep(1,constants$n0),rep(NA,constants$M - constants$n0)) data$zeros = c(rep(NA,constants$n0),rep(0,constants$M - constants$n0)) # define all initial values inits = list(sigma = runif(2, 0.250, 0.350), s = d_list$s.st,alpha0 = 3, p0 = runif(1, 0.05, 0.15),psi_sex = runif(1,0.5, 0.7), sex=c(rep(NA,constants$n0),rbinom(constants$M-constants$n0,1,0.6)), z=c(rep(NA,constants$n0),rep(0,constants$M-constants$n0))) # parameters to monitor params = c("sigma","psi","p0","N","D","EN","psi_sex","alpha0","s","z") # get model discrete_model = get_discrete(type="marked",dim_y = 2, enc_dist = "binomial", sex_sigma = TRUE,trapsClustered=FALSE) # run model library(tictoc) tic() # track time elapsed out = run_discrete(model = discrete_model, data=data, constants=constants, inits=inits, params = params,niter = 10000, nburnin=1000, thin=1, nchains=2, parallel=TRUE, RNGseed = 500) toc() nimSummary(out, exclude = c("s","z")) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.