diff --git a/.Rbuildignore b/.Rbuildignore index 506d650..7ac7f87 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -6,7 +6,9 @@ ^\.github$ ^Dockerfile$ ^_pkgdown\.yml$ +^\.devcontainer$ ^docs$ ^pkgdown$ ^CNAME$ ^air\.toml$ +^src/StanExports_* \ No newline at end of file diff --git a/.gitignore b/.gitignore index 5ec3713..22e2954 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,7 @@ data-raw .DS_Store docs .vscode +src/RcppExports.cpp +src/stanExports_* +*.o +*.so \ No newline at end of file diff --git a/R/COA_Tag_Integrated.R b/R/COA_Tag_Integrated.R index bdf9683..f807091 100644 --- a/R/COA_Tag_Integrated.R +++ b/R/COA_Tag_Integrated.R @@ -2,20 +2,22 @@ #' Fits a test-tag integrated Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan #' -#' @param nind Number of tagged individuals -#' @param nrec Number of receivers -#' @param ntime Number of time steps -#' @param ntest Number of test tags -#' @param ntrans Number of expected transmissions per tag per time interval -#' @param y Array of detection data, where row = individual, column = time step, and matrix = receiver -#' @param test Array of test tag detection data, where row = individual tag, column = time step, and matrix = receiver -#' @param recX Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency) -#' @param recY Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency) -#' @param xlim East-west boundaries of spatial extent (receiver array + buffer) -#' @param ylim North-south boundaries of spatial extent (receiver array + buffer) -#' @param testX Test tag coordinates in the east-west direction (should be projected and scaled for computational efficiency) -#' @param testY Test tag coordinates in the north-south direction (should be projected and scaled for computational efficiency) -#' @param ndraws to be passed to `generated_quantities`. Changes the number of draws. Default is 10. +#' @param n_ind Number of tagged individuals +#' @param n_rec Number of receivers +#' @param n_time Number of time steps +#' @param n_test Number of test tags +#' @param n_trans Number of expected transmissions per tag per time interval +#' @param n_trans_test Number of expected transmissions per sentinel tag per time interval +#' @param det Array of detection data, where row = individual, column = receiver, and matrix = time step +#' @param det_test Array of test tag detection data, where row = individual tag, column = receiver, and matrix = time step +#' @param rec_x Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency) +#' @param rec_y Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency) +#' @param x_lim East-west boundaries of spatial extent (receiver array + buffer) +#' @param y_lim North-south boundaries of spatial extent (receiver array + buffer) +#' @param test_x Test tag coordinates in the east-west direction (should be projected and scaled for computational efficiency) +#' @param test_y Test tag coordinates in the north-south direction (should be projected and scaled for computational efficiency) +#' @param decay desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian". +#' @param n_draws to be passed to `generated_quantities`. Changes the number of draws. Default is 10. #' @param ... Additional arguments passed to `sampling` from `rstan`. #' This can include setting `chains`, `iter`, `warmup`, and `control`. Please see #' `rstan::sampling()` for more info. @@ -26,42 +28,46 @@ #' #' @seealso [rstan::sampling()] #' @export + COA_TagInt <- function( - nind, - nrec, - ntime, - ntest, - ntrans, - y, - test, - recX, - recY, - xlim, - ylim, - testX, - testY, - ndraws = NULL, + n_ind, + n_rec, + n_time, + n_test, + n_trans, + n_trans_test, + det, + det_test, + rec_x, + rec_y, + x_lim, + y_lim, + test_x, + test_y, + decay = "gaussian", + n_draws = NULL, ... ) { # First move everything into a list standata <- list( - nind = nind, - nrec = nrec, - ntime = ntime, - ntest = ntest, - ntrans = ntrans, - y = y, - test = test, - recX = recX, - recY = recY, - xlim = xlim, - ylim = ylim, - testX = testX, - testY = testY + n_ind = n_ind, + n_rec = n_rec, + n_time = n_time, + n_test = n_test, + n_trans = n_trans, + n_trans_test = n_trans_test, + det = det, + det_test = det_test, + rec_x = rec_x, + rec_y = rec_y, + x_lim = x_lim, + y_lim = y_lim, + test_x = test_x, + test_y = test_y ) # validate this list prior to sending it to the model - exp_len <- expected_lengths(recX = recX, recY = recY, ntest_len = ntest) + exp_len <- expected_lengths(rec_x = rec_x, rec_y = rec_y, n_test_len = n_test) validate_standata(standata, exp_len) # set rstan options @@ -70,53 +76,68 @@ COA_TagInt <- function( options(mc.cores = parallel::detectCores()) # fit model - fit_model <- rstan::sampling( - stanmodels$COA_Tag_Integrated, - data = standata, - ... - ) + if (decay == "gaussian") { + fit_model <- rstan::sampling( + stanmodels$COA_Tag_Integrated_gaussian, + data = standata, + ... + ) + } else if (decay == "logistic") { + fit_model <- rstan::sampling( + stanmodels$COA_Tag_Integrated_logistic, + data = standata, + ... + ) + } else { + stop("decay parameter must be one of \"gaussian\" or \"logistic\".") + } # Save chains after discarding warmup - fit_estimates <- as.data.frame(fit_model) # Note this returns parameters and latent states/derived values + fit_estimates <- as.data.frame(fit_model) + # Note this returns parameters and latent states/derived values # Summary statistics and convergence diagnostics - fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary - #fit_summary <- fit_sum$summary + if (decay == "gaussian") { + fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary + } else if (decay == "logistic") { + fit_summary <- rstan::summary(fit_model, pars = c("p0"))$summary + } # How much time did fitting take (in minutes)? fit_time <- sum(print(rstan::get_elapsed_time(fit_model))) / 60 - # # calculate generated quantities + + # calculate generated quantities fit_generated_quantities <- generated_quantities( model = fit_model, standata = standata, - ndraws = ndraws + n_draws = n_draws ) # transform gq into matrix tran_fit_gq <- transform_gq(fit_generated_quantities) # Extract COA estimates - coas <- array(NA, dim = c(ntime, 7, nind)) + coas <- array(NA, dim = c(n_time, 7, n_ind)) dimnames(coas)[[2]] <- c( - 'time', - 'x', - 'y', - 'x_lower', - 'x_upper', - 'y_lower', - 'y_upper' + "time", + "x", + "y", + "x_lower", + "x_upper", + "y_lower", + "y_upper" ) ew <- NULL ns <- NULL - for (i in 1:nind) { - coas[, 1, i] <- seq(1, ntime, 1) + for (i in 1:n_ind) { + coas[, 1, i] <- seq(1, n_time, 1) ew <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("sx[", i, ",", sep = '')) + dplyr::starts_with(paste("x[", i, ",", sep = "")) ) ns <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("sy[", i, ",", sep = '')) + dplyr::starts_with(paste("y[", i, ",", sep = "")) ) coas[, 2, i] <- apply(ew, 2, stats::median) coas[, 3, i] <- apply(ns, 2, stats::median) @@ -128,15 +149,15 @@ COA_TagInt <- function( coas <- as.data.frame(coas[,, 1]) # Extract time-varying detection probability estimates - d_probs <- array(NA, dim = c(nrec, ntime)) + d_probs <- array(NA, dim = c(n_rec, n_time)) p0_est <- NULL - for (i in 1:ntime) { + for (i in 1:n_time) { p0_est <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("p0[", i, ",", sep = '')) + dplyr::starts_with(paste("p0[", i, ",", sep = "")) ) - for (j in 1:nrec) { + for (j in 1:n_rec) { d_probs[j, i] <- stats::median(p0_est[, j]) } } @@ -152,13 +173,13 @@ COA_TagInt <- function( tran_fit_gq ) names(model_results) <- c( - 'model', - 'summary', - 'time', - 'coas', - 'detection_probs', - 'all_estimates', - 'generated_quantities' + "model", + "summary", + "time", + "coas", + "detection_probs", + "all_estimates", + "generated_quantities" ) return(model_results) } diff --git a/R/COA_TimeVarying.R b/R/COA_TimeVarying.R index d8375ad..d7e3293 100644 --- a/R/COA_TimeVarying.R +++ b/R/COA_TimeVarying.R @@ -1,54 +1,53 @@ -# Save this file as `R/COA_TimeVarying.R` - #' Fits a time-varying Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan #' - -#' @param nind Number of tagged individuals -#' @param nrec Number of receivers -#' @param ntime Number of time steps -#' @param ntrans Number of expected transmissions per tag per time interval -#' @param y Array of detection data, where row = individual, column = time step, and matrix = receiver -#' @param recX Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency) -#' @param recY Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency) -#' @param xlim East-west boundaries of spatial extent (receiver array + buffer) -#' @param ylim North-south boundaries of spatial extent (receiver array + buffer) -#' @param ndraws to be passed to `generated_quantities`. Changes the number of draws. Default is 10. +#' @param n_ind Number of tagged individuals +#' @param n_rec Number of receivers +#' @param n_time Number of time steps +#' @param n_trans Number of expected transmissions per tag per time interval +#' @param det Array of detection data, where row = individual, column = time step, and matrix = receiver +#' @param rec_x Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency) +#' @param rec_y Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency) +#' @param x_lim East-west boundaries of spatial extent (receiver array + buffer) +#' @param y_lim North-south boundaries of spatial extent (receiver array + buffer) +#' @param decay desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian". +#' @param n_draws to be passed to `generated_quantities`. Changes the number of draws. Default is 10. #' @param ... Additional arguments passed to `sampling` from `rstan`. #' This can include setting `chains`, `iter`, `warmup`, and `control`. Please see #' `rstan::sampling()` for more info. #' -#' @return COA_TimeVarying returns an object of class `stanfit` returned by `rstan::sampling`. See the 'rstan' package documentation for details. +#' @return COA_TimeVarying returns an object of class `stanfit` returned by `rstan::sampling`. See the "rstan" package documentation for details. #' @return This function returns a list containing the following components: 1) a summary of the detection function parameters; 2) the time required for model fitting; 3) time-varying detection probabilites for each receiver; 4) the estimated COAs for each individual in each time step and 95 percent credible interval; and 5) a dataframe containing values for each parameter and latent parameter from chain iterations. These can be used to plot posterior distributions and the credible interval around each estimated COA. #' #' @export -#' + COA_TimeVarying <- function( - nind, - nrec, - ntime, - ntrans, - y, - recX, - recY, - xlim, - ylim, - ndraws = NULL, + n_ind, + n_rec, + n_time, + n_trans, + det, + rec_x, + rec_y, + x_lim, + y_lim, + decay = "gaussian", + n_draws = NULL, ... ) { # First move everything into a list standata <- list( - nind = nind, - nrec = nrec, - ntime = ntime, - ntrans = ntrans, - y = y, - recX = recX, - recY = recY, - xlim = xlim, - ylim = ylim + n_ind = n_ind, + n_rec = n_rec, + n_time = n_time, + n_trans = n_trans, + det = det, + rec_x = rec_x, + rec_y = rec_y, + x_lim = x_lim, + y_lim = y_lim ) # validate this list prior to sending it to the model - exp_len <- expected_lengths(recX = recX, recY = recY) + exp_len <- expected_lengths(rec_x = rec_x, rec_y = rec_y) validate_standata(standata, exp_len) @@ -58,14 +57,32 @@ COA_TimeVarying <- function( options(mc.cores = parallel::detectCores()) # fit model - fit_model <- rstan::sampling(stanmodels$COA_TimeVarying, data = standata, ...) + if (decay == "gaussian") { + fit_model <- rstan::sampling( + stanmodels$COA_TimeVarying_gaussian, + data = standata, + ... + ) + } else if (decay == "logistic") { + fit_model <- rstan::sampling( + stanmodels$COA_TimeVarying_logistic, + data = standata, + ... + ) + } else { + stop("decay parameter must be one of \"gaussian\" or \"logistic\".") + } # Save chains after discarding warmup - fit_estimates <- as.data.frame(fit_model) # Note this returns parameters and latent states/derived values + fit_estimates <- as.data.frame(fit_model) + # Note this returns parameters and latent states/derived values # Summary statistics and convergence diagnostics - fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary - #fit_summary <- fit_sum$summary + if (decay == "gaussian") { + fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary + } else if (decay == "logistic") { + fit_summary <- rstan::summary(fit_model, pars = c("p0"))$summary + } # How much time did fitting take (in minutes)? fit_time <- sum(print(rstan::get_elapsed_time(fit_model))) / 60 @@ -74,33 +91,33 @@ COA_TimeVarying <- function( fit_generated_quantities <- generated_quantities( model = fit_model, standata = standata, - ndraws = ndraws + n_draws = n_draws ) # transform gq into matrix tran_fit_gq <- transform_gq(fit_generated_quantities) # Extract COA estimates - coas <- array(NA, dim = c(ntime, 7, nind)) + coas <- array(NA, dim = c(n_time, 7, n_ind)) dimnames(coas)[[2]] <- c( - 'time', - 'x', - 'y', - 'x_lower', - 'x_upper', - 'y_lower', - 'y_upper' + "time", + "x", + "y", + "x_lower", + "x_upper", + "y_lower", + "y_upper" ) ew <- NULL ns <- NULL - for (i in 1:nind) { - coas[, 1, i] <- seq(1, ntime, 1) + for (i in 1:n_ind) { + coas[, 1, i] <- seq(1, n_time, 1) ew <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("sx[", i, ",", sep = '')) + dplyr::starts_with(paste("x[", i, ",", sep = "")) ) ns <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("sy[", i, ",", sep = '')) + dplyr::starts_with(paste("y[", i, ",", sep = "")) ) coas[, 2, i] <- apply(ew, 2, stats::median) coas[, 3, i] <- apply(ns, 2, stats::median) @@ -111,15 +128,15 @@ COA_TimeVarying <- function( } coas <- as.data.frame(coas[,, 1]) # Extract time-varying detection probability estimates - d_probs <- array(NA, dim = c(nrec, ntime)) + d_probs <- array(NA, dim = c(n_rec, n_time)) p0est <- NULL - for (i in 1:ntime) { + for (i in 1:n_time) { p0est <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("p0[", i, ",", sep = '')) + dplyr::starts_with(paste("p0[", i, ",", sep = "")) ) - for (j in 1:nrec) { + for (j in 1:n_rec) { d_probs[j, i] <- stats::median(p0est[, j]) } } @@ -135,13 +152,13 @@ COA_TimeVarying <- function( tran_fit_gq ) names(model_results) <- c( - 'model', - 'summary', - 'time', - 'coas', - 'detection_probs', - 'all_estimates', - 'generated_quantities' + "model", + "summary", + "time", + "coas", + "detection_probs", + "all_estimates", + "generated_quantities" ) return(model_results) } diff --git a/R/COA_standard.R b/R/COA_standard.R index acad438..0b67c0c 100644 --- a/R/COA_standard.R +++ b/R/COA_standard.R @@ -1,18 +1,17 @@ -# Save this file as `R/COA_standard.R` - #' Fits a Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan #' -#' @param nind Number of tagged individuals -#' @param nrec Number of receivers -#' @param ntime Number of time steps -#' @param ntrans Number of expected transmissions per tag per time interval -#' @param y Array of detection data, where row = individual, column = time step, and matrix = receiver -#' @param recX Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency) -#' @param recY Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency) -#' @param xlim East-west boundaries of spatial extent (receiver array + buffer) -#' @param ylim North-south boundaries of spatial extent (receiver array + buffer). -#' @param ndraws to be passed to `generated_quantities`. Changes the number of draws. Default is 10. +#' @param n_ind Number of tagged individuals +#' @param n_rec Number of receivers +#' @param n_time Number of time steps +#' @param n_trans Number of expected transmissions per tag per time interval +#' @param det Array of detection data, where row = individual, column = time step, and matrix = receiver +#' @param rec_x Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency) +#' @param rec_y Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency) +#' @param x_lim East-west boundaries of spatial extent (receiver array + buffer) +#' @param y_lim North-south boundaries of spatial extent (receiver array + buffer). +#' @param decay desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian". +#' @param n_draws to be passed to `generated_quantities`. Changes the number of draws. Default is 10. #' @param ... Additional arguments passed to `sampling` from `rstan`. #' This can include setting `chains`, `iter`, `warmup`, and `control`. Please see #' `rstan::sampling()` for more info. @@ -24,32 +23,33 @@ #' #' @export COA_Standard <- function( - nind, - nrec, - ntime, - ntrans, - y, - recX, - recY, - xlim, - ylim, - ndraws = NULL, + n_ind, + n_rec, + n_time, + n_trans, + det, + rec_x, + rec_y, + x_lim, + y_lim, + decay = "gaussian", + n_draws = NULL, ... ) { # First move everything into a list standata <- list( - nind = nind, - nrec = nrec, - ntime = ntime, - ntrans = ntrans, - y = y, - recX = recX, - recY = recY, - xlim = xlim, - ylim = ylim + n_ind = n_ind, + n_rec = n_rec, + n_time = n_time, + n_trans = n_trans, + det = det, + rec_x = rec_x, + rec_y = rec_y, + x_lim = x_lim, + y_lim = y_lim ) # validate this list prior to sending it to the model - exp_len <- expected_lengths(recX = recX, recY = recY) + exp_len <- expected_lengths(rec_x = rec_x, rec_y = rec_y) validate_standata(standata, exp_len) @@ -59,15 +59,32 @@ COA_Standard <- function( options(mc.cores = parallel::detectCores()) # fit model - fit_model <- rstan::sampling(stanmodels$COA_Standard, data = standata, ...) + if (decay == "gaussian") { + fit_model <- rstan::sampling( + stanmodels$COA_Standard_gaussian, + data = standata, + ... + ) + } else if (decay == "logistic") { + fit_model <- rstan::sampling( + stanmodels$COA_Standard_logistic, + data = standata, + ... + ) + } else { + stop("decay parameter must be one of \"gaussian\" or \"logistic\".") + } # Save chains after discarding warmup fit_estimates <- as.data.frame(fit_model) # Note this returns parameters and latent states/derived values # Summary statistics and convergence diagnostics - fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary - #fit_summary <- fit_sum$summary + if (decay == "gaussian") { + fit_summary <- rstan::summary(fit_model, pars = c("p0", "sigma"))$summary + } else if (decay == "logistic") { + fit_summary <- rstan::summary(fit_model, pars = c("p0"))$summary + } # How much time did fitting take? fit_time <- sum(print(rstan::get_elapsed_time(fit_model))) / 60 @@ -76,33 +93,33 @@ COA_Standard <- function( fit_generated_quantities <- generated_quantities( model = fit_model, standata = standata, - ndraws = ndraws + n_draws = n_draws ) # transform gq into matrix tran_fit_gq <- transform_gq(fit_generated_quantities) # Extract COA estimates - coas <- array(NA, dim = c(ntime, 7, nind)) + coas <- array(NA, dim = c(n_time, 7, n_ind)) dimnames(coas)[[2]] <- c( - 'time', - 'x', - 'y', - 'x_lower', - 'x_upper', - 'y_lower', - 'y_upper' + "time", + "x", + "y", + "x_lower", + "x_upper", + "y_lower", + "y_upper" ) ew <- NULL ns <- NULL - for (i in 1:nind) { - coas[, 1, i] <- seq(1, ntime, 1) + for (i in 1:n_ind) { + coas[, 1, i] <- seq(1, n_time, 1) ew <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("sx[", i, ",", sep = '')) + dplyr::starts_with(paste("x[", i, ",", sep = "")) ) ns <- dplyr::select( fit_estimates, - dplyr::starts_with(paste("sy[", i, ",", sep = '')) + dplyr::starts_with(paste("y[", i, ",", sep = "")) ) coas[, 2, i] <- apply(ew, 2, stats::median) coas[, 3, i] <- apply(ns, 2, stats::median) @@ -123,12 +140,12 @@ COA_Standard <- function( tran_fit_gq ) names(model_results) <- c( - 'model', - 'summary', - 'time', - 'coas', - 'all_estimates', - 'generated_quantities' + "model", + "summary", + "time", + "coas", + "all_estimates", + "generated_quantities" ) return(model_results) } diff --git a/R/generated_quantities.R b/R/generated_quantities.R index 6afb461..e70ac5c 100644 --- a/R/generated_quantities.R +++ b/R/generated_quantities.R @@ -4,24 +4,24 @@ #' #' @param model Stan model object #' @param standata Data fed to Stan model -#' @param ndraws is the number of draws to take. Default to 10. +#' @param n_draws is the number of draws to take. Default to 10. #' #' @return generated quantities from the model #' @keywords internal #' @name generated_quantities -generated_quantities <- function(model, standata, ndraws = NULL) { +generated_quantities <- function(model, standata, n_draws = NULL) { # check stan object check_stan_object(model) - # check to see if ntest is in standata this will allow for gq's to be made + # check to see if n_test is in standata this will allow for gq's to be made # for test tag - check_test_tag <- "ntest" %in% names(standata) + check_test_tag <- "n_test" %in% names(standata) # Set default number of draws - if (is.null(ndraws)) { - ndraws <- 10 + if (is.null(n_draws)) { + n_draws <- 10 } - check_num_vec_len(ndraws, vec_length = 1) + check_num_vec_len(n_draws, vec_length = 1) # extract posteriors post <- rstan::extract(model) @@ -32,57 +32,57 @@ generated_quantities <- function(model, standata, ndraws = NULL) { list2env(standata, envir = environment()) # create vector to loop over - ndraw <- seq_len(ndraws) + n_draw <- seq_len(n_draws) # list to dump stuf into - yrep_list <- vector("list", length = ndraws) + y_rep_list <- vector("list", length = n_draws) if (check_test_tag) { - yrep_test_list <- vector("list", length = ndraws) + y_rep_test_list <- vector("list", length = n_draws) } - for (k in seq_along(ndraw)) { + for (k in seq_along(n_draw)) { # grab extracted values for ndarws - draw <- ndraw[k] + draw <- n_draw[k] a1 <- alpha1[k] if (length(dim(alpha0)) %in% 3) { - # p0 has shape [ndraws, ntime, nrec] + # p0 has shape [n_draws, n_time, n_rec] p0 <- stats::plogis(alpha0[draw, , ]) } else { p0 <- stats::plogis(alpha0[draw]) } # create blank array with the name of eveyrhting - yrep <- array( + y_rep <- array( NA, - c(nind, nrec, ntime), + c(n_ind, n_rec, n_time), dimnames = list( - tag = seq_len(nind), - rec = seq_len(nrec), - time = seq_len(ntime) + tag = seq_len(n_ind), + rec = seq_len(n_rec), + time = seq_len(n_time) ) ) if (check_test_tag) { - yrep_test <- array( + y_rep_test <- array( NA, - c(ntest, nrec, ntime), + c(n_test, n_rec, n_time), dimnames = list( - tag = seq_len(ntest), - rec = seq_len(nrec), - time = seq_len(ntime) + tag = seq_len(n_test), + rec = seq_len(n_rec), + time = seq_len(n_time) ) ) } # ----- generate quantities ------ # First for number of detections for each tagged individual - for (t in 1:ntime) { - for (i in 1:nind) { - for (j in 1:nrec) { + for (t in 1:n_time) { + for (i in 1:n_ind) { + for (j in 1:n_rec) { # create distances - d <- sqrt( - (sx[draw, i, t] - recX[j])^2 + - (sy[draw, i, t] - recY[j])^2 + dist <- sqrt( + (x[draw, i, t] - rec_x[j])^2 + + (y[draw, i, t] - rec_y[j])^2 ) # make this work for when p0 is dimensions if (is.matrix(p0)) { @@ -90,38 +90,40 @@ generated_quantities <- function(model, standata, ndraws = NULL) { } else { base <- p0 } - p <- base * exp(-a1 * d^2) + p <- base * exp(-a1 * dist^2) # make sure the pobablity is above 0 p <- min(max(p, 1e-9), 1 - 1e-9) # then run int using a the iteration of transmission by probability # to get the number of detections - yrep[i, j, t] <- stats::rbinom(1, ntrans, p) + y_rep[i, j, t] <- stats::rbinom(1, n_trans, p) } } } - yrep_list[[k]] <- yrep + y_rep_list[[k]] <- y_rep - # ----- run generated quantities for ntest ------ + # ----- run generated quantities for n_test ------ if (check_test_tag) { - for (l in 1:ntime) { - for (m in 1:nrec) { - for (s in 1:ntest) { + for (l in 1:n_time) { + for (m in 1:n_rec) { + for (s in 1:n_test) { # Euclidean distance between test tag s and receiver m - td <- sqrt((testX[s] - recX[m])^2 + (testY[s] - recY[m])^2) + test_dist <- sqrt( + (test_x[s] - rec_x[m])^2 + (test_y[s] - rec_y[m])^2 + ) # Probability - ptest <- p0[l, m] * exp(-a1 * td^2) - ptest <- min(max(ptest, 1e-9), 1 - 1e-9) + p_test <- p0[l, m] * exp(-a1 * test_dist^2) + p_test <- min(max(p_test, 1e-9), 1 - 1e-9) # Simulate detection - yrep_test[s, m, l] <- stats::rbinom(1, ntrans, ptest) + y_rep_test[s, m, l] <- stats::rbinom(1, n_trans_test, p_test) } } } - yrep_test_list[[k]] <- yrep_test + y_rep_test_list[[k]] <- y_rep_test } } if (check_test_tag) { - return(list(yrep = yrep_list, testrep = yrep_test_list)) + return(list(y_rep = y_rep_list, test_rep = y_rep_test_list)) } else { - return(list(yrep = yrep_list)) + return(list(y_rep = y_rep_list)) } } diff --git a/R/stanmodels.R b/R/stanmodels.R index d29d19c..76d7c92 100644 --- a/R/stanmodels.R +++ b/R/stanmodels.R @@ -1,12 +1,15 @@ # Generated by rstantools. Do not edit by hand. # names of stan models -stanmodels <- c("COA_Standard", "COA_Tag_Integrated", "COA_TimeVarying") +stanmodels <- c("COA_Standard_gaussian", "COA_Standard_logistic", "COA_Tag_Integrated_gaussian", "COA_Tag_Integrated_logistic", "COA_TimeVarying_gaussian", "COA_TimeVarying_logistic") # load each stan module -Rcpp::loadModule("stan_fit4COA_Standard_mod", what = TRUE) -Rcpp::loadModule("stan_fit4COA_Tag_Integrated_mod", what = TRUE) -Rcpp::loadModule("stan_fit4COA_TimeVarying_mod", what = TRUE) +Rcpp::loadModule("stan_fit4COA_Standard_gaussian_mod", what = TRUE) +Rcpp::loadModule("stan_fit4COA_Standard_logistic_mod", what = TRUE) +Rcpp::loadModule("stan_fit4COA_Tag_Integrated_gaussian_mod", what = TRUE) +Rcpp::loadModule("stan_fit4COA_Tag_Integrated_logistic_mod", what = TRUE) +Rcpp::loadModule("stan_fit4COA_TimeVarying_gaussian_mod", what = TRUE) +Rcpp::loadModule("stan_fit4COA_TimeVarying_logistic_mod", what = TRUE) # instantiate each stanmodel object stanmodels <- sapply(stanmodels, function(model_name) { diff --git a/R/utils.R b/R/utils.R index 8553f52..b53e70e 100644 --- a/R/utils.R +++ b/R/utils.R @@ -87,32 +87,33 @@ check_stan_object <- function(x, arg_name = NULL) { #' Expected lengths of variables in `standata` #' -#' @param recX is the receiver or station x coordinates (e.g, lon). -#' @param recY is the receiver or station y coordinates (e.g., lat). -#' @param ntest_len is the number of reference tags which is used as length +#' @param rec_x is the receiver or station x coordinates (e.g, lon). +#' @param rec_y is the receiver or station y coordinates (e.g., lat). +#' @param n_test_len is the number of reference tags which is used as length #' by `testX` and `testY`. #' #' #' @keywords internal #' @name expected_lengths -expected_lengths <- function(recX = NULL, recY = NULL, ntest_len = NULL) { - if (!is.null(ntest_len)) { - check_num_vec_len(ntest_len, vec_length = 1, arg_name = "ntest") +expected_lengths <- function(rec_x = NULL, rec_y = NULL, n_test_len = NULL) { + if (!is.null(n_test_len)) { + check_num_vec_len(n_test_len, vec_length = 1, arg_name = "n_test") } lengths <- list( - nind = 1, - nrec = 1, - ntime = 1, - ntrans = 1, - ntest = 1, - recX = length(recX), - recY = length(recY), - xlim = 2, - ylim = 2, - testX = ntest_len, - testY = ntest_len + n_ind = 1, + n_rec = 1, + n_time = 1, + n_trans = 1, + n_trans_test = 1, + n_test = 1, + rec_x = length(rec_x), + rec_y = length(rec_y), + x_lim = 2, + y_lim = 2, + test_x = n_test_len, + test_y = n_test_len ) return(lengths) } @@ -127,11 +128,14 @@ expected_lengths <- function(recX = NULL, recY = NULL, ntest_len = NULL) { #' @name vaidate_standata validate_standata <- function(standata, lengths) { - array_vars <- intersect(c("y", "test", "testX", "testY"), names(standata)) + array_vars <- intersect( + c("det", "det_test", "test_x", "test_y"), + names(standata) + ) for (var in array_vars) { # check station locations - if (var %in% c("testX", "testY")) { + if (var %in% c("test_x", "test_y")) { check_array_tag(standata[[var]], len = lengths[[var]], arg_name = var) } else { # Check 3d array used for counts diff --git a/R/zzz.R b/R/zzz.R index 286862d..7ae0029 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -11,16 +11,16 @@ utils::globalVariables( c( "alpha0", "alpha1", - "nind", - "nrec", - "ntest", - "ntime", - "ntrans", - "recX", - "recY", - "sx", - "sy", - "testX", - "testY" + "n_ind", + "n_rec", + "n_test", + "n_time", + "n_trans", + "rec_x", + "rec_y", + "x", + "y", + "test_x", + "test_y" ) ) diff --git a/inst/stan/COA_Standard.stan b/inst/stan/COA_Standard.stan deleted file mode 100644 index adcf85a..0000000 --- a/inst/stan/COA_Standard.stan +++ /dev/null @@ -1,51 +0,0 @@ -// Declare data -data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntrans; // number of trials/expected number of transmissions per time step - array[nind, ntime, nrec] int y; // number of detections for each individual at each receiver in each time step - vector[nrec] recX; // receiver locations in east-west direction - vector[nrec] recY; // receiver locations in north-south direction - vector[2] xlim; // area bounds east-west - vector[2] ylim; // area boundes north-south -} - -// Declare parameters -parameters { - // fixed effects - real alpha0; // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - matrix [nind, ntime] sx; // E-W center of activity coordinate - bounds reflect spatial extent - matrix [nind, ntime] sy; // N-S center of activity coordinate - bounds reflect spatial extent -} - -model { - // priors - alpha0 ~ cauchy(0, 2.5); - alpha1 ~ cauchy(0, 2.5); - - for (i in 1:nind) { - for (t in 1:ntime) { - // Calculate squared distance (d2) - vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]); - - // Run binomial on logit scale - y[i, t] ~ binomial_logit(ntrans, alpha0 - (alpha1 * d2)); - } - } -} - -generated quantities { - // Detection probability at a distance of 0 - // Inverse logit of alpha0 - constrains probability b/tw 0 and 1 - real p0 = inv_logit(alpha0); - - // Standard deviation of the distance-decay function - // Derived from coefficient specifying distance-related decay in detection prob. - // (this is 1 / 2 * sigma^2 = a1) solved to equal sigma - this then is used - // in the full model - real sigma = sqrt(1.0 / (2.0 * alpha1)); -} diff --git a/inst/stan/COA_Standard_gaussian.stan b/inst/stan/COA_Standard_gaussian.stan new file mode 100644 index 0000000..d363110 --- /dev/null +++ b/inst/stan/COA_Standard_gaussian.stan @@ -0,0 +1,51 @@ +// Declare data +data { + int n_ind; // number of individuals + int n_rec; // number of receivers + int n_time; // number of time steps + int n_trans; // number of trials/expected number of transmissions per time step + array[n_ind, n_time, n_rec] int det; // number of detections for each individual at each receiver in each time step + vector[n_rec] rec_x; // receiver locations in east-west direction + vector[n_rec] rec_y; // receiver locations in north-south direction + vector[2] x_lim; // area bounds east-west + vector[2] y_lim; // area boundes north-south +} + +// Declare parameters +parameters { + // fixed effects + real alpha0; // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space + real alpha1; // coef. for decline in detection probability with distance + + // latent variables + matrix [n_ind, n_time] x; // E-W center of activity coordinate - bounds reflect spatial extent + matrix [n_ind, n_time] y; // N-S center of activity coordinate - bounds reflect spatial extent +} + +model { + // priors + alpha0 ~ cauchy(0, 2.5); + alpha1 ~ cauchy(0, 2.5); + + for (i in 1:n_ind) { + for (t in 1:n_time) { + // Calculate squared distance (d2) + vector[n_rec] d2 = square(rec_x - x[i, t]) + square(rec_y - y[i, t]); + + // Run binomial on logit scale + det[i, t] ~ binomial_logit(n_trans, alpha0 - (alpha1 * d2)); + } + } +} + +generated quantities { + // Detection probability at a distance of 0 + // Inverse logit of alpha0 - constrains probability b/tw 0 and 1 + real p0 = inv_logit(alpha0); + + // Standard deviation of the distance-decay function + // Derived from coefficient specifying distance-related decay in detection prob. + // (this is 1 / 2 * sigma^2 = a1) solved to equal sigma - this then is used + // in the full model + real sigma = sqrt(1.0 / (2.0 * alpha1)); +} diff --git a/inst/stan/COA_Standard_logistic.stan b/inst/stan/COA_Standard_logistic.stan new file mode 100644 index 0000000..6a34e8f --- /dev/null +++ b/inst/stan/COA_Standard_logistic.stan @@ -0,0 +1,45 @@ +// Declare data +data { + int n_ind; // number of individuals + int n_rec; // number of receivers + int n_time; // number of time steps + int n_trans; // number of trials/expected number of transmissions per time step + array[n_ind, n_time, n_rec] int det; // number of detections for each individual at each receiver in each time step + vector[n_rec] rec_x; // receiver locations in east-west direction + vector[n_rec] rec_y; // receiver locations in north-south direction + vector[2] x_lim; // area bounds east-west + vector[2] y_lim; // area boundes north-south +} + +// Declare parameters +parameters { + // fixed effects + real alpha0; // detection probability intercept on the logit scale - bounds are to ensure only searching reasonable parameter space + real alpha1; // coef. for decline in detection probability with distance + + // latent variables + matrix [n_ind, n_time] x; // E-W center of activity coordinate - bounds reflect spatial extent + matrix [n_ind, n_time] y; // N-S center of activity coordinate - bounds reflect spatial extent +} + +model { + // priors + alpha0 ~ cauchy(0, 2.5); + alpha1 ~ cauchy(0, 2.5); + + for (i in 1:n_ind) { + for (t in 1:n_time) { + // Calculate euclidean distance (d) + vector[n_rec] d = sqrt(square(rec_x - x[i, t]) + square(rec_y - y[i, t])); + + // Run binomial on logit scale + det[i, t] ~ binomial_logit(n_trans, alpha0 - (alpha1 * d)); + } + } +} + +generated quantities { + // Detection probability at a distance of 0 + // Inverse logit of alpha0 - constrains probability b/tw 0 and 1 + real p0 = inv_logit(alpha0); +} diff --git a/inst/stan/COA_Tag_Integrated.stan b/inst/stan/COA_Tag_Integrated.stan deleted file mode 100644 index e0d52fe..0000000 --- a/inst/stan/COA_Tag_Integrated.stan +++ /dev/null @@ -1,77 +0,0 @@ -// Declare data -data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntest; // number of test tags - - // number of trials/expected number of transmissions per time step - int ntrans; - // number of detections for each individual at each receiver in each time step - array[nind, ntime, nrec] int y; - // number of detections from each test tag at each receiver in each time step - array[ntest, ntime, nrec] int test; - - vector[nrec] recX; // trap locations in east-west direction - vector[nrec] recY; // trap locations in north-south direction - vector[2] xlim; // area bounds east-west - vector[2] ylim; // area boundes north-south - vector[ntest] testX; // test tag locations east-west - vector[ntest] testY; // test tag locations north-south -} - -transformed data { - // Pre-calculate squared distances from receiver for fixed test tags - matrix[ntest, nrec] td2; - for (s in 1:ntest) { - td2[s, ] = to_row_vector(square(recX - testX[s]) + square(recY - testY[s])); - } -} - -// Declare parameters -parameters { - // fixed effects - // detection probability intercept - max of ~1 - matrix[ntime, nrec] alpha0; // time effect - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - // E-W center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sx; - // N-S center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sy; -} - -// Model specification -model { - // priors - to_vector(alpha0) ~ cauchy(0, 2.5); - alpha1 ~ cauchy(0, 2.5); - - // Likelihood for test tags (fixed locations) - for (s in 1:ntest) { // For each test tag - // decay over distance portion of the binomial model - vector[nrec] dist_decay = alpha1 * to_vector(td2[s, ]); - - for (t in 1:ntime) { // Run binomial on logit scale - // row(p0, t) pulls the alpha0 vector for all receivers at time t - test[s, t] ~ binomial_logit(ntrans, row(alpha0, t)' - dist_decay); - } - } - - // Likelihood for individual COAs (estimated locations) - for (i in 1:nind) { - for (t in 1:ntime) { - // Distance squared from each COA to each receiver at each time - vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]); - - // Run binomial on logit scale - y[i, t] ~ binomial_logit(ntrans, row(alpha0, t)' - (alpha1 * d2)); - } - } -} - -generated quantities { - matrix[ntime, nrec] p0 = inv_logit(alpha0); - real sigma = sqrt(1.0 / (2.0 * alpha1)); -} diff --git a/inst/stan/COA_Tag_Integrated_gaussian.stan b/inst/stan/COA_Tag_Integrated_gaussian.stan new file mode 100644 index 0000000..00c6faf --- /dev/null +++ b/inst/stan/COA_Tag_Integrated_gaussian.stan @@ -0,0 +1,82 @@ +// Declare data +data { + int n_ind; // number of individuals + int n_rec; // number of receivers + int n_time; // number of time steps + int n_test; // number of test tags + + // number of trials/expected number of transmissions per time step + int n_trans; + // number of trials/expected number of transmissions per time step for test tag + int n_trans_test; + // number of detections for each individual at each receiver in each time step + array[n_ind, n_time, n_rec] int det; + // number of detections from each test tag at each receiver in each time step + array[n_test, n_time, n_rec] int det_test; + + vector[n_rec] rec_x; // trap locations in east-west direction + vector[n_rec] rec_y; // trap locations in north-south direction + vector[2] x_lim; // area bounds east-west + vector[2] y_lim; // area boundes north-south + vector[n_test] test_x; // test tag locations east-west + vector[n_test] test_y; // test tag locations north-south +} + +transformed data { + // Pre-calculate squared distances from receiver for fixed test tags + matrix[n_test, n_rec] td2; + for (s in 1:n_test) { + td2[s, ] = to_row_vector(square(rec_x - test_x[s]) + square(rec_y - test_y[s])); + } +} + +// Declare parameters +parameters { + // fixed effects + // detection probability intercept - max of ~1 + matrix[n_time, n_rec] alpha0; // time effect + real alpha1; // coef. for decline in detection probability with distance + + // latent variables + // E-W center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] x; + // N-S center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] y; +} + +// Model specification +model { + // priors + to_vector(alpha0) ~ cauchy(0, 2.5); + alpha1 ~ cauchy(0, 2.5); + + // Likelihood for test tags (fixed locations) + for (s in 1:n_test) { // For each test tag + // decay over distance portion of the binomial model + vector[n_rec] dist_decay = alpha1 * td2[s, ]'; + + for (t in 1:n_time) { // Run binomial on logit scale + // row(p0, t) pulls the alpha0 vector for all receivers at time t + det_test[s, t] ~ binomial_logit(n_trans_test, row(alpha0, t)' - dist_decay); + } + } + + // Likelihood for individual COAs (estimated locations) + for (i in 1:n_ind) { + for (t in 1:n_time) { + // Distance squared from each COA to each receiver at each time + vector[n_rec] d2 = square(rec_x - x[i, t]) + square(rec_y - y[i, t]); + + // Calculate linear predictor + vector[n_rec] lp = row(alpha0, t)' - (alpha1 * d2); + + // Run binomial on logit scale + det[i, t] ~ binomial_logit(n_trans, lp); + } + } +} + +generated quantities { + matrix[n_time, n_rec] p0 = inv_logit(alpha0); + real sigma = sqrt(1.0 / (2.0 * alpha1)); +} diff --git a/inst/stan/COA_Tag_Integrated_logistic.stan b/inst/stan/COA_Tag_Integrated_logistic.stan new file mode 100644 index 0000000..aaec225 --- /dev/null +++ b/inst/stan/COA_Tag_Integrated_logistic.stan @@ -0,0 +1,81 @@ +// Declare data +data { + int n_ind; // number of individuals + int n_rec; // number of receivers + int n_time; // number of time steps + int n_test; // number of test tags + + // number of trials/expected number of transmissions per time step + int n_trans; + // number of trials/expected number of transmissions per time step for test tag + int n_trans_test; + // number of detections for each individual at each receiver in each time step + array[n_ind, n_time, n_rec] int det; + // number of detections from each test tag at each receiver in each time step + array[n_test, n_time, n_rec] int det_test; + + vector[n_rec] rec_x; // trap locations in east-west direction + vector[n_rec] rec_y; // trap locations in north-south direction + vector[2] x_lim; // area bounds east-west + vector[2] y_lim; // area boundes north-south + vector[n_test] test_x; // test tag locations east-west + vector[n_test] test_y; // test tag locations north-south +} + +transformed data { + // Pre-calculate squared distances from receiver for fixed test tags + matrix[n_test, n_rec] td; + for (s in 1:n_test) { + td[s, ] = to_row_vector(sqrt(square(rec_x - test_x[s]) + square(rec_y - test_y[s]))); + } +} + +// Declare parameters +parameters { + // fixed effects + // detection probability intercept - max of ~1 + matrix[n_time, n_rec] alpha0; // time effect + real alpha1; // coef. for decline in detection probability with distance + + // latent variables + // E-W center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] x; + // N-S center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] y; +} + +// Model specification +model { + // priors + to_vector(alpha0) ~ cauchy(0, 2.5); + alpha1 ~ cauchy(0, 2.5); + + // Likelihood for test tags (fixed locations) + for (s in 1:n_test) { // For each test tag + // decay over distance portion of the binomial model + vector[n_rec] dist_decay = alpha1 * td[s, ]'; + + for (t in 1:n_time) { // Run binomial on logit scale + // row(p0, t) pulls the alpha0 vector for all receivers at time t + det_test[s, t] ~ binomial_logit(n_trans_test, row(alpha0, t)' - dist_decay); + } + } + + // Likelihood for individual COAs (estimated locations) + for (i in 1:n_ind) { + for (t in 1:n_time) { + // Distance squared from each COA to each receiver at each time + vector[n_rec] d = sqrt(square(rec_x - x[i, t]) + square(rec_y - y[i, t])); + + // Calculate linear predictor + vector[n_rec] lp = row(alpha0, t)' - (alpha1 * d); + + // Run binomial on logit scale + det[i, t] ~ binomial_logit(n_trans, lp); + } + } +} + +generated quantities { + matrix[n_time, n_rec] p0 = inv_logit(alpha0); +} diff --git a/inst/stan/COA_TimeVarying.stan b/inst/stan/COA_TimeVarying.stan deleted file mode 100644 index d138cf3..0000000 --- a/inst/stan/COA_TimeVarying.stan +++ /dev/null @@ -1,52 +0,0 @@ -// Declare data -data { - int nind; // number of individuals - int nrec; // number of receivers - int ntime; // number of time steps - int ntrans; // number of trials/expected number of transmissions per time step - array[nind, ntime, nrec] int y; // number of detections for each individual at each receiver in each time step - vector[nrec] recX; // trap locations in east-west direction - vector[nrec] recY; // trap locations in north-south direction - vector[2] xlim; // area bounds east-west - vector[2] ylim; // area boundes north-south -} - -// Declare parameters -parameters { - // fixed effects - matrix[ntime, nrec] alpha0; // time effect - real alpha1; // coef. for decline in detection probability with distance - - // latent variables - // E-W center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sx; - // N-S center of activity coordinate - bounds reflect spatial extent - matrix[nind, ntime] sy; -} - -// Model specification -model { - // priors - to_vector(alpha0) ~ cauchy(0, 2.5); - alpha1 ~ cauchy(0, 2.5); - - // likelihood - for (i in 1:nind) { - for (t in 1:ntime) { - // Calculate squared distance - vector[nrec] d2 = square(recX - sx[i, t]) + square(recY - sy[i, t]); - - // Run binomial on logit scale - // row(alpha0, t) is a row_vector; ' converts to column vector - y[i, t] ~ binomial_logit(ntrans, row(alpha0, t)' - (alpha1 * d2)); - } - } -} - -generated quantities { - // Detection probability at a distance of 0 - time-varying - matrix[ntime, nrec] p0 = inv_logit(alpha0); - // Standard deviation of the distance-decay function - assume constant - real sigma = sqrt(1.0 / (2.0 * alpha1)); -} - diff --git a/inst/stan/COA_TimeVarying_gaussian.stan b/inst/stan/COA_TimeVarying_gaussian.stan new file mode 100644 index 0000000..e8ea4d0 --- /dev/null +++ b/inst/stan/COA_TimeVarying_gaussian.stan @@ -0,0 +1,52 @@ +// Declare data +data { + int n_ind; // number of individuals + int n_rec; // number of receivers + int n_time; // number of time steps + int n_trans; // number of trials/expected number of transmissions per time step + array[n_ind, n_time, n_rec] int det; // number of detections for each individual at each receiver in each time step + vector[n_rec] rec_x; // trap locations in east-west direction + vector[n_rec] rec_y; // trap locations in north-south direction + vector[2] x_lim; // area bounds east-west + vector[2] y_lim; // area boundes north-south +} + +// Declare parameters +parameters { + // fixed effects + matrix[n_time, n_rec] alpha0; // time effect + real alpha1; // coef. for decline in detection probability with distance + + // latent variables + // E-W center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] x; + // N-S center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] y; +} + +// Model specification +model { + // priors + to_vector(alpha0) ~ cauchy(0, 2.5); + alpha1 ~ cauchy(0, 2.5); + + // likelihood + for (i in 1:n_ind) { + for (t in 1:n_time) { + // Calculate squared distance + vector[n_rec] d2 = square(rec_x - x[i, t]) + square(rec_y - y[i, t]); + + // Run binomial on logit scale + // row(alpha0, t) is a row_vector; ' converts to column vector + det[i, t] ~ binomial_logit(n_trans, row(alpha0, t)' - (alpha1 * d2)); + } + } +} + +generated quantities { + // Detection probability at a distance of 0 - time-varying + matrix[n_time, n_rec] p0 = inv_logit(alpha0); + // Standard deviation of the distance-decay function - assume constant + real sigma = sqrt(1.0 / (2.0 * alpha1)); +} + diff --git a/inst/stan/COA_TimeVarying_logistic.stan b/inst/stan/COA_TimeVarying_logistic.stan new file mode 100644 index 0000000..9d56362 --- /dev/null +++ b/inst/stan/COA_TimeVarying_logistic.stan @@ -0,0 +1,50 @@ +// Declare data +data { + int n_ind; // number of individuals + int n_rec; // number of receivers + int n_time; // number of time steps + int n_trans; // number of trials/expected number of transmissions per time step + array[n_ind, n_time, n_rec] int det; // number of detections for each individual at each receiver in each time step + vector[n_rec] rec_x; // trap locations in east-west direction + vector[n_rec] rec_y; // trap locations in north-south direction + vector[2] x_lim; // area bounds east-west + vector[2] y_lim; // area boundes north-south +} + +// Declare parameters +parameters { + // fixed effects + matrix[n_time, n_rec] alpha0; // time effect + real alpha1; // coef. for decline in detection probability with distance + + // latent variables + // E-W center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] x; + // N-S center of activity coordinate - bounds reflect spatial extent + matrix[n_ind, n_time] y; +} + +// Model specification +model { + // priors + to_vector(alpha0) ~ cauchy(0, 2.5); + alpha1 ~ cauchy(0, 2.5); + + // likelihood + for (i in 1:n_ind) { + for (t in 1:n_time) { + // Calculate squared distance + vector[n_rec] d = sqrt(square(rec_x - x[i, t]) + square(rec_y - y[i, t])); + + // Run binomial on logit scale + // row(alpha0, t) is a row_vector; ' converts to column vector + det[i, t] ~ binomial_logit(n_trans, row(alpha0, t)' - (alpha1 * d)); + } + } +} + +generated quantities { + // Detection probability at a distance of 0 - time-varying + matrix[n_time, n_rec] p0 = inv_logit(alpha0); +} + diff --git a/man/COA_Standard.Rd b/man/COA_Standard.Rd index c4e1be3..913a1a1 100644 --- a/man/COA_Standard.Rd +++ b/man/COA_Standard.Rd @@ -5,39 +5,42 @@ \title{Fits a Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan} \usage{ COA_Standard( - nind, - nrec, - ntime, - ntrans, - y, - recX, - recY, - xlim, - ylim, - ndraws = NULL, + n_ind, + n_rec, + n_time, + n_trans, + det, + rec_x, + rec_y, + x_lim, + y_lim, + decay = "gaussian", + n_draws = NULL, ... ) } \arguments{ -\item{nind}{Number of tagged individuals} +\item{n_ind}{Number of tagged individuals} -\item{nrec}{Number of receivers} +\item{n_rec}{Number of receivers} -\item{ntime}{Number of time steps} +\item{n_time}{Number of time steps} -\item{ntrans}{Number of expected transmissions per tag per time interval} +\item{n_trans}{Number of expected transmissions per tag per time interval} -\item{y}{Array of detection data, where row = individual, column = time step, and matrix = receiver} +\item{det}{Array of detection data, where row = individual, column = time step, and matrix = receiver} -\item{recX}{Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)} +\item{rec_x}{Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)} -\item{recY}{Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)} +\item{rec_y}{Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)} -\item{xlim}{East-west boundaries of spatial extent (receiver array + buffer)} +\item{x_lim}{East-west boundaries of spatial extent (receiver array + buffer)} -\item{ylim}{North-south boundaries of spatial extent (receiver array + buffer).} +\item{y_lim}{North-south boundaries of spatial extent (receiver array + buffer).} -\item{ndraws}{to be passed to \code{generated_quantities}. Changes the number of draws. Default is 10.} +\item{decay}{desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian".} + +\item{n_draws}{to be passed to \code{generated_quantities}. Changes the number of draws. Default is 10.} \item{...}{Additional arguments passed to \code{sampling} from \code{rstan}. This can include setting \code{chains}, \code{iter}, \code{warmup}, and \code{control}. Please see diff --git a/man/COA_TagInt.Rd b/man/COA_TagInt.Rd index 1a331e9..81ecba3 100644 --- a/man/COA_TagInt.Rd +++ b/man/COA_TagInt.Rd @@ -5,51 +5,57 @@ \title{Fits a test-tag integrated Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan} \usage{ COA_TagInt( - nind, - nrec, - ntime, - ntest, - ntrans, - y, - test, - recX, - recY, - xlim, - ylim, - testX, - testY, - ndraws = NULL, + n_ind, + n_rec, + n_time, + n_test, + n_trans, + n_trans_test, + det, + det_test, + rec_x, + rec_y, + x_lim, + y_lim, + test_x, + test_y, + decay = "gaussian", + n_draws = NULL, ... ) } \arguments{ -\item{nind}{Number of tagged individuals} +\item{n_ind}{Number of tagged individuals} -\item{nrec}{Number of receivers} +\item{n_rec}{Number of receivers} -\item{ntime}{Number of time steps} +\item{n_time}{Number of time steps} -\item{ntest}{Number of test tags} +\item{n_test}{Number of test tags} -\item{ntrans}{Number of expected transmissions per tag per time interval} +\item{n_trans}{Number of expected transmissions per tag per time interval} -\item{y}{Array of detection data, where row = individual, column = time step, and matrix = receiver} +\item{n_trans_test}{Number of expected transmissions per sentinel tag per time interval} -\item{test}{Array of test tag detection data, where row = individual tag, column = time step, and matrix = receiver} +\item{det}{Array of detection data, where row = individual, column = receiver, and matrix = time step} -\item{recX}{Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)} +\item{det_test}{Array of test tag detection data, where row = individual tag, column = receiver, and matrix = time step} -\item{recY}{Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)} +\item{rec_x}{Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)} -\item{xlim}{East-west boundaries of spatial extent (receiver array + buffer)} +\item{rec_y}{Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)} -\item{ylim}{North-south boundaries of spatial extent (receiver array + buffer)} +\item{x_lim}{East-west boundaries of spatial extent (receiver array + buffer)} -\item{testX}{Test tag coordinates in the east-west direction (should be projected and scaled for computational efficiency)} +\item{y_lim}{North-south boundaries of spatial extent (receiver array + buffer)} -\item{testY}{Test tag coordinates in the north-south direction (should be projected and scaled for computational efficiency)} +\item{test_x}{Test tag coordinates in the east-west direction (should be projected and scaled for computational efficiency)} -\item{ndraws}{to be passed to \code{generated_quantities}. Changes the number of draws. Default is 10.} +\item{test_y}{Test tag coordinates in the north-south direction (should be projected and scaled for computational efficiency)} + +\item{decay}{desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian".} + +\item{n_draws}{to be passed to \code{generated_quantities}. Changes the number of draws. Default is 10.} \item{...}{Additional arguments passed to \code{sampling} from \code{rstan}. This can include setting \code{chains}, \code{iter}, \code{warmup}, and \code{control}. Please see diff --git a/man/COA_TimeVarying.Rd b/man/COA_TimeVarying.Rd index 10172b5..4b0b4b9 100644 --- a/man/COA_TimeVarying.Rd +++ b/man/COA_TimeVarying.Rd @@ -5,46 +5,49 @@ \title{Fits a time-varying Bayesian Spatial Point Process model to estimate individual centers of activity from acoustic telemetry data using Stan} \usage{ COA_TimeVarying( - nind, - nrec, - ntime, - ntrans, - y, - recX, - recY, - xlim, - ylim, - ndraws = NULL, + n_ind, + n_rec, + n_time, + n_trans, + det, + rec_x, + rec_y, + x_lim, + y_lim, + decay = "gaussian", + n_draws = NULL, ... ) } \arguments{ -\item{nind}{Number of tagged individuals} +\item{n_ind}{Number of tagged individuals} -\item{nrec}{Number of receivers} +\item{n_rec}{Number of receivers} -\item{ntime}{Number of time steps} +\item{n_time}{Number of time steps} -\item{ntrans}{Number of expected transmissions per tag per time interval} +\item{n_trans}{Number of expected transmissions per tag per time interval} -\item{y}{Array of detection data, where row = individual, column = time step, and matrix = receiver} +\item{det}{Array of detection data, where row = individual, column = time step, and matrix = receiver} -\item{recX}{Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)} +\item{rec_x}{Receiver coordinates in the east-west direction (should be projected and scaled for computational efficiency)} -\item{recY}{Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)} +\item{rec_y}{Receiver coordinates in the north-south direction (should be projected and scaled for computational efficiency)} -\item{xlim}{East-west boundaries of spatial extent (receiver array + buffer)} +\item{x_lim}{East-west boundaries of spatial extent (receiver array + buffer)} -\item{ylim}{North-south boundaries of spatial extent (receiver array + buffer)} +\item{y_lim}{North-south boundaries of spatial extent (receiver array + buffer)} -\item{ndraws}{to be passed to \code{generated_quantities}. Changes the number of draws. Default is 10.} +\item{decay}{desired decay function. Currently one of "gaussian" or "logistic". Default is "gaussian".} + +\item{n_draws}{to be passed to \code{generated_quantities}. Changes the number of draws. Default is 10.} \item{...}{Additional arguments passed to \code{sampling} from \code{rstan}. This can include setting \code{chains}, \code{iter}, \code{warmup}, and \code{control}. Please see \code{rstan::sampling()} for more info.} } \value{ -COA_TimeVarying returns an object of class \code{stanfit} returned by \code{rstan::sampling}. See the 'rstan' package documentation for details. +COA_TimeVarying returns an object of class \code{stanfit} returned by \code{rstan::sampling}. See the "rstan" package documentation for details. This function returns a list containing the following components: 1) a summary of the detection function parameters; 2) the time required for model fitting; 3) time-varying detection probabilites for each receiver; 4) the estimated COAs for each individual in each time step and 95 percent credible interval; and 5) a dataframe containing values for each parameter and latent parameter from chain iterations. These can be used to plot posterior distributions and the credible interval around each estimated COA. } diff --git a/man/expected_lengths.Rd b/man/expected_lengths.Rd index 3ea40f1..198e47c 100644 --- a/man/expected_lengths.Rd +++ b/man/expected_lengths.Rd @@ -4,14 +4,14 @@ \alias{expected_lengths} \title{Expected lengths of variables in \code{standata}} \usage{ -expected_lengths(recX = NULL, recY = NULL, ntest_len = NULL) +expected_lengths(rec_x = NULL, rec_y = NULL, n_test_len = NULL) } \arguments{ -\item{recX}{is the receiver or station x coordinates (e.g, lon).} +\item{rec_x}{is the receiver or station x coordinates (e.g, lon).} -\item{recY}{is the receiver or station y coordinates (e.g., lat).} +\item{rec_y}{is the receiver or station y coordinates (e.g., lat).} -\item{ntest_len}{is the number of reference tags which is used as length +\item{n_test_len}{is the number of reference tags which is used as length by \code{testX} and \code{testY}.} } \description{ diff --git a/man/generated_quantities.Rd b/man/generated_quantities.Rd index 125e208..7cd36a3 100644 --- a/man/generated_quantities.Rd +++ b/man/generated_quantities.Rd @@ -4,14 +4,14 @@ \alias{generated_quantities} \title{Generate Quantities} \usage{ -generated_quantities(model, standata, ndraws = NULL) +generated_quantities(model, standata, n_draws = NULL) } \arguments{ \item{model}{Stan model object} \item{standata}{Data fed to Stan model} -\item{ndraws}{is the number of draws to take. Default to 10.} +\item{n_draws}{is the number of draws to take. Default to 10.} } \value{ generated quantities from the model diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3eaf92d..85b871c 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,14 +12,20 @@ Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif -RcppExport SEXP _rcpp_module_boot_stan_fit4COA_Standard_mod(); -RcppExport SEXP _rcpp_module_boot_stan_fit4COA_Tag_Integrated_mod(); -RcppExport SEXP _rcpp_module_boot_stan_fit4COA_TimeVarying_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4COA_Standard_gaussian_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4COA_Standard_logistic_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4COA_Tag_Integrated_gaussian_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4COA_Tag_Integrated_logistic_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4COA_TimeVarying_gaussian_mod(); +RcppExport SEXP _rcpp_module_boot_stan_fit4COA_TimeVarying_logistic_mod(); static const R_CallMethodDef CallEntries[] = { - {"_rcpp_module_boot_stan_fit4COA_Standard_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_Standard_mod, 0}, - {"_rcpp_module_boot_stan_fit4COA_Tag_Integrated_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_Tag_Integrated_mod, 0}, - {"_rcpp_module_boot_stan_fit4COA_TimeVarying_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_TimeVarying_mod, 0}, + {"_rcpp_module_boot_stan_fit4COA_Standard_gaussian_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_Standard_gaussian_mod, 0}, + {"_rcpp_module_boot_stan_fit4COA_Standard_logistic_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_Standard_logistic_mod, 0}, + {"_rcpp_module_boot_stan_fit4COA_Tag_Integrated_gaussian_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_Tag_Integrated_gaussian_mod, 0}, + {"_rcpp_module_boot_stan_fit4COA_Tag_Integrated_logistic_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_Tag_Integrated_logistic_mod, 0}, + {"_rcpp_module_boot_stan_fit4COA_TimeVarying_gaussian_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_TimeVarying_gaussian_mod, 0}, + {"_rcpp_module_boot_stan_fit4COA_TimeVarying_logistic_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4COA_TimeVarying_logistic_mod, 0}, {NULL, NULL, 0} }; diff --git a/src/RcppExports.o b/src/RcppExports.o deleted file mode 100644 index 56322bc..0000000 Binary files a/src/RcppExports.o and /dev/null differ diff --git a/src/TelemetrySpace.so b/src/TelemetrySpace.so deleted file mode 100755 index 175c658..0000000 Binary files a/src/TelemetrySpace.so and /dev/null differ diff --git a/src/stanExports_COA_Standard.cc b/src/stanExports_COA_Standard.cc deleted file mode 100644 index e039998..0000000 --- a/src/stanExports_COA_Standard.cc +++ /dev/null @@ -1,32 +0,0 @@ -// Generated by rstantools. Do not edit by hand. - -#include -using namespace Rcpp ; -#include "stanExports_COA_Standard.h" - -RCPP_MODULE(stan_fit4COA_Standard_mod) { - - - class_ >("rstantools_model_COA_Standard") - - .constructor() - - - .method("call_sampler", &rstan::stan_fit ::call_sampler) - .method("param_names", &rstan::stan_fit ::param_names) - .method("param_names_oi", &rstan::stan_fit ::param_names_oi) - .method("param_fnames_oi", &rstan::stan_fit ::param_fnames_oi) - .method("param_dims", &rstan::stan_fit ::param_dims) - .method("param_dims_oi", &rstan::stan_fit ::param_dims_oi) - .method("update_param_oi", &rstan::stan_fit ::update_param_oi) - .method("param_oi_tidx", &rstan::stan_fit ::param_oi_tidx) - .method("grad_log_prob", &rstan::stan_fit ::grad_log_prob) - .method("log_prob", &rstan::stan_fit ::log_prob) - .method("unconstrain_pars", &rstan::stan_fit ::unconstrain_pars) - .method("constrain_pars", &rstan::stan_fit ::constrain_pars) - .method("num_pars_unconstrained", &rstan::stan_fit ::num_pars_unconstrained) - .method("unconstrained_param_names", &rstan::stan_fit ::unconstrained_param_names) - .method("constrained_param_names", &rstan::stan_fit ::constrained_param_names) - .method("standalone_gqs", &rstan::stan_fit ::standalone_gqs) - ; -} diff --git a/src/stanExports_COA_Standard.h b/src/stanExports_COA_Standard.h deleted file mode 100644 index 46e8cf2..0000000 --- a/src/stanExports_COA_Standard.h +++ /dev/null @@ -1,773 +0,0 @@ -// Generated by rstantools. Do not edit by hand. - -/* - TelemetrySpace is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - TelemetrySpace is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with TelemetrySpace. If not, see . -*/ -#ifndef MODELS_HPP -#define MODELS_HPP -#define STAN__SERVICES__COMMAND_HPP -#ifndef USE_STANC3 -#define USE_STANC3 -#endif -#include -// Code generated by stanc v2.32.2 -#include -namespace model_COA_Standard_namespace { -using stan::model::model_base_crtp; -using namespace stan::math; -stan::math::profile_map profiles__; -static constexpr std::array locations_array__ = - {" (found before start of program)", - " (in 'COA_Standard', line 16, column 2 to column 37)", - " (in 'COA_Standard', line 17, column 2 to column 25)", - " (in 'COA_Standard', line 19, column 2 to column 60)", - " (in 'COA_Standard', line 20, column 2 to column 60)", - " (in 'COA_Standard', line 39, column 2 to column 30)", - " (in 'COA_Standard', line 44, column 2 to column 42)", - " (in 'COA_Standard', line 24, column 2 to column 26)", - " (in 'COA_Standard', line 25, column 2 to column 26)", - " (in 'COA_Standard', line 29, column 13 to column 17)", - " (in 'COA_Standard', line 29, column 6 to column 74)", - " (in 'COA_Standard', line 32, column 6 to column 63)", - " (in 'COA_Standard', line 27, column 23 to line 33, column 5)", - " (in 'COA_Standard', line 27, column 4 to line 33, column 5)", - " (in 'COA_Standard', line 26, column 20 to line 34, column 3)", - " (in 'COA_Standard', line 26, column 2 to line 34, column 3)", - " (in 'COA_Standard', line 3, column 2 to column 22)", - " (in 'COA_Standard', line 4, column 2 to column 22)", - " (in 'COA_Standard', line 5, column 2 to column 23)", - " (in 'COA_Standard', line 6, column 2 to column 24)", - " (in 'COA_Standard', line 7, column 8 to column 12)", - " (in 'COA_Standard', line 7, column 14 to column 19)", - " (in 'COA_Standard', line 7, column 21 to column 25)", - " (in 'COA_Standard', line 7, column 2 to column 44)", - " (in 'COA_Standard', line 8, column 9 to column 13)", - " (in 'COA_Standard', line 8, column 2 to column 20)", - " (in 'COA_Standard', line 9, column 9 to column 13)", - " (in 'COA_Standard', line 9, column 2 to column 20)", - " (in 'COA_Standard', line 10, column 2 to column 17)", - " (in 'COA_Standard', line 11, column 2 to column 17)", - " (in 'COA_Standard', line 19, column 44 to column 48)", - " (in 'COA_Standard', line 19, column 50 to column 55)", - " (in 'COA_Standard', line 20, column 44 to column 48)", - " (in 'COA_Standard', line 20, column 50 to column 55)"}; -#include -class model_COA_Standard final : public model_base_crtp { -private: - int nind; - int nrec; - int ntime; - int ntrans; - std::vector>> y; - Eigen::Matrix recX_data__; - Eigen::Matrix recY_data__; - Eigen::Matrix xlim_data__; - Eigen::Matrix ylim_data__; - Eigen::Map> recX{nullptr, 0}; - Eigen::Map> recY{nullptr, 0}; - Eigen::Map> xlim{nullptr, 0}; - Eigen::Map> ylim{nullptr, 0}; -public: - ~model_COA_Standard() {} - model_COA_Standard(stan::io::var_context& context__, unsigned int - random_seed__ = 0, std::ostream* pstream__ = nullptr) - : model_base_crtp(0) { - int current_statement__ = 0; - using local_scalar_t__ = double; - boost::ecuyer1988 base_rng__ = - stan::services::util::create_rng(random_seed__, 0); - // suppress unused var warning - (void) base_rng__; - static constexpr const char* function__ = - "model_COA_Standard_namespace::model_COA_Standard"; - // suppress unused var warning - (void) function__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - int pos__ = std::numeric_limits::min(); - pos__ = 1; - current_statement__ = 16; - context__.validate_dims("data initialization", "nind", "int", - std::vector{}); - nind = std::numeric_limits::min(); - current_statement__ = 16; - nind = context__.vals_i("nind")[(1 - 1)]; - current_statement__ = 16; - stan::math::check_greater_or_equal(function__, "nind", nind, 0); - current_statement__ = 17; - context__.validate_dims("data initialization", "nrec", "int", - std::vector{}); - nrec = std::numeric_limits::min(); - current_statement__ = 17; - nrec = context__.vals_i("nrec")[(1 - 1)]; - current_statement__ = 17; - stan::math::check_greater_or_equal(function__, "nrec", nrec, 0); - current_statement__ = 18; - context__.validate_dims("data initialization", "ntime", "int", - std::vector{}); - ntime = std::numeric_limits::min(); - current_statement__ = 18; - ntime = context__.vals_i("ntime")[(1 - 1)]; - current_statement__ = 18; - stan::math::check_greater_or_equal(function__, "ntime", ntime, 0); - current_statement__ = 19; - context__.validate_dims("data initialization", "ntrans", "int", - std::vector{}); - ntrans = std::numeric_limits::min(); - current_statement__ = 19; - ntrans = context__.vals_i("ntrans")[(1 - 1)]; - current_statement__ = 19; - stan::math::check_greater_or_equal(function__, "ntrans", ntrans, 0); - current_statement__ = 20; - stan::math::validate_non_negative_index("y", "nind", nind); - current_statement__ = 21; - stan::math::validate_non_negative_index("y", "ntime", ntime); - current_statement__ = 22; - stan::math::validate_non_negative_index("y", "nrec", nrec); - current_statement__ = 23; - context__.validate_dims("data initialization", "y", "int", - std::vector{static_cast(nind), - static_cast(ntime), static_cast(nrec)}); - y = std::vector>>(nind, - std::vector>(ntime, - std::vector(nrec, std::numeric_limits::min()))); - { - std::vector y_flat__; - current_statement__ = 23; - y_flat__ = context__.vals_i("y"); - current_statement__ = 23; - pos__ = 1; - current_statement__ = 23; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 23; - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - current_statement__ = 23; - for (int sym3__ = 1; sym3__ <= nind; ++sym3__) { - current_statement__ = 23; - stan::model::assign(y, y_flat__[(pos__ - 1)], - "assigning variable y", stan::model::index_uni(sym3__), - stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 23; - pos__ = (pos__ + 1); - } - } - } - } - current_statement__ = 23; - stan::math::check_greater_or_equal(function__, "y", y, 0); - current_statement__ = 24; - stan::math::validate_non_negative_index("recX", "nrec", nrec); - current_statement__ = 25; - context__.validate_dims("data initialization", "recX", "double", - std::vector{static_cast(nrec)}); - recX_data__ = Eigen::Matrix::Constant(nrec, - std::numeric_limits::quiet_NaN()); - new (&recX) Eigen::Map>(recX_data__.data(), - nrec); - { - std::vector recX_flat__; - current_statement__ = 25; - recX_flat__ = context__.vals_r("recX"); - current_statement__ = 25; - pos__ = 1; - current_statement__ = 25; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 25; - stan::model::assign(recX, recX_flat__[(pos__ - 1)], - "assigning variable recX", stan::model::index_uni(sym1__)); - current_statement__ = 25; - pos__ = (pos__ + 1); - } - } - current_statement__ = 26; - stan::math::validate_non_negative_index("recY", "nrec", nrec); - current_statement__ = 27; - context__.validate_dims("data initialization", "recY", "double", - std::vector{static_cast(nrec)}); - recY_data__ = Eigen::Matrix::Constant(nrec, - std::numeric_limits::quiet_NaN()); - new (&recY) Eigen::Map>(recY_data__.data(), - nrec); - { - std::vector recY_flat__; - current_statement__ = 27; - recY_flat__ = context__.vals_r("recY"); - current_statement__ = 27; - pos__ = 1; - current_statement__ = 27; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 27; - stan::model::assign(recY, recY_flat__[(pos__ - 1)], - "assigning variable recY", stan::model::index_uni(sym1__)); - current_statement__ = 27; - pos__ = (pos__ + 1); - } - } - current_statement__ = 28; - context__.validate_dims("data initialization", "xlim", "double", - std::vector{static_cast(2)}); - xlim_data__ = Eigen::Matrix::Constant(2, - std::numeric_limits::quiet_NaN()); - new (&xlim) Eigen::Map>(xlim_data__.data(), - 2); - { - std::vector xlim_flat__; - current_statement__ = 28; - xlim_flat__ = context__.vals_r("xlim"); - current_statement__ = 28; - pos__ = 1; - current_statement__ = 28; - for (int sym1__ = 1; sym1__ <= 2; ++sym1__) { - current_statement__ = 28; - stan::model::assign(xlim, xlim_flat__[(pos__ - 1)], - "assigning variable xlim", stan::model::index_uni(sym1__)); - current_statement__ = 28; - pos__ = (pos__ + 1); - } - } - current_statement__ = 29; - context__.validate_dims("data initialization", "ylim", "double", - std::vector{static_cast(2)}); - ylim_data__ = Eigen::Matrix::Constant(2, - std::numeric_limits::quiet_NaN()); - new (&ylim) Eigen::Map>(ylim_data__.data(), - 2); - { - std::vector ylim_flat__; - current_statement__ = 29; - ylim_flat__ = context__.vals_r("ylim"); - current_statement__ = 29; - pos__ = 1; - current_statement__ = 29; - for (int sym1__ = 1; sym1__ <= 2; ++sym1__) { - current_statement__ = 29; - stan::model::assign(ylim, ylim_flat__[(pos__ - 1)], - "assigning variable ylim", stan::model::index_uni(sym1__)); - current_statement__ = 29; - pos__ = (pos__ + 1); - } - } - current_statement__ = 30; - stan::math::validate_non_negative_index("sx", "nind", nind); - current_statement__ = 31; - stan::math::validate_non_negative_index("sx", "ntime", ntime); - current_statement__ = 32; - stan::math::validate_non_negative_index("sy", "nind", nind); - current_statement__ = 33; - stan::math::validate_non_negative_index("sy", "ntime", ntime); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - num_params_r__ = 1 + 1 + (nind * ntime) + (nind * ntime); - } - inline std::string model_name() const final { - return "model_COA_Standard"; - } - inline std::vector model_compile_info() const noexcept { - return std::vector{"stanc_version = stanc3 v2.32.2", - "stancflags = --allow-undefined"}; - } - template * = nullptr, - stan::require_vector_like_vt* = nullptr> - inline stan::scalar_type_t - log_prob_impl(VecR& params_r__, VecI& params_i__, std::ostream* - pstream__ = nullptr) const { - using T__ = stan::scalar_type_t; - using local_scalar_t__ = T__; - T__ lp__(0.0); - stan::math::accumulator lp_accum__; - stan::io::deserializer in__(params_r__, params_i__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - static constexpr const char* function__ = - "model_COA_Standard_namespace::log_prob"; - // suppress unused var warning - (void) function__; - try { - local_scalar_t__ alpha0 = DUMMY_VAR__; - current_statement__ = 1; - alpha0 = in__.template read_constrain_lub(-7, 7, lp__); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = in__.template read_constrain_lb(0, lp__); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 3; - sx = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), - lp__, nind, ntime); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 4; - sy = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), - lp__, nind, ntime); - { - current_statement__ = 7; - lp_accum__.add(stan::math::cauchy_lpdf(alpha0, 0, 2.5)); - current_statement__ = 8; - lp_accum__.add(stan::math::cauchy_lpdf(alpha1, 0, 2.5)); - current_statement__ = 15; - for (int i = 1; i <= nind; ++i) { - current_statement__ = 13; - for (int t = 1; t <= ntime; ++t) { - current_statement__ = 9; - stan::math::validate_non_negative_index("d2", "nrec", nrec); - Eigen::Matrix d2 = - Eigen::Matrix::Constant(nrec, - DUMMY_VAR__); - current_statement__ = 10; - stan::model::assign(d2, - stan::math::add( - stan::math::square( - stan::math::subtract(recX, - stan::model::rvalue(sx, "sx", stan::model::index_uni(i), - stan::model::index_uni(t)))), - stan::math::square( - stan::math::subtract(recY, - stan::model::rvalue(sy, "sy", stan::model::index_uni(i), - stan::model::index_uni(t))))), "assigning variable d2"); - current_statement__ = 11; - lp_accum__.add(stan::math::binomial_logit_lpmf( - stan::model::rvalue(y, "y", - stan::model::index_uni(i), - stan::model::index_uni(t)), ntrans, - stan::math::subtract(alpha0, - stan::math::multiply(alpha1, d2)))); - } - } - } - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - lp_accum__.add(lp__); - return lp_accum__.sum(); - } - template * = nullptr, stan::require_vector_like_vt* = nullptr, stan::require_vector_vt* = nullptr> - inline void - write_array_impl(RNG& base_rng__, VecR& params_r__, VecI& params_i__, - VecVar& vars__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true, std::ostream* - pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::deserializer in__(params_r__, params_i__); - stan::io::serializer out__(vars__); - static constexpr bool propto__ = true; - // suppress unused var warning - (void) propto__; - double lp__ = 0.0; - // suppress unused var warning - (void) lp__; - int current_statement__ = 0; - stan::math::accumulator lp_accum__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - constexpr bool jacobian__ = false; - static constexpr const char* function__ = - "model_COA_Standard_namespace::write_array"; - // suppress unused var warning - (void) function__; - try { - double alpha0 = std::numeric_limits::quiet_NaN(); - current_statement__ = 1; - alpha0 = in__.template read_constrain_lub(-7, 7, lp__); - double alpha1 = std::numeric_limits::quiet_NaN(); - current_statement__ = 2; - alpha1 = in__.template read_constrain_lb(0, lp__); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - std::numeric_limits::quiet_NaN()); - current_statement__ = 3; - sx = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), - lp__, nind, ntime); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - std::numeric_limits::quiet_NaN()); - current_statement__ = 4; - sy = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), - lp__, nind, ntime); - out__.write(alpha0); - out__.write(alpha1); - out__.write(sx); - out__.write(sy); - if (stan::math::logical_negation( - (stan::math::primitive_value(emit_transformed_parameters__) || - stan::math::primitive_value(emit_generated_quantities__)))) { - return ; - } - if (stan::math::logical_negation(emit_generated_quantities__)) { - return ; - } - double p0 = std::numeric_limits::quiet_NaN(); - current_statement__ = 5; - p0 = stan::math::inv_logit(alpha0); - double sigma = std::numeric_limits::quiet_NaN(); - current_statement__ = 6; - sigma = stan::math::sqrt((1.0 / (2.0 * alpha1))); - out__.write(p0); - out__.write(sigma); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - template * = nullptr, - stan::require_vector_like_vt* = nullptr> - inline void - unconstrain_array_impl(const VecVar& params_r__, const VecI& params_i__, - VecVar& vars__, std::ostream* pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::deserializer in__(params_r__, params_i__); - stan::io::serializer out__(vars__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - int pos__ = std::numeric_limits::min(); - pos__ = 1; - local_scalar_t__ alpha0 = DUMMY_VAR__; - current_statement__ = 1; - alpha0 = in__.read(); - out__.write_free_lub(-7, 7, alpha0); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = in__.read(); - out__.write_free_lb(0, alpha1); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 3; - stan::model::assign(sx, - in__.read>(nind, ntime), - "assigning variable sx"); - out__.write_free_lub(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), sx); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 4; - stan::model::assign(sy, - in__.read>(nind, ntime), - "assigning variable sy"); - out__.write_free_lub(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), sy); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - template * = nullptr> - inline void - transform_inits_impl(const stan::io::var_context& context__, VecVar& - vars__, std::ostream* pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::serializer out__(vars__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - current_statement__ = 1; - context__.validate_dims("parameter initialization", "alpha0", "double", - std::vector{}); - current_statement__ = 2; - context__.validate_dims("parameter initialization", "alpha1", "double", - std::vector{}); - current_statement__ = 3; - context__.validate_dims("parameter initialization", "sx", "double", - std::vector{static_cast(nind), - static_cast(ntime)}); - current_statement__ = 4; - context__.validate_dims("parameter initialization", "sy", "double", - std::vector{static_cast(nind), - static_cast(ntime)}); - int pos__ = std::numeric_limits::min(); - pos__ = 1; - local_scalar_t__ alpha0 = DUMMY_VAR__; - current_statement__ = 1; - alpha0 = context__.vals_r("alpha0")[(1 - 1)]; - out__.write_free_lub(-7, 7, alpha0); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = context__.vals_r("alpha1")[(1 - 1)]; - out__.write_free_lb(0, alpha1); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - { - std::vector sx_flat__; - current_statement__ = 3; - sx_flat__ = context__.vals_r("sx"); - current_statement__ = 3; - pos__ = 1; - current_statement__ = 3; - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - current_statement__ = 3; - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - current_statement__ = 3; - stan::model::assign(sx, sx_flat__[(pos__ - 1)], - "assigning variable sx", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 3; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), sx); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - { - std::vector sy_flat__; - current_statement__ = 4; - sy_flat__ = context__.vals_r("sy"); - current_statement__ = 4; - pos__ = 1; - current_statement__ = 4; - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - current_statement__ = 4; - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - current_statement__ = 4; - stan::model::assign(sy, sy_flat__[(pos__ - 1)], - "assigning variable sy", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 4; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), sy); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - inline void - get_param_names(std::vector& names__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true) const { - names__ = std::vector{"alpha0", "alpha1", "sx", "sy"}; - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - std::vector temp{"p0", "sigma"}; - names__.reserve(names__.size() + temp.size()); - names__.insert(names__.end(), temp.begin(), temp.end()); - } - } - inline void - get_dims(std::vector>& dimss__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true) const { - dimss__ = std::vector>{std::vector{}, - std::vector{}, - std::vector{static_cast(nind), - static_cast(ntime)}, - std::vector{static_cast(nind), - static_cast(ntime)}}; - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - std::vector> - temp{std::vector{}, std::vector{}}; - dimss__.reserve(dimss__.size() + temp.size()); - dimss__.insert(dimss__.end(), temp.begin(), temp.end()); - } - } - inline void - constrained_param_names(std::vector& param_names__, bool - emit_transformed_parameters__ = true, bool - emit_generated_quantities__ = true) const final { - param_names__.emplace_back(std::string() + "alpha0"); - param_names__.emplace_back(std::string() + "alpha1"); - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sx" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sy" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - param_names__.emplace_back(std::string() + "p0"); - param_names__.emplace_back(std::string() + "sigma"); - } - } - inline void - unconstrained_param_names(std::vector& param_names__, bool - emit_transformed_parameters__ = true, bool - emit_generated_quantities__ = true) const final { - param_names__.emplace_back(std::string() + "alpha0"); - param_names__.emplace_back(std::string() + "alpha1"); - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sx" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sy" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - param_names__.emplace_back(std::string() + "p0"); - param_names__.emplace_back(std::string() + "sigma"); - } - } - inline std::string get_constrained_sizedtypes() const { - return std::string("[{\"name\":\"alpha0\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"alpha1\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"sx\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"sy\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"p0\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"},{\"name\":\"sigma\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]"); - } - inline std::string get_unconstrained_sizedtypes() const { - return std::string("[{\"name\":\"alpha0\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"alpha1\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"sx\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"sy\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"p0\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"},{\"name\":\"sigma\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]"); - } - // Begin method overload boilerplate - template inline void - write_array(RNG& base_rng, Eigen::Matrix& params_r, - Eigen::Matrix& vars, const bool - emit_transformed_parameters = true, const bool - emit_generated_quantities = true, std::ostream* - pstream = nullptr) const { - const size_t num_params__ = (((1 + 1) + (nind * ntime)) + (nind * ntime)); - const size_t num_transformed = emit_transformed_parameters * (0); - const size_t num_gen_quantities = emit_generated_quantities * ((1 + 1)); - const size_t num_to_write = num_params__ + num_transformed + - num_gen_quantities; - std::vector params_i; - vars = Eigen::Matrix::Constant(num_to_write, - std::numeric_limits::quiet_NaN()); - write_array_impl(base_rng, params_r, params_i, vars, - emit_transformed_parameters, emit_generated_quantities, pstream); - } - template inline void - write_array(RNG& base_rng, std::vector& params_r, std::vector& - params_i, std::vector& vars, bool - emit_transformed_parameters = true, bool - emit_generated_quantities = true, std::ostream* - pstream = nullptr) const { - const size_t num_params__ = (((1 + 1) + (nind * ntime)) + (nind * ntime)); - const size_t num_transformed = emit_transformed_parameters * (0); - const size_t num_gen_quantities = emit_generated_quantities * ((1 + 1)); - const size_t num_to_write = num_params__ + num_transformed + - num_gen_quantities; - vars = std::vector(num_to_write, - std::numeric_limits::quiet_NaN()); - write_array_impl(base_rng, params_r, params_i, vars, - emit_transformed_parameters, emit_generated_quantities, pstream); - } - template inline T_ - log_prob(Eigen::Matrix& params_r, std::ostream* pstream = nullptr) const { - Eigen::Matrix params_i; - return log_prob_impl(params_r, params_i, pstream); - } - template inline T_ - log_prob(std::vector& params_r, std::vector& params_i, - std::ostream* pstream = nullptr) const { - return log_prob_impl(params_r, params_i, pstream); - } - inline void - transform_inits(const stan::io::var_context& context, - Eigen::Matrix& params_r, std::ostream* - pstream = nullptr) const final { - std::vector params_r_vec(params_r.size()); - std::vector params_i; - transform_inits(context, params_i, params_r_vec, pstream); - params_r = Eigen::Map>(params_r_vec.data(), - params_r_vec.size()); - } - inline void - transform_inits(const stan::io::var_context& context, std::vector& - params_i, std::vector& vars, std::ostream* - pstream__ = nullptr) const { - vars.resize(num_params_r__); - transform_inits_impl(context, vars, pstream__); - } - inline void - unconstrain_array(const std::vector& params_constrained, - std::vector& params_unconstrained, std::ostream* - pstream = nullptr) const { - const std::vector params_i; - params_unconstrained = std::vector(num_params_r__, - std::numeric_limits::quiet_NaN()); - unconstrain_array_impl(params_constrained, params_i, - params_unconstrained, pstream); - } - inline void - unconstrain_array(const Eigen::Matrix& params_constrained, - Eigen::Matrix& params_unconstrained, - std::ostream* pstream = nullptr) const { - const std::vector params_i; - params_unconstrained = Eigen::Matrix::Constant(num_params_r__, - std::numeric_limits::quiet_NaN()); - unconstrain_array_impl(params_constrained, params_i, - params_unconstrained, pstream); - } -}; -} -using stan_model = model_COA_Standard_namespace::model_COA_Standard; -#ifndef USING_R -// Boilerplate -stan::model::model_base& -new_model(stan::io::var_context& data_context, unsigned int seed, - std::ostream* msg_stream) { - stan_model* m = new stan_model(data_context, seed, msg_stream); - return *m; -} -stan::math::profile_map& get_stan_profile_data() { - return model_COA_Standard_namespace::profiles__; -} -#endif -#endif diff --git a/src/stanExports_COA_Standard.o b/src/stanExports_COA_Standard.o deleted file mode 100644 index 127e0a7..0000000 Binary files a/src/stanExports_COA_Standard.o and /dev/null differ diff --git a/src/stanExports_COA_Tag_Integrated.cc b/src/stanExports_COA_Tag_Integrated.cc deleted file mode 100644 index f4ab7ab..0000000 --- a/src/stanExports_COA_Tag_Integrated.cc +++ /dev/null @@ -1,32 +0,0 @@ -// Generated by rstantools. Do not edit by hand. - -#include -using namespace Rcpp ; -#include "stanExports_COA_Tag_Integrated.h" - -RCPP_MODULE(stan_fit4COA_Tag_Integrated_mod) { - - - class_ >("rstantools_model_COA_Tag_Integrated") - - .constructor() - - - .method("call_sampler", &rstan::stan_fit ::call_sampler) - .method("param_names", &rstan::stan_fit ::param_names) - .method("param_names_oi", &rstan::stan_fit ::param_names_oi) - .method("param_fnames_oi", &rstan::stan_fit ::param_fnames_oi) - .method("param_dims", &rstan::stan_fit ::param_dims) - .method("param_dims_oi", &rstan::stan_fit ::param_dims_oi) - .method("update_param_oi", &rstan::stan_fit ::update_param_oi) - .method("param_oi_tidx", &rstan::stan_fit ::param_oi_tidx) - .method("grad_log_prob", &rstan::stan_fit ::grad_log_prob) - .method("log_prob", &rstan::stan_fit ::log_prob) - .method("unconstrain_pars", &rstan::stan_fit ::unconstrain_pars) - .method("constrain_pars", &rstan::stan_fit ::constrain_pars) - .method("num_pars_unconstrained", &rstan::stan_fit ::num_pars_unconstrained) - .method("unconstrained_param_names", &rstan::stan_fit ::unconstrained_param_names) - .method("constrained_param_names", &rstan::stan_fit ::constrained_param_names) - .method("standalone_gqs", &rstan::stan_fit ::standalone_gqs) - ; -} diff --git a/src/stanExports_COA_Tag_Integrated.h b/src/stanExports_COA_Tag_Integrated.h deleted file mode 100644 index 9727569..0000000 --- a/src/stanExports_COA_Tag_Integrated.h +++ /dev/null @@ -1,1023 +0,0 @@ -// Generated by rstantools. Do not edit by hand. - -/* - TelemetrySpace is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - TelemetrySpace is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with TelemetrySpace. If not, see . -*/ -#ifndef MODELS_HPP -#define MODELS_HPP -#define STAN__SERVICES__COMMAND_HPP -#ifndef USE_STANC3 -#define USE_STANC3 -#endif -#include -// Code generated by stanc v2.32.2 -#include -namespace model_COA_Tag_Integrated_namespace { -using stan::model::model_base_crtp; -using namespace stan::math; -stan::math::profile_map profiles__; -static constexpr std::array locations_array__ = - {" (found before start of program)", - " (in 'COA_Tag_Integrated', line 31, column 2 to column 52)", - " (in 'COA_Tag_Integrated', line 32, column 2 to column 25)", - " (in 'COA_Tag_Integrated', line 35, column 2 to column 59)", - " (in 'COA_Tag_Integrated', line 37, column 2 to column 59)", - " (in 'COA_Tag_Integrated', line 65, column 3 to column 46)", - " (in 'COA_Tag_Integrated', line 66, column 3 to column 43)", - " (in 'COA_Tag_Integrated', line 42, column 2 to column 37)", - " (in 'COA_Tag_Integrated', line 43, column 2 to column 26)", - " (in 'COA_Tag_Integrated', line 47, column 11 to column 15)", - " (in 'COA_Tag_Integrated', line 47, column 4 to column 59)", - " (in 'COA_Tag_Integrated', line 50, column 6 to column 72)", - " (in 'COA_Tag_Integrated', line 48, column 23 to line 51, column 5)", - " (in 'COA_Tag_Integrated', line 48, column 4 to line 51, column 5)", - " (in 'COA_Tag_Integrated', line 45, column 21 to line 52, column 3)", - " (in 'COA_Tag_Integrated', line 45, column 2 to line 52, column 3)", - " (in 'COA_Tag_Integrated', line 57, column 13 to column 17)", - " (in 'COA_Tag_Integrated', line 57, column 6 to column 74)", - " (in 'COA_Tag_Integrated', line 60, column 6 to column 72)", - " (in 'COA_Tag_Integrated', line 55, column 23 to line 61, column 5)", - " (in 'COA_Tag_Integrated', line 55, column 4 to line 61, column 5)", - " (in 'COA_Tag_Integrated', line 54, column 20 to line 62, column 3)", - " (in 'COA_Tag_Integrated', line 54, column 2 to line 62, column 3)", - " (in 'COA_Tag_Integrated', line 3, column 2 to column 22)", - " (in 'COA_Tag_Integrated', line 4, column 2 to column 22)", - " (in 'COA_Tag_Integrated', line 5, column 2 to column 23)", - " (in 'COA_Tag_Integrated', line 6, column 2 to column 23)", - " (in 'COA_Tag_Integrated', line 8, column 2 to column 24)", - " (in 'COA_Tag_Integrated', line 10, column 8 to column 12)", - " (in 'COA_Tag_Integrated', line 10, column 14 to column 19)", - " (in 'COA_Tag_Integrated', line 10, column 21 to column 25)", - " (in 'COA_Tag_Integrated', line 10, column 2 to column 44)", - " (in 'COA_Tag_Integrated', line 12, column 8 to column 13)", - " (in 'COA_Tag_Integrated', line 12, column 15 to column 20)", - " (in 'COA_Tag_Integrated', line 12, column 22 to column 26)", - " (in 'COA_Tag_Integrated', line 12, column 2 to column 48)", - " (in 'COA_Tag_Integrated', line 13, column 9 to column 13)", - " (in 'COA_Tag_Integrated', line 13, column 2 to column 20)", - " (in 'COA_Tag_Integrated', line 14, column 9 to column 13)", - " (in 'COA_Tag_Integrated', line 14, column 2 to column 20)", - " (in 'COA_Tag_Integrated', line 15, column 2 to column 17)", - " (in 'COA_Tag_Integrated', line 16, column 2 to column 17)", - " (in 'COA_Tag_Integrated', line 17, column 9 to column 14)", - " (in 'COA_Tag_Integrated', line 17, column 2 to column 22)", - " (in 'COA_Tag_Integrated', line 18, column 9 to column 14)", - " (in 'COA_Tag_Integrated', line 18, column 2 to column 22)", - " (in 'COA_Tag_Integrated', line 22, column 9 to column 14)", - " (in 'COA_Tag_Integrated', line 22, column 16 to column 20)", - " (in 'COA_Tag_Integrated', line 22, column 2 to column 26)", - " (in 'COA_Tag_Integrated', line 24, column 4 to column 80)", - " (in 'COA_Tag_Integrated', line 23, column 21 to line 25, column 3)", - " (in 'COA_Tag_Integrated', line 23, column 2 to line 25, column 3)", - " (in 'COA_Tag_Integrated', line 31, column 32 to column 37)", - " (in 'COA_Tag_Integrated', line 31, column 39 to column 43)", - " (in 'COA_Tag_Integrated', line 35, column 43 to column 47)", - " (in 'COA_Tag_Integrated', line 35, column 49 to column 54)", - " (in 'COA_Tag_Integrated', line 37, column 43 to column 47)", - " (in 'COA_Tag_Integrated', line 37, column 49 to column 54)", - " (in 'COA_Tag_Integrated', line 65, column 10 to column 15)", - " (in 'COA_Tag_Integrated', line 65, column 17 to column 21)"}; -#include -class model_COA_Tag_Integrated final : public model_base_crtp { -private: - int nind; - int nrec; - int ntime; - int ntest; - int ntrans; - std::vector>> y; - std::vector>> test; - Eigen::Matrix recX_data__; - Eigen::Matrix recY_data__; - Eigen::Matrix xlim_data__; - Eigen::Matrix ylim_data__; - Eigen::Matrix testX_data__; - Eigen::Matrix testY_data__; - Eigen::Matrix td2_data__; - Eigen::Map> recX{nullptr, 0}; - Eigen::Map> recY{nullptr, 0}; - Eigen::Map> xlim{nullptr, 0}; - Eigen::Map> ylim{nullptr, 0}; - Eigen::Map> testX{nullptr, 0}; - Eigen::Map> testY{nullptr, 0}; - Eigen::Map> td2{nullptr, 0, 0}; -public: - ~model_COA_Tag_Integrated() {} - model_COA_Tag_Integrated(stan::io::var_context& context__, unsigned int - random_seed__ = 0, std::ostream* - pstream__ = nullptr) : model_base_crtp(0) { - int current_statement__ = 0; - using local_scalar_t__ = double; - boost::ecuyer1988 base_rng__ = - stan::services::util::create_rng(random_seed__, 0); - // suppress unused var warning - (void) base_rng__; - static constexpr const char* function__ = - "model_COA_Tag_Integrated_namespace::model_COA_Tag_Integrated"; - // suppress unused var warning - (void) function__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - int pos__ = std::numeric_limits::min(); - pos__ = 1; - current_statement__ = 23; - context__.validate_dims("data initialization", "nind", "int", - std::vector{}); - nind = std::numeric_limits::min(); - current_statement__ = 23; - nind = context__.vals_i("nind")[(1 - 1)]; - current_statement__ = 23; - stan::math::check_greater_or_equal(function__, "nind", nind, 0); - current_statement__ = 24; - context__.validate_dims("data initialization", "nrec", "int", - std::vector{}); - nrec = std::numeric_limits::min(); - current_statement__ = 24; - nrec = context__.vals_i("nrec")[(1 - 1)]; - current_statement__ = 24; - stan::math::check_greater_or_equal(function__, "nrec", nrec, 0); - current_statement__ = 25; - context__.validate_dims("data initialization", "ntime", "int", - std::vector{}); - ntime = std::numeric_limits::min(); - current_statement__ = 25; - ntime = context__.vals_i("ntime")[(1 - 1)]; - current_statement__ = 25; - stan::math::check_greater_or_equal(function__, "ntime", ntime, 0); - current_statement__ = 26; - context__.validate_dims("data initialization", "ntest", "int", - std::vector{}); - ntest = std::numeric_limits::min(); - current_statement__ = 26; - ntest = context__.vals_i("ntest")[(1 - 1)]; - current_statement__ = 26; - stan::math::check_greater_or_equal(function__, "ntest", ntest, 0); - current_statement__ = 27; - context__.validate_dims("data initialization", "ntrans", "int", - std::vector{}); - ntrans = std::numeric_limits::min(); - current_statement__ = 27; - ntrans = context__.vals_i("ntrans")[(1 - 1)]; - current_statement__ = 27; - stan::math::check_greater_or_equal(function__, "ntrans", ntrans, 0); - current_statement__ = 28; - stan::math::validate_non_negative_index("y", "nind", nind); - current_statement__ = 29; - stan::math::validate_non_negative_index("y", "ntime", ntime); - current_statement__ = 30; - stan::math::validate_non_negative_index("y", "nrec", nrec); - current_statement__ = 31; - context__.validate_dims("data initialization", "y", "int", - std::vector{static_cast(nind), - static_cast(ntime), static_cast(nrec)}); - y = std::vector>>(nind, - std::vector>(ntime, - std::vector(nrec, std::numeric_limits::min()))); - { - std::vector y_flat__; - current_statement__ = 31; - y_flat__ = context__.vals_i("y"); - current_statement__ = 31; - pos__ = 1; - current_statement__ = 31; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 31; - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - current_statement__ = 31; - for (int sym3__ = 1; sym3__ <= nind; ++sym3__) { - current_statement__ = 31; - stan::model::assign(y, y_flat__[(pos__ - 1)], - "assigning variable y", stan::model::index_uni(sym3__), - stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 31; - pos__ = (pos__ + 1); - } - } - } - } - current_statement__ = 31; - stan::math::check_greater_or_equal(function__, "y", y, 0); - current_statement__ = 32; - stan::math::validate_non_negative_index("test", "ntest", ntest); - current_statement__ = 33; - stan::math::validate_non_negative_index("test", "ntime", ntime); - current_statement__ = 34; - stan::math::validate_non_negative_index("test", "nrec", nrec); - current_statement__ = 35; - context__.validate_dims("data initialization", "test", "int", - std::vector{static_cast(ntest), - static_cast(ntime), static_cast(nrec)}); - test = std::vector>>(ntest, - std::vector>(ntime, - std::vector(nrec, std::numeric_limits::min()))); - { - std::vector test_flat__; - current_statement__ = 35; - test_flat__ = context__.vals_i("test"); - current_statement__ = 35; - pos__ = 1; - current_statement__ = 35; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 35; - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - current_statement__ = 35; - for (int sym3__ = 1; sym3__ <= ntest; ++sym3__) { - current_statement__ = 35; - stan::model::assign(test, test_flat__[(pos__ - 1)], - "assigning variable test", stan::model::index_uni(sym3__), - stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 35; - pos__ = (pos__ + 1); - } - } - } - } - current_statement__ = 35; - stan::math::check_greater_or_equal(function__, "test", test, 0); - current_statement__ = 36; - stan::math::validate_non_negative_index("recX", "nrec", nrec); - current_statement__ = 37; - context__.validate_dims("data initialization", "recX", "double", - std::vector{static_cast(nrec)}); - recX_data__ = Eigen::Matrix::Constant(nrec, - std::numeric_limits::quiet_NaN()); - new (&recX) Eigen::Map>(recX_data__.data(), - nrec); - { - std::vector recX_flat__; - current_statement__ = 37; - recX_flat__ = context__.vals_r("recX"); - current_statement__ = 37; - pos__ = 1; - current_statement__ = 37; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 37; - stan::model::assign(recX, recX_flat__[(pos__ - 1)], - "assigning variable recX", stan::model::index_uni(sym1__)); - current_statement__ = 37; - pos__ = (pos__ + 1); - } - } - current_statement__ = 38; - stan::math::validate_non_negative_index("recY", "nrec", nrec); - current_statement__ = 39; - context__.validate_dims("data initialization", "recY", "double", - std::vector{static_cast(nrec)}); - recY_data__ = Eigen::Matrix::Constant(nrec, - std::numeric_limits::quiet_NaN()); - new (&recY) Eigen::Map>(recY_data__.data(), - nrec); - { - std::vector recY_flat__; - current_statement__ = 39; - recY_flat__ = context__.vals_r("recY"); - current_statement__ = 39; - pos__ = 1; - current_statement__ = 39; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 39; - stan::model::assign(recY, recY_flat__[(pos__ - 1)], - "assigning variable recY", stan::model::index_uni(sym1__)); - current_statement__ = 39; - pos__ = (pos__ + 1); - } - } - current_statement__ = 40; - context__.validate_dims("data initialization", "xlim", "double", - std::vector{static_cast(2)}); - xlim_data__ = Eigen::Matrix::Constant(2, - std::numeric_limits::quiet_NaN()); - new (&xlim) Eigen::Map>(xlim_data__.data(), - 2); - { - std::vector xlim_flat__; - current_statement__ = 40; - xlim_flat__ = context__.vals_r("xlim"); - current_statement__ = 40; - pos__ = 1; - current_statement__ = 40; - for (int sym1__ = 1; sym1__ <= 2; ++sym1__) { - current_statement__ = 40; - stan::model::assign(xlim, xlim_flat__[(pos__ - 1)], - "assigning variable xlim", stan::model::index_uni(sym1__)); - current_statement__ = 40; - pos__ = (pos__ + 1); - } - } - current_statement__ = 41; - context__.validate_dims("data initialization", "ylim", "double", - std::vector{static_cast(2)}); - ylim_data__ = Eigen::Matrix::Constant(2, - std::numeric_limits::quiet_NaN()); - new (&ylim) Eigen::Map>(ylim_data__.data(), - 2); - { - std::vector ylim_flat__; - current_statement__ = 41; - ylim_flat__ = context__.vals_r("ylim"); - current_statement__ = 41; - pos__ = 1; - current_statement__ = 41; - for (int sym1__ = 1; sym1__ <= 2; ++sym1__) { - current_statement__ = 41; - stan::model::assign(ylim, ylim_flat__[(pos__ - 1)], - "assigning variable ylim", stan::model::index_uni(sym1__)); - current_statement__ = 41; - pos__ = (pos__ + 1); - } - } - current_statement__ = 42; - stan::math::validate_non_negative_index("testX", "ntest", ntest); - current_statement__ = 43; - context__.validate_dims("data initialization", "testX", "double", - std::vector{static_cast(ntest)}); - testX_data__ = Eigen::Matrix::Constant(ntest, - std::numeric_limits::quiet_NaN()); - new (&testX) - Eigen::Map>(testX_data__.data(), ntest); - { - std::vector testX_flat__; - current_statement__ = 43; - testX_flat__ = context__.vals_r("testX"); - current_statement__ = 43; - pos__ = 1; - current_statement__ = 43; - for (int sym1__ = 1; sym1__ <= ntest; ++sym1__) { - current_statement__ = 43; - stan::model::assign(testX, testX_flat__[(pos__ - 1)], - "assigning variable testX", stan::model::index_uni(sym1__)); - current_statement__ = 43; - pos__ = (pos__ + 1); - } - } - current_statement__ = 44; - stan::math::validate_non_negative_index("testY", "ntest", ntest); - current_statement__ = 45; - context__.validate_dims("data initialization", "testY", "double", - std::vector{static_cast(ntest)}); - testY_data__ = Eigen::Matrix::Constant(ntest, - std::numeric_limits::quiet_NaN()); - new (&testY) - Eigen::Map>(testY_data__.data(), ntest); - { - std::vector testY_flat__; - current_statement__ = 45; - testY_flat__ = context__.vals_r("testY"); - current_statement__ = 45; - pos__ = 1; - current_statement__ = 45; - for (int sym1__ = 1; sym1__ <= ntest; ++sym1__) { - current_statement__ = 45; - stan::model::assign(testY, testY_flat__[(pos__ - 1)], - "assigning variable testY", stan::model::index_uni(sym1__)); - current_statement__ = 45; - pos__ = (pos__ + 1); - } - } - current_statement__ = 46; - stan::math::validate_non_negative_index("td2", "ntest", ntest); - current_statement__ = 47; - stan::math::validate_non_negative_index("td2", "nrec", nrec); - current_statement__ = 48; - td2_data__ = Eigen::Matrix::Constant(ntest, nrec, - std::numeric_limits::quiet_NaN()); - new (&td2) Eigen::Map>(td2_data__.data(), - ntest, nrec); - current_statement__ = 51; - for (int s = 1; s <= ntest; ++s) { - current_statement__ = 49; - stan::model::assign(td2, - stan::math::to_row_vector( - stan::math::add( - stan::math::square( - stan::math::subtract(recX, - stan::model::rvalue(testX, "testX", - stan::model::index_uni(s)))), - stan::math::square( - stan::math::subtract(recY, - stan::model::rvalue(testY, "testY", - stan::model::index_uni(s)))))), "assigning variable td2", - stan::model::index_uni(s), stan::model::index_omni()); - } - current_statement__ = 52; - stan::math::validate_non_negative_index("alpha0", "ntime", ntime); - current_statement__ = 53; - stan::math::validate_non_negative_index("alpha0", "nrec", nrec); - current_statement__ = 54; - stan::math::validate_non_negative_index("sx", "nind", nind); - current_statement__ = 55; - stan::math::validate_non_negative_index("sx", "ntime", ntime); - current_statement__ = 56; - stan::math::validate_non_negative_index("sy", "nind", nind); - current_statement__ = 57; - stan::math::validate_non_negative_index("sy", "ntime", ntime); - current_statement__ = 58; - stan::math::validate_non_negative_index("p0", "ntime", ntime); - current_statement__ = 59; - stan::math::validate_non_negative_index("p0", "nrec", nrec); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - num_params_r__ = (ntime * nrec) + 1 + (nind * ntime) + (nind * ntime); - } - inline std::string model_name() const final { - return "model_COA_Tag_Integrated"; - } - inline std::vector model_compile_info() const noexcept { - return std::vector{"stanc_version = stanc3 v2.32.2", - "stancflags = --allow-undefined"}; - } - template * = nullptr, - stan::require_vector_like_vt* = nullptr> - inline stan::scalar_type_t - log_prob_impl(VecR& params_r__, VecI& params_i__, std::ostream* - pstream__ = nullptr) const { - using T__ = stan::scalar_type_t; - using local_scalar_t__ = T__; - T__ lp__(0.0); - stan::math::accumulator lp_accum__; - stan::io::deserializer in__(params_r__, params_i__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - static constexpr const char* function__ = - "model_COA_Tag_Integrated_namespace::log_prob"; - // suppress unused var warning - (void) function__; - try { - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - DUMMY_VAR__); - current_statement__ = 1; - alpha0 = in__.template read_constrain_lub< - Eigen::Matrix, jacobian__>(-5, 5, - lp__, ntime, nrec); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = in__.template read_constrain_lb(0, lp__); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 3; - sx = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), - lp__, nind, ntime); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 4; - sy = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), - lp__, nind, ntime); - { - current_statement__ = 7; - lp_accum__.add(stan::math::cauchy_lpdf( - stan::math::to_vector(alpha0), 0, 2.5)); - current_statement__ = 8; - lp_accum__.add(stan::math::cauchy_lpdf(alpha1, 0, 2.5)); - current_statement__ = 15; - for (int s = 1; s <= ntest; ++s) { - current_statement__ = 9; - stan::math::validate_non_negative_index("dist_decay", "nrec", nrec); - Eigen::Matrix dist_decay = - Eigen::Matrix::Constant(nrec, DUMMY_VAR__); - current_statement__ = 10; - stan::model::assign(dist_decay, - stan::math::multiply(alpha1, - stan::math::to_vector( - stan::model::rvalue(td2, "td2", stan::model::index_uni(s), - stan::model::index_omni()))), - "assigning variable dist_decay"); - current_statement__ = 13; - for (int t = 1; t <= ntime; ++t) { - current_statement__ = 11; - lp_accum__.add(stan::math::binomial_logit_lpmf( - stan::model::rvalue(test, "test", - stan::model::index_uni(s), - stan::model::index_uni(t)), ntrans, - stan::math::subtract( - stan::math::transpose( - stan::math::row(alpha0, t)), dist_decay))); - } - } - current_statement__ = 22; - for (int i = 1; i <= nind; ++i) { - current_statement__ = 20; - for (int t = 1; t <= ntime; ++t) { - current_statement__ = 16; - stan::math::validate_non_negative_index("d2", "nrec", nrec); - Eigen::Matrix d2 = - Eigen::Matrix::Constant(nrec, - DUMMY_VAR__); - current_statement__ = 17; - stan::model::assign(d2, - stan::math::add( - stan::math::square( - stan::math::subtract(recX, - stan::model::rvalue(sx, "sx", stan::model::index_uni(i), - stan::model::index_uni(t)))), - stan::math::square( - stan::math::subtract(recY, - stan::model::rvalue(sy, "sy", stan::model::index_uni(i), - stan::model::index_uni(t))))), "assigning variable d2"); - current_statement__ = 18; - lp_accum__.add(stan::math::binomial_logit_lpmf( - stan::model::rvalue(y, "y", - stan::model::index_uni(i), - stan::model::index_uni(t)), ntrans, - stan::math::subtract( - stan::math::transpose( - stan::math::row(alpha0, t)), - stan::math::multiply(alpha1, d2)))); - } - } - } - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - lp_accum__.add(lp__); - return lp_accum__.sum(); - } - template * = nullptr, stan::require_vector_like_vt* = nullptr, stan::require_vector_vt* = nullptr> - inline void - write_array_impl(RNG& base_rng__, VecR& params_r__, VecI& params_i__, - VecVar& vars__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true, std::ostream* - pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::deserializer in__(params_r__, params_i__); - stan::io::serializer out__(vars__); - static constexpr bool propto__ = true; - // suppress unused var warning - (void) propto__; - double lp__ = 0.0; - // suppress unused var warning - (void) lp__; - int current_statement__ = 0; - stan::math::accumulator lp_accum__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - constexpr bool jacobian__ = false; - static constexpr const char* function__ = - "model_COA_Tag_Integrated_namespace::write_array"; - // suppress unused var warning - (void) function__; - try { - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - std::numeric_limits::quiet_NaN()); - current_statement__ = 1; - alpha0 = in__.template read_constrain_lub< - Eigen::Matrix, jacobian__>(-5, 5, - lp__, ntime, nrec); - double alpha1 = std::numeric_limits::quiet_NaN(); - current_statement__ = 2; - alpha1 = in__.template read_constrain_lb(0, lp__); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - std::numeric_limits::quiet_NaN()); - current_statement__ = 3; - sx = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), - lp__, nind, ntime); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - std::numeric_limits::quiet_NaN()); - current_statement__ = 4; - sy = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), - lp__, nind, ntime); - out__.write(alpha0); - out__.write(alpha1); - out__.write(sx); - out__.write(sy); - if (stan::math::logical_negation( - (stan::math::primitive_value(emit_transformed_parameters__) || - stan::math::primitive_value(emit_generated_quantities__)))) { - return ; - } - if (stan::math::logical_negation(emit_generated_quantities__)) { - return ; - } - Eigen::Matrix p0 = - Eigen::Matrix::Constant(ntime, nrec, - std::numeric_limits::quiet_NaN()); - current_statement__ = 5; - stan::model::assign(p0, stan::math::inv_logit(alpha0), - "assigning variable p0"); - double sigma = std::numeric_limits::quiet_NaN(); - current_statement__ = 6; - sigma = stan::math::sqrt((1.0 / (2.0 * alpha1))); - out__.write(p0); - out__.write(sigma); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - template * = nullptr, - stan::require_vector_like_vt* = nullptr> - inline void - unconstrain_array_impl(const VecVar& params_r__, const VecI& params_i__, - VecVar& vars__, std::ostream* pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::deserializer in__(params_r__, params_i__); - stan::io::serializer out__(vars__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - int pos__ = std::numeric_limits::min(); - pos__ = 1; - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - DUMMY_VAR__); - current_statement__ = 1; - stan::model::assign(alpha0, - in__.read>(ntime, nrec), - "assigning variable alpha0"); - out__.write_free_lub(-5, 5, alpha0); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = in__.read(); - out__.write_free_lb(0, alpha1); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 3; - stan::model::assign(sx, - in__.read>(nind, ntime), - "assigning variable sx"); - out__.write_free_lub(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), sx); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 4; - stan::model::assign(sy, - in__.read>(nind, ntime), - "assigning variable sy"); - out__.write_free_lub(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), sy); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - template * = nullptr> - inline void - transform_inits_impl(const stan::io::var_context& context__, VecVar& - vars__, std::ostream* pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::serializer out__(vars__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - current_statement__ = 1; - context__.validate_dims("parameter initialization", "alpha0", "double", - std::vector{static_cast(ntime), - static_cast(nrec)}); - current_statement__ = 2; - context__.validate_dims("parameter initialization", "alpha1", "double", - std::vector{}); - current_statement__ = 3; - context__.validate_dims("parameter initialization", "sx", "double", - std::vector{static_cast(nind), - static_cast(ntime)}); - current_statement__ = 4; - context__.validate_dims("parameter initialization", "sy", "double", - std::vector{static_cast(nind), - static_cast(ntime)}); - int pos__ = std::numeric_limits::min(); - pos__ = 1; - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - DUMMY_VAR__); - { - std::vector alpha0_flat__; - current_statement__ = 1; - alpha0_flat__ = context__.vals_r("alpha0"); - current_statement__ = 1; - pos__ = 1; - current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 1; - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - current_statement__ = 1; - stan::model::assign(alpha0, alpha0_flat__[(pos__ - 1)], - "assigning variable alpha0", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 1; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(-5, 5, alpha0); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = context__.vals_r("alpha1")[(1 - 1)]; - out__.write_free_lb(0, alpha1); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - { - std::vector sx_flat__; - current_statement__ = 3; - sx_flat__ = context__.vals_r("sx"); - current_statement__ = 3; - pos__ = 1; - current_statement__ = 3; - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - current_statement__ = 3; - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - current_statement__ = 3; - stan::model::assign(sx, sx_flat__[(pos__ - 1)], - "assigning variable sx", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 3; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), sx); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - { - std::vector sy_flat__; - current_statement__ = 4; - sy_flat__ = context__.vals_r("sy"); - current_statement__ = 4; - pos__ = 1; - current_statement__ = 4; - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - current_statement__ = 4; - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - current_statement__ = 4; - stan::model::assign(sy, sy_flat__[(pos__ - 1)], - "assigning variable sy", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 4; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), sy); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - inline void - get_param_names(std::vector& names__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true) const { - names__ = std::vector{"alpha0", "alpha1", "sx", "sy"}; - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - std::vector temp{"p0", "sigma"}; - names__.reserve(names__.size() + temp.size()); - names__.insert(names__.end(), temp.begin(), temp.end()); - } - } - inline void - get_dims(std::vector>& dimss__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true) const { - dimss__ = std::vector>{std::vector{static_cast< - size_t>( - ntime), - static_cast(nrec)}, - std::vector{}, - std::vector{static_cast(nind), - static_cast(ntime)}, - std::vector{static_cast(nind), - static_cast(ntime)}}; - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - std::vector> - temp{std::vector{static_cast(ntime), - static_cast(nrec)}, std::vector{}}; - dimss__.reserve(dimss__.size() + temp.size()); - dimss__.insert(dimss__.end(), temp.begin(), temp.end()); - } - } - inline void - constrained_param_names(std::vector& param_names__, bool - emit_transformed_parameters__ = true, bool - emit_generated_quantities__ = true) const final { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "alpha0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "alpha1"); - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sx" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sy" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "p0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "sigma"); - } - } - inline void - unconstrained_param_names(std::vector& param_names__, bool - emit_transformed_parameters__ = true, bool - emit_generated_quantities__ = true) const final { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "alpha0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "alpha1"); - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sx" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sy" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "p0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "sigma"); - } - } - inline std::string get_constrained_sizedtypes() const { - return std::string("[{\"name\":\"alpha0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"parameters\"},{\"name\":\"alpha1\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"sx\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"sy\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"p0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"generated_quantities\"},{\"name\":\"sigma\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]"); - } - inline std::string get_unconstrained_sizedtypes() const { - return std::string("[{\"name\":\"alpha0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"parameters\"},{\"name\":\"alpha1\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"sx\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"sy\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"p0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"generated_quantities\"},{\"name\":\"sigma\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]"); - } - // Begin method overload boilerplate - template inline void - write_array(RNG& base_rng, Eigen::Matrix& params_r, - Eigen::Matrix& vars, const bool - emit_transformed_parameters = true, const bool - emit_generated_quantities = true, std::ostream* - pstream = nullptr) const { - const size_t num_params__ = ((((ntime * nrec) + 1) + (nind * ntime)) + - (nind * ntime)); - const size_t num_transformed = emit_transformed_parameters * (0); - const size_t num_gen_quantities = emit_generated_quantities * (((ntime * - nrec) + 1)); - const size_t num_to_write = num_params__ + num_transformed + - num_gen_quantities; - std::vector params_i; - vars = Eigen::Matrix::Constant(num_to_write, - std::numeric_limits::quiet_NaN()); - write_array_impl(base_rng, params_r, params_i, vars, - emit_transformed_parameters, emit_generated_quantities, pstream); - } - template inline void - write_array(RNG& base_rng, std::vector& params_r, std::vector& - params_i, std::vector& vars, bool - emit_transformed_parameters = true, bool - emit_generated_quantities = true, std::ostream* - pstream = nullptr) const { - const size_t num_params__ = ((((ntime * nrec) + 1) + (nind * ntime)) + - (nind * ntime)); - const size_t num_transformed = emit_transformed_parameters * (0); - const size_t num_gen_quantities = emit_generated_quantities * (((ntime * - nrec) + 1)); - const size_t num_to_write = num_params__ + num_transformed + - num_gen_quantities; - vars = std::vector(num_to_write, - std::numeric_limits::quiet_NaN()); - write_array_impl(base_rng, params_r, params_i, vars, - emit_transformed_parameters, emit_generated_quantities, pstream); - } - template inline T_ - log_prob(Eigen::Matrix& params_r, std::ostream* pstream = nullptr) const { - Eigen::Matrix params_i; - return log_prob_impl(params_r, params_i, pstream); - } - template inline T_ - log_prob(std::vector& params_r, std::vector& params_i, - std::ostream* pstream = nullptr) const { - return log_prob_impl(params_r, params_i, pstream); - } - inline void - transform_inits(const stan::io::var_context& context, - Eigen::Matrix& params_r, std::ostream* - pstream = nullptr) const final { - std::vector params_r_vec(params_r.size()); - std::vector params_i; - transform_inits(context, params_i, params_r_vec, pstream); - params_r = Eigen::Map>(params_r_vec.data(), - params_r_vec.size()); - } - inline void - transform_inits(const stan::io::var_context& context, std::vector& - params_i, std::vector& vars, std::ostream* - pstream__ = nullptr) const { - vars.resize(num_params_r__); - transform_inits_impl(context, vars, pstream__); - } - inline void - unconstrain_array(const std::vector& params_constrained, - std::vector& params_unconstrained, std::ostream* - pstream = nullptr) const { - const std::vector params_i; - params_unconstrained = std::vector(num_params_r__, - std::numeric_limits::quiet_NaN()); - unconstrain_array_impl(params_constrained, params_i, - params_unconstrained, pstream); - } - inline void - unconstrain_array(const Eigen::Matrix& params_constrained, - Eigen::Matrix& params_unconstrained, - std::ostream* pstream = nullptr) const { - const std::vector params_i; - params_unconstrained = Eigen::Matrix::Constant(num_params_r__, - std::numeric_limits::quiet_NaN()); - unconstrain_array_impl(params_constrained, params_i, - params_unconstrained, pstream); - } -}; -} -using stan_model = model_COA_Tag_Integrated_namespace::model_COA_Tag_Integrated; -#ifndef USING_R -// Boilerplate -stan::model::model_base& -new_model(stan::io::var_context& data_context, unsigned int seed, - std::ostream* msg_stream) { - stan_model* m = new stan_model(data_context, seed, msg_stream); - return *m; -} -stan::math::profile_map& get_stan_profile_data() { - return model_COA_Tag_Integrated_namespace::profiles__; -} -#endif -#endif diff --git a/src/stanExports_COA_Tag_Integrated.o b/src/stanExports_COA_Tag_Integrated.o deleted file mode 100644 index 70f285b..0000000 Binary files a/src/stanExports_COA_Tag_Integrated.o and /dev/null differ diff --git a/src/stanExports_COA_TimeVarying.cc b/src/stanExports_COA_TimeVarying.cc deleted file mode 100644 index a9476ab..0000000 --- a/src/stanExports_COA_TimeVarying.cc +++ /dev/null @@ -1,32 +0,0 @@ -// Generated by rstantools. Do not edit by hand. - -#include -using namespace Rcpp ; -#include "stanExports_COA_TimeVarying.h" - -RCPP_MODULE(stan_fit4COA_TimeVarying_mod) { - - - class_ >("rstantools_model_COA_TimeVarying") - - .constructor() - - - .method("call_sampler", &rstan::stan_fit ::call_sampler) - .method("param_names", &rstan::stan_fit ::param_names) - .method("param_names_oi", &rstan::stan_fit ::param_names_oi) - .method("param_fnames_oi", &rstan::stan_fit ::param_fnames_oi) - .method("param_dims", &rstan::stan_fit ::param_dims) - .method("param_dims_oi", &rstan::stan_fit ::param_dims_oi) - .method("update_param_oi", &rstan::stan_fit ::update_param_oi) - .method("param_oi_tidx", &rstan::stan_fit ::param_oi_tidx) - .method("grad_log_prob", &rstan::stan_fit ::grad_log_prob) - .method("log_prob", &rstan::stan_fit ::log_prob) - .method("unconstrain_pars", &rstan::stan_fit ::unconstrain_pars) - .method("constrain_pars", &rstan::stan_fit ::constrain_pars) - .method("num_pars_unconstrained", &rstan::stan_fit ::num_pars_unconstrained) - .method("unconstrained_param_names", &rstan::stan_fit ::unconstrained_param_names) - .method("constrained_param_names", &rstan::stan_fit ::constrained_param_names) - .method("standalone_gqs", &rstan::stan_fit ::standalone_gqs) - ; -} diff --git a/src/stanExports_COA_TimeVarying.h b/src/stanExports_COA_TimeVarying.h deleted file mode 100644 index 4ceb802..0000000 --- a/src/stanExports_COA_TimeVarying.h +++ /dev/null @@ -1,849 +0,0 @@ -// Generated by rstantools. Do not edit by hand. - -/* - TelemetrySpace is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - TelemetrySpace is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with TelemetrySpace. If not, see . -*/ -#ifndef MODELS_HPP -#define MODELS_HPP -#define STAN__SERVICES__COMMAND_HPP -#ifndef USE_STANC3 -#define USE_STANC3 -#endif -#include -// Code generated by stanc v2.32.2 -#include -namespace model_COA_TimeVarying_namespace { -using stan::model::model_base_crtp; -using namespace stan::math; -stan::math::profile_map profiles__; -static constexpr std::array locations_array__ = - {" (found before start of program)", - " (in 'COA_TimeVarying', line 16, column 2 to column 52)", - " (in 'COA_TimeVarying', line 17, column 2 to column 25)", - " (in 'COA_TimeVarying', line 20, column 2 to column 59)", - " (in 'COA_TimeVarying', line 22, column 2 to column 59)", - " (in 'COA_TimeVarying', line 42, column 2 to column 45)", - " (in 'COA_TimeVarying', line 44, column 2 to column 42)", - " (in 'COA_TimeVarying', line 27, column 2 to column 37)", - " (in 'COA_TimeVarying', line 28, column 2 to column 26)", - " (in 'COA_TimeVarying', line 33, column 13 to column 17)", - " (in 'COA_TimeVarying', line 33, column 6 to column 74)", - " (in 'COA_TimeVarying', line 36, column 6 to column 72)", - " (in 'COA_TimeVarying', line 31, column 23 to line 37, column 5)", - " (in 'COA_TimeVarying', line 31, column 4 to line 37, column 5)", - " (in 'COA_TimeVarying', line 30, column 20 to line 38, column 3)", - " (in 'COA_TimeVarying', line 30, column 2 to line 38, column 3)", - " (in 'COA_TimeVarying', line 3, column 2 to column 22)", - " (in 'COA_TimeVarying', line 4, column 2 to column 22)", - " (in 'COA_TimeVarying', line 5, column 2 to column 23)", - " (in 'COA_TimeVarying', line 6, column 2 to column 24)", - " (in 'COA_TimeVarying', line 7, column 8 to column 12)", - " (in 'COA_TimeVarying', line 7, column 14 to column 19)", - " (in 'COA_TimeVarying', line 7, column 21 to column 25)", - " (in 'COA_TimeVarying', line 7, column 2 to column 44)", - " (in 'COA_TimeVarying', line 8, column 9 to column 13)", - " (in 'COA_TimeVarying', line 8, column 2 to column 20)", - " (in 'COA_TimeVarying', line 9, column 9 to column 13)", - " (in 'COA_TimeVarying', line 9, column 2 to column 20)", - " (in 'COA_TimeVarying', line 10, column 2 to column 17)", - " (in 'COA_TimeVarying', line 11, column 2 to column 17)", - " (in 'COA_TimeVarying', line 16, column 32 to column 37)", - " (in 'COA_TimeVarying', line 16, column 39 to column 43)", - " (in 'COA_TimeVarying', line 20, column 43 to column 47)", - " (in 'COA_TimeVarying', line 20, column 49 to column 54)", - " (in 'COA_TimeVarying', line 22, column 43 to column 47)", - " (in 'COA_TimeVarying', line 22, column 49 to column 54)", - " (in 'COA_TimeVarying', line 42, column 9 to column 14)", - " (in 'COA_TimeVarying', line 42, column 16 to column 20)"}; -#include -class model_COA_TimeVarying final : public model_base_crtp { -private: - int nind; - int nrec; - int ntime; - int ntrans; - std::vector>> y; - Eigen::Matrix recX_data__; - Eigen::Matrix recY_data__; - Eigen::Matrix xlim_data__; - Eigen::Matrix ylim_data__; - Eigen::Map> recX{nullptr, 0}; - Eigen::Map> recY{nullptr, 0}; - Eigen::Map> xlim{nullptr, 0}; - Eigen::Map> ylim{nullptr, 0}; -public: - ~model_COA_TimeVarying() {} - model_COA_TimeVarying(stan::io::var_context& context__, unsigned int - random_seed__ = 0, std::ostream* pstream__ = nullptr) - : model_base_crtp(0) { - int current_statement__ = 0; - using local_scalar_t__ = double; - boost::ecuyer1988 base_rng__ = - stan::services::util::create_rng(random_seed__, 0); - // suppress unused var warning - (void) base_rng__; - static constexpr const char* function__ = - "model_COA_TimeVarying_namespace::model_COA_TimeVarying"; - // suppress unused var warning - (void) function__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - int pos__ = std::numeric_limits::min(); - pos__ = 1; - current_statement__ = 16; - context__.validate_dims("data initialization", "nind", "int", - std::vector{}); - nind = std::numeric_limits::min(); - current_statement__ = 16; - nind = context__.vals_i("nind")[(1 - 1)]; - current_statement__ = 16; - stan::math::check_greater_or_equal(function__, "nind", nind, 0); - current_statement__ = 17; - context__.validate_dims("data initialization", "nrec", "int", - std::vector{}); - nrec = std::numeric_limits::min(); - current_statement__ = 17; - nrec = context__.vals_i("nrec")[(1 - 1)]; - current_statement__ = 17; - stan::math::check_greater_or_equal(function__, "nrec", nrec, 0); - current_statement__ = 18; - context__.validate_dims("data initialization", "ntime", "int", - std::vector{}); - ntime = std::numeric_limits::min(); - current_statement__ = 18; - ntime = context__.vals_i("ntime")[(1 - 1)]; - current_statement__ = 18; - stan::math::check_greater_or_equal(function__, "ntime", ntime, 0); - current_statement__ = 19; - context__.validate_dims("data initialization", "ntrans", "int", - std::vector{}); - ntrans = std::numeric_limits::min(); - current_statement__ = 19; - ntrans = context__.vals_i("ntrans")[(1 - 1)]; - current_statement__ = 19; - stan::math::check_greater_or_equal(function__, "ntrans", ntrans, 0); - current_statement__ = 20; - stan::math::validate_non_negative_index("y", "nind", nind); - current_statement__ = 21; - stan::math::validate_non_negative_index("y", "ntime", ntime); - current_statement__ = 22; - stan::math::validate_non_negative_index("y", "nrec", nrec); - current_statement__ = 23; - context__.validate_dims("data initialization", "y", "int", - std::vector{static_cast(nind), - static_cast(ntime), static_cast(nrec)}); - y = std::vector>>(nind, - std::vector>(ntime, - std::vector(nrec, std::numeric_limits::min()))); - { - std::vector y_flat__; - current_statement__ = 23; - y_flat__ = context__.vals_i("y"); - current_statement__ = 23; - pos__ = 1; - current_statement__ = 23; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 23; - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - current_statement__ = 23; - for (int sym3__ = 1; sym3__ <= nind; ++sym3__) { - current_statement__ = 23; - stan::model::assign(y, y_flat__[(pos__ - 1)], - "assigning variable y", stan::model::index_uni(sym3__), - stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 23; - pos__ = (pos__ + 1); - } - } - } - } - current_statement__ = 23; - stan::math::check_greater_or_equal(function__, "y", y, 0); - current_statement__ = 24; - stan::math::validate_non_negative_index("recX", "nrec", nrec); - current_statement__ = 25; - context__.validate_dims("data initialization", "recX", "double", - std::vector{static_cast(nrec)}); - recX_data__ = Eigen::Matrix::Constant(nrec, - std::numeric_limits::quiet_NaN()); - new (&recX) Eigen::Map>(recX_data__.data(), - nrec); - { - std::vector recX_flat__; - current_statement__ = 25; - recX_flat__ = context__.vals_r("recX"); - current_statement__ = 25; - pos__ = 1; - current_statement__ = 25; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 25; - stan::model::assign(recX, recX_flat__[(pos__ - 1)], - "assigning variable recX", stan::model::index_uni(sym1__)); - current_statement__ = 25; - pos__ = (pos__ + 1); - } - } - current_statement__ = 26; - stan::math::validate_non_negative_index("recY", "nrec", nrec); - current_statement__ = 27; - context__.validate_dims("data initialization", "recY", "double", - std::vector{static_cast(nrec)}); - recY_data__ = Eigen::Matrix::Constant(nrec, - std::numeric_limits::quiet_NaN()); - new (&recY) Eigen::Map>(recY_data__.data(), - nrec); - { - std::vector recY_flat__; - current_statement__ = 27; - recY_flat__ = context__.vals_r("recY"); - current_statement__ = 27; - pos__ = 1; - current_statement__ = 27; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 27; - stan::model::assign(recY, recY_flat__[(pos__ - 1)], - "assigning variable recY", stan::model::index_uni(sym1__)); - current_statement__ = 27; - pos__ = (pos__ + 1); - } - } - current_statement__ = 28; - context__.validate_dims("data initialization", "xlim", "double", - std::vector{static_cast(2)}); - xlim_data__ = Eigen::Matrix::Constant(2, - std::numeric_limits::quiet_NaN()); - new (&xlim) Eigen::Map>(xlim_data__.data(), - 2); - { - std::vector xlim_flat__; - current_statement__ = 28; - xlim_flat__ = context__.vals_r("xlim"); - current_statement__ = 28; - pos__ = 1; - current_statement__ = 28; - for (int sym1__ = 1; sym1__ <= 2; ++sym1__) { - current_statement__ = 28; - stan::model::assign(xlim, xlim_flat__[(pos__ - 1)], - "assigning variable xlim", stan::model::index_uni(sym1__)); - current_statement__ = 28; - pos__ = (pos__ + 1); - } - } - current_statement__ = 29; - context__.validate_dims("data initialization", "ylim", "double", - std::vector{static_cast(2)}); - ylim_data__ = Eigen::Matrix::Constant(2, - std::numeric_limits::quiet_NaN()); - new (&ylim) Eigen::Map>(ylim_data__.data(), - 2); - { - std::vector ylim_flat__; - current_statement__ = 29; - ylim_flat__ = context__.vals_r("ylim"); - current_statement__ = 29; - pos__ = 1; - current_statement__ = 29; - for (int sym1__ = 1; sym1__ <= 2; ++sym1__) { - current_statement__ = 29; - stan::model::assign(ylim, ylim_flat__[(pos__ - 1)], - "assigning variable ylim", stan::model::index_uni(sym1__)); - current_statement__ = 29; - pos__ = (pos__ + 1); - } - } - current_statement__ = 30; - stan::math::validate_non_negative_index("alpha0", "ntime", ntime); - current_statement__ = 31; - stan::math::validate_non_negative_index("alpha0", "nrec", nrec); - current_statement__ = 32; - stan::math::validate_non_negative_index("sx", "nind", nind); - current_statement__ = 33; - stan::math::validate_non_negative_index("sx", "ntime", ntime); - current_statement__ = 34; - stan::math::validate_non_negative_index("sy", "nind", nind); - current_statement__ = 35; - stan::math::validate_non_negative_index("sy", "ntime", ntime); - current_statement__ = 36; - stan::math::validate_non_negative_index("p0", "ntime", ntime); - current_statement__ = 37; - stan::math::validate_non_negative_index("p0", "nrec", nrec); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - num_params_r__ = (ntime * nrec) + 1 + (nind * ntime) + (nind * ntime); - } - inline std::string model_name() const final { - return "model_COA_TimeVarying"; - } - inline std::vector model_compile_info() const noexcept { - return std::vector{"stanc_version = stanc3 v2.32.2", - "stancflags = --allow-undefined"}; - } - template * = nullptr, - stan::require_vector_like_vt* = nullptr> - inline stan::scalar_type_t - log_prob_impl(VecR& params_r__, VecI& params_i__, std::ostream* - pstream__ = nullptr) const { - using T__ = stan::scalar_type_t; - using local_scalar_t__ = T__; - T__ lp__(0.0); - stan::math::accumulator lp_accum__; - stan::io::deserializer in__(params_r__, params_i__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - static constexpr const char* function__ = - "model_COA_TimeVarying_namespace::log_prob"; - // suppress unused var warning - (void) function__; - try { - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - DUMMY_VAR__); - current_statement__ = 1; - alpha0 = in__.template read_constrain_lub< - Eigen::Matrix, jacobian__>(-7, 7, - lp__, ntime, nrec); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = in__.template read_constrain_lb(0, lp__); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 3; - sx = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), - lp__, nind, ntime); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 4; - sy = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), - lp__, nind, ntime); - { - current_statement__ = 7; - lp_accum__.add(stan::math::cauchy_lpdf( - stan::math::to_vector(alpha0), 0, 2.5)); - current_statement__ = 8; - lp_accum__.add(stan::math::cauchy_lpdf(alpha1, 0, 2.5)); - current_statement__ = 15; - for (int i = 1; i <= nind; ++i) { - current_statement__ = 13; - for (int t = 1; t <= ntime; ++t) { - current_statement__ = 9; - stan::math::validate_non_negative_index("d2", "nrec", nrec); - Eigen::Matrix d2 = - Eigen::Matrix::Constant(nrec, - DUMMY_VAR__); - current_statement__ = 10; - stan::model::assign(d2, - stan::math::add( - stan::math::square( - stan::math::subtract(recX, - stan::model::rvalue(sx, "sx", stan::model::index_uni(i), - stan::model::index_uni(t)))), - stan::math::square( - stan::math::subtract(recY, - stan::model::rvalue(sy, "sy", stan::model::index_uni(i), - stan::model::index_uni(t))))), "assigning variable d2"); - current_statement__ = 11; - lp_accum__.add(stan::math::binomial_logit_lpmf( - stan::model::rvalue(y, "y", - stan::model::index_uni(i), - stan::model::index_uni(t)), ntrans, - stan::math::subtract( - stan::math::transpose( - stan::math::row(alpha0, t)), - stan::math::multiply(alpha1, d2)))); - } - } - } - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - lp_accum__.add(lp__); - return lp_accum__.sum(); - } - template * = nullptr, stan::require_vector_like_vt* = nullptr, stan::require_vector_vt* = nullptr> - inline void - write_array_impl(RNG& base_rng__, VecR& params_r__, VecI& params_i__, - VecVar& vars__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true, std::ostream* - pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::deserializer in__(params_r__, params_i__); - stan::io::serializer out__(vars__); - static constexpr bool propto__ = true; - // suppress unused var warning - (void) propto__; - double lp__ = 0.0; - // suppress unused var warning - (void) lp__; - int current_statement__ = 0; - stan::math::accumulator lp_accum__; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - constexpr bool jacobian__ = false; - static constexpr const char* function__ = - "model_COA_TimeVarying_namespace::write_array"; - // suppress unused var warning - (void) function__; - try { - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - std::numeric_limits::quiet_NaN()); - current_statement__ = 1; - alpha0 = in__.template read_constrain_lub< - Eigen::Matrix, jacobian__>(-7, 7, - lp__, ntime, nrec); - double alpha1 = std::numeric_limits::quiet_NaN(); - current_statement__ = 2; - alpha1 = in__.template read_constrain_lb(0, lp__); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - std::numeric_limits::quiet_NaN()); - current_statement__ = 3; - sx = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), - lp__, nind, ntime); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - std::numeric_limits::quiet_NaN()); - current_statement__ = 4; - sy = in__.template read_constrain_lub< - Eigen::Matrix, - jacobian__>(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), - lp__, nind, ntime); - out__.write(alpha0); - out__.write(alpha1); - out__.write(sx); - out__.write(sy); - if (stan::math::logical_negation( - (stan::math::primitive_value(emit_transformed_parameters__) || - stan::math::primitive_value(emit_generated_quantities__)))) { - return ; - } - if (stan::math::logical_negation(emit_generated_quantities__)) { - return ; - } - Eigen::Matrix p0 = - Eigen::Matrix::Constant(ntime, nrec, - std::numeric_limits::quiet_NaN()); - current_statement__ = 5; - stan::model::assign(p0, stan::math::inv_logit(alpha0), - "assigning variable p0"); - double sigma = std::numeric_limits::quiet_NaN(); - current_statement__ = 6; - sigma = stan::math::sqrt((1.0 / (2.0 * alpha1))); - out__.write(p0); - out__.write(sigma); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - template * = nullptr, - stan::require_vector_like_vt* = nullptr> - inline void - unconstrain_array_impl(const VecVar& params_r__, const VecI& params_i__, - VecVar& vars__, std::ostream* pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::deserializer in__(params_r__, params_i__); - stan::io::serializer out__(vars__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - int pos__ = std::numeric_limits::min(); - pos__ = 1; - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - DUMMY_VAR__); - current_statement__ = 1; - stan::model::assign(alpha0, - in__.read>(ntime, nrec), - "assigning variable alpha0"); - out__.write_free_lub(-7, 7, alpha0); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = in__.read(); - out__.write_free_lb(0, alpha1); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 3; - stan::model::assign(sx, - in__.read>(nind, ntime), - "assigning variable sx"); - out__.write_free_lub(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), sx); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - current_statement__ = 4; - stan::model::assign(sy, - in__.read>(nind, ntime), - "assigning variable sy"); - out__.write_free_lub(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), sy); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - template * = nullptr> - inline void - transform_inits_impl(const stan::io::var_context& context__, VecVar& - vars__, std::ostream* pstream__ = nullptr) const { - using local_scalar_t__ = double; - stan::io::serializer out__(vars__); - int current_statement__ = 0; - local_scalar_t__ DUMMY_VAR__(std::numeric_limits::quiet_NaN()); - // suppress unused var warning - (void) DUMMY_VAR__; - try { - current_statement__ = 1; - context__.validate_dims("parameter initialization", "alpha0", "double", - std::vector{static_cast(ntime), - static_cast(nrec)}); - current_statement__ = 2; - context__.validate_dims("parameter initialization", "alpha1", "double", - std::vector{}); - current_statement__ = 3; - context__.validate_dims("parameter initialization", "sx", "double", - std::vector{static_cast(nind), - static_cast(ntime)}); - current_statement__ = 4; - context__.validate_dims("parameter initialization", "sy", "double", - std::vector{static_cast(nind), - static_cast(ntime)}); - int pos__ = std::numeric_limits::min(); - pos__ = 1; - Eigen::Matrix alpha0 = - Eigen::Matrix::Constant(ntime, nrec, - DUMMY_VAR__); - { - std::vector alpha0_flat__; - current_statement__ = 1; - alpha0_flat__ = context__.vals_r("alpha0"); - current_statement__ = 1; - pos__ = 1; - current_statement__ = 1; - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - current_statement__ = 1; - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - current_statement__ = 1; - stan::model::assign(alpha0, alpha0_flat__[(pos__ - 1)], - "assigning variable alpha0", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 1; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(-7, 7, alpha0); - local_scalar_t__ alpha1 = DUMMY_VAR__; - current_statement__ = 2; - alpha1 = context__.vals_r("alpha1")[(1 - 1)]; - out__.write_free_lb(0, alpha1); - Eigen::Matrix sx = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - { - std::vector sx_flat__; - current_statement__ = 3; - sx_flat__ = context__.vals_r("sx"); - current_statement__ = 3; - pos__ = 1; - current_statement__ = 3; - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - current_statement__ = 3; - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - current_statement__ = 3; - stan::model::assign(sx, sx_flat__[(pos__ - 1)], - "assigning variable sx", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 3; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(stan::model::rvalue(xlim, "xlim", - stan::model::index_uni(1)), - stan::model::rvalue(xlim, "xlim", stan::model::index_uni(2)), sx); - Eigen::Matrix sy = - Eigen::Matrix::Constant(nind, ntime, - DUMMY_VAR__); - { - std::vector sy_flat__; - current_statement__ = 4; - sy_flat__ = context__.vals_r("sy"); - current_statement__ = 4; - pos__ = 1; - current_statement__ = 4; - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - current_statement__ = 4; - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - current_statement__ = 4; - stan::model::assign(sy, sy_flat__[(pos__ - 1)], - "assigning variable sy", stan::model::index_uni(sym2__), - stan::model::index_uni(sym1__)); - current_statement__ = 4; - pos__ = (pos__ + 1); - } - } - } - out__.write_free_lub(stan::model::rvalue(ylim, "ylim", - stan::model::index_uni(1)), - stan::model::rvalue(ylim, "ylim", stan::model::index_uni(2)), sy); - } catch (const std::exception& e) { - stan::lang::rethrow_located(e, locations_array__[current_statement__]); - } - } - inline void - get_param_names(std::vector& names__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true) const { - names__ = std::vector{"alpha0", "alpha1", "sx", "sy"}; - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - std::vector temp{"p0", "sigma"}; - names__.reserve(names__.size() + temp.size()); - names__.insert(names__.end(), temp.begin(), temp.end()); - } - } - inline void - get_dims(std::vector>& dimss__, const bool - emit_transformed_parameters__ = true, const bool - emit_generated_quantities__ = true) const { - dimss__ = std::vector>{std::vector{static_cast< - size_t>( - ntime), - static_cast(nrec)}, - std::vector{}, - std::vector{static_cast(nind), - static_cast(ntime)}, - std::vector{static_cast(nind), - static_cast(ntime)}}; - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - std::vector> - temp{std::vector{static_cast(ntime), - static_cast(nrec)}, std::vector{}}; - dimss__.reserve(dimss__.size() + temp.size()); - dimss__.insert(dimss__.end(), temp.begin(), temp.end()); - } - } - inline void - constrained_param_names(std::vector& param_names__, bool - emit_transformed_parameters__ = true, bool - emit_generated_quantities__ = true) const final { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "alpha0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "alpha1"); - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sx" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sy" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "p0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "sigma"); - } - } - inline void - unconstrained_param_names(std::vector& param_names__, bool - emit_transformed_parameters__ = true, bool - emit_generated_quantities__ = true) const final { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "alpha0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "alpha1"); - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sx" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - for (int sym1__ = 1; sym1__ <= ntime; ++sym1__) { - for (int sym2__ = 1; sym2__ <= nind; ++sym2__) { - param_names__.emplace_back(std::string() + "sy" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - if (emit_transformed_parameters__) {} - if (emit_generated_quantities__) { - for (int sym1__ = 1; sym1__ <= nrec; ++sym1__) { - for (int sym2__ = 1; sym2__ <= ntime; ++sym2__) { - param_names__.emplace_back(std::string() + "p0" + '.' + - std::to_string(sym2__) + '.' + std::to_string(sym1__)); - } - } - param_names__.emplace_back(std::string() + "sigma"); - } - } - inline std::string get_constrained_sizedtypes() const { - return std::string("[{\"name\":\"alpha0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"parameters\"},{\"name\":\"alpha1\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"sx\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"sy\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"p0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"generated_quantities\"},{\"name\":\"sigma\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]"); - } - inline std::string get_unconstrained_sizedtypes() const { - return std::string("[{\"name\":\"alpha0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"parameters\"},{\"name\":\"alpha1\",\"type\":{\"name\":\"real\"},\"block\":\"parameters\"},{\"name\":\"sx\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"sy\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(nind) + ",\"cols\":" + std::to_string(ntime) + "},\"block\":\"parameters\"},{\"name\":\"p0\",\"type\":{\"name\":\"matrix\",\"rows\":" + std::to_string(ntime) + ",\"cols\":" + std::to_string(nrec) + "},\"block\":\"generated_quantities\"},{\"name\":\"sigma\",\"type\":{\"name\":\"real\"},\"block\":\"generated_quantities\"}]"); - } - // Begin method overload boilerplate - template inline void - write_array(RNG& base_rng, Eigen::Matrix& params_r, - Eigen::Matrix& vars, const bool - emit_transformed_parameters = true, const bool - emit_generated_quantities = true, std::ostream* - pstream = nullptr) const { - const size_t num_params__ = ((((ntime * nrec) + 1) + (nind * ntime)) + - (nind * ntime)); - const size_t num_transformed = emit_transformed_parameters * (0); - const size_t num_gen_quantities = emit_generated_quantities * (((ntime * - nrec) + 1)); - const size_t num_to_write = num_params__ + num_transformed + - num_gen_quantities; - std::vector params_i; - vars = Eigen::Matrix::Constant(num_to_write, - std::numeric_limits::quiet_NaN()); - write_array_impl(base_rng, params_r, params_i, vars, - emit_transformed_parameters, emit_generated_quantities, pstream); - } - template inline void - write_array(RNG& base_rng, std::vector& params_r, std::vector& - params_i, std::vector& vars, bool - emit_transformed_parameters = true, bool - emit_generated_quantities = true, std::ostream* - pstream = nullptr) const { - const size_t num_params__ = ((((ntime * nrec) + 1) + (nind * ntime)) + - (nind * ntime)); - const size_t num_transformed = emit_transformed_parameters * (0); - const size_t num_gen_quantities = emit_generated_quantities * (((ntime * - nrec) + 1)); - const size_t num_to_write = num_params__ + num_transformed + - num_gen_quantities; - vars = std::vector(num_to_write, - std::numeric_limits::quiet_NaN()); - write_array_impl(base_rng, params_r, params_i, vars, - emit_transformed_parameters, emit_generated_quantities, pstream); - } - template inline T_ - log_prob(Eigen::Matrix& params_r, std::ostream* pstream = nullptr) const { - Eigen::Matrix params_i; - return log_prob_impl(params_r, params_i, pstream); - } - template inline T_ - log_prob(std::vector& params_r, std::vector& params_i, - std::ostream* pstream = nullptr) const { - return log_prob_impl(params_r, params_i, pstream); - } - inline void - transform_inits(const stan::io::var_context& context, - Eigen::Matrix& params_r, std::ostream* - pstream = nullptr) const final { - std::vector params_r_vec(params_r.size()); - std::vector params_i; - transform_inits(context, params_i, params_r_vec, pstream); - params_r = Eigen::Map>(params_r_vec.data(), - params_r_vec.size()); - } - inline void - transform_inits(const stan::io::var_context& context, std::vector& - params_i, std::vector& vars, std::ostream* - pstream__ = nullptr) const { - vars.resize(num_params_r__); - transform_inits_impl(context, vars, pstream__); - } - inline void - unconstrain_array(const std::vector& params_constrained, - std::vector& params_unconstrained, std::ostream* - pstream = nullptr) const { - const std::vector params_i; - params_unconstrained = std::vector(num_params_r__, - std::numeric_limits::quiet_NaN()); - unconstrain_array_impl(params_constrained, params_i, - params_unconstrained, pstream); - } - inline void - unconstrain_array(const Eigen::Matrix& params_constrained, - Eigen::Matrix& params_unconstrained, - std::ostream* pstream = nullptr) const { - const std::vector params_i; - params_unconstrained = Eigen::Matrix::Constant(num_params_r__, - std::numeric_limits::quiet_NaN()); - unconstrain_array_impl(params_constrained, params_i, - params_unconstrained, pstream); - } -}; -} -using stan_model = model_COA_TimeVarying_namespace::model_COA_TimeVarying; -#ifndef USING_R -// Boilerplate -stan::model::model_base& -new_model(stan::io::var_context& data_context, unsigned int seed, - std::ostream* msg_stream) { - stan_model* m = new stan_model(data_context, seed, msg_stream); - return *m; -} -stan::math::profile_map& get_stan_profile_data() { - return model_COA_TimeVarying_namespace::profiles__; -} -#endif -#endif diff --git a/src/stanExports_COA_TimeVarying.o b/src/stanExports_COA_TimeVarying.o deleted file mode 100644 index 960bdf3..0000000 Binary files a/src/stanExports_COA_TimeVarying.o and /dev/null differ diff --git a/tests/testthat/setup-test-env.R b/tests/testthat/setup-test-env.R index adea9dd..fb2f16f 100644 --- a/tests/testthat/setup-test-env.R +++ b/tests/testthat/setup-test-env.R @@ -1,30 +1,31 @@ -nsentinel <- 1 +n_sentinel <- 1 # ----- standata 1 ------ standata <- list( - nind = model_param_ex$nind, # number of individuals - nrec = model_param_ex$nrec, # number of receivers - ntime = model_param_ex$tsteps, # number of time steps - ntrans = model_param_ex$ntrans, # number of expected transmissions per tag per time interval - y = Y, # array of detections - recX = rlocs$east, # E-W receiver coordinates - recY = rlocs$north, # N-S receiver coordinates - xlim = example_extent$xlim, # E-W boundary of spatial extent (receiver array + buffer) - ylim = example_extent$ylim + n_ind = model_param_ex$nind, # number of individuals + n_rec = model_param_ex$nrec, # number of receivers + n_time = model_param_ex$tsteps, # number of time steps + n_trans = model_param_ex$ntrans, # number of expected transmissions per tag per time interval + det = Y, # array of detections + rec_x = rlocs$east, # E-W receiver coordinates + rec_y = rlocs$north, # N-S receiver coordinates + x_lim = example_extent$xlim, # E-W boundary of spatial extent (receiver array + buffer) + y_lim = example_extent$ylim ) -standata_1 <- list( - nind = model_param_ex$nind, # number of individuals - nrec = model_param_ex$nrec, # number of receivers - ntime = model_param_ex$tsteps, # number of time steps - ntrans = model_param_ex$ntrans, # number of expected transmissions per tag per time interval - y = Y, # array of detections - recX = rlocs$east, # E-W receiver coordinates - recY = rlocs$north, # N-S receiver coordinates - xlim = example_extent$xlim, # E-W boundary of spatial extent (receiver array + buffer) - ylim = example_extent$ylim, - ntest = nsentinel, - test = testY, - testX = array(testloc$east, dim = c(nsentinel)), - testY = array(testloc$north, dim = c(nsentinel)) # N-S b +standata_test_tag <- list( + n_ind = model_param_ex$nind, # number of individuals + n_rec = model_param_ex$nrec, # number of receivers + n_time = model_param_ex$tsteps, # number of time steps + n_trans = model_param_ex$ntrans, # number of expected transmissions per tag per time interval + n_trans_test = model_param_ex$ntrans, # number of expected transmissions per tag per time interval + det = Y, # array of detections + rec_x = rlocs$east, # E-W receiver coordinates + rec_y = rlocs$north, # N-S receiver coordinates + x_lim = example_extent$xlim, # E-W boundary of spatial extent (receiver array + buffer) + y_lim = example_extent$ylim, + n_test = n_sentinel, + det_test = testY, + test_x = array(testloc$east, dim = c(n_sentinel)), + test_y = array(testloc$north, dim = c(n_sentinel)) # N-S b ) init_fun <- function() { @@ -44,24 +45,40 @@ init_fun <- function() { # ----- run each model ------ # ----- standard coa ------ -set.seed(8675309) -model_coa_standard <- do.call( +standard_gaussian <- do.call( COA_Standard, c( standata, list( chains = 2, - warmup = 1000, + warmup = 1500, iter = 2000, control = list(adapt_delta = 0.95), seed = 4, + n_draws = 11, + init = init_fun, + decay = "gaussian" + ) + ) +) +standard_logistic <- do.call( + COA_Standard, + c( + standata, + list( + chains = 2, + warmup = 2000, + iter = 2500, + control = list(adapt_delta = 0.95), + seed = 4, ndraws = 11, - init = init_fun + init = init_fun, + decay = "logistic" ) ) ) # ----- time integrated ----- -model_coa_time_vary <- do.call( +time_vary_gaussian <- do.call( COA_TimeVarying, c( standata, @@ -72,16 +89,49 @@ model_coa_time_vary <- do.call( control = list(adapt_delta = 0.95), seed = 4, ndraws = 11, - init = init_fun + init = init_fun, + decay = "gaussian" + ) + ) +) +time_vary_logistic <- do.call( + COA_TimeVarying, + c( + standata, + list( + chains = 2, + warmup = 3000, + iter = 7000, + control = list(adapt_delta = 0.95), + seed = 4, + ndraws = 11, + init = init_fun, + decay = "logistic" ) ) ) # ----- tag integraged ----- -model_coa_tag_int <- do.call( +tag_int_gaussian <- do.call( + COA_TagInt, + c( + standata_test_tag, + list( + chains = 2, + warmup = 4000, + iter = 8000, + control = list(adapt_delta = 0.95), + seed = 4, + ndraws = 11, + init = init_fun, + decay = "gaussian" + ) + ) +) +tag_int_logistic <- do.call( COA_TagInt, c( - standata_1, + standata_test_tag, list( chains = 2, warmup = 4000, @@ -89,7 +139,8 @@ model_coa_tag_int <- do.call( control = list(adapt_delta = 0.95), seed = 4, ndraws = 11, - init = init_fun + init = init_fun, + decay = "logistic" ) ) ) diff --git a/tests/testthat/test-COA_Standard.R b/tests/testthat/test-COA_Standard.R index 3985ea8..792a30a 100644 --- a/tests/testthat/test-COA_Standard.R +++ b/tests/testthat/test-COA_Standard.R @@ -4,15 +4,15 @@ # Base arguments for COA_Standard coa_args <- list( - nind = model_param_ex$nind, - nrec = model_param_ex$nrec, - ntime = model_param_ex$tsteps, - ntrans = model_param_ex$ntrans, - y = Y, - recX = rlocs$east, - recY = rlocs$north, - xlim = example_extent$xlim, - ylim = example_extent$ylim, + n_ind = model_param_ex$nind, + n_rec = model_param_ex$nrec, + n_time = model_param_ex$tsteps, + n_trans = model_param_ex$ntrans, + det = Y, + rec_x = rlocs$east, + rec_y = rlocs$north, + x_lim = example_extent$xlim, + y_lim = example_extent$ylim, chains = 2, warmup = 1000, iter = 2000, @@ -27,54 +27,49 @@ call_coa <- function(overrides) { params_table <- list( list( - param = "nind", + param = "n_ind", bad = list("a", NA, c(1, 2)), - regex = "`nind` must be a numeric vector that has a length of 1." + regex = "`n_ind` must be a numeric vector that has a length of 1." ), list( - param = "nrec", + param = "n_rec", bad = list("a", NA, c(1, 2)), - regex = "`nrec` must be a numeric vector that has a length of 1." + regex = "`n_rec` must be a numeric vector that has a length of 1." ), list( - param = "ntime", + param = "n_time", bad = list("a", NA, c(1, 2)), - regex = "`ntime` must be a numeric vector that has a length of 1." + regex = "`n_time` must be a numeric vector that has a length of 1." ), list( - param = "ntrans", + param = "n_trans", bad = list(c(model_param_ex$ntrans, model_param_ex$ntrans), "1"), - regex = "`ntrans` must be a numeric vector that has a length of 1." + regex = "`n_trans` must be a numeric vector that has a length of 1." ), list( - param = "y", + param = "det", bad = list(c(1, 2, 3), "a"), - regex = "`y` must be a 3-dimensional numeric array." + regex = "`det` must be a 3-dimensional numeric array." ), list( - param = "recX", + param = "rec_x", bad = list("a", NA), - regex = "`recX` must be a numeric vector that has a length of 1." + regex = "`rec_x` must be a numeric vector that has a length of 1." ), list( - param = "recY", + param = "rec_y", bad = list("a", NA), - regex = "`recY` must be a numeric vector that has a length of 1." + regex = "`rec_y` must be a numeric vector that has a length of 1." ), list( - param = "xlim", + param = "x_lim", bad = list("a", c(1, 2, 3)), - regex = "`xlim` must be a numeric vector that has a length of 2." + regex = "`x_lim` must be a numeric vector that has a length of 2." ), list( - param = "ylim", + param = "y_lim", bad = list("a", c(1, 2, 3)), - regex = "`ylim` must be a numeric vector that has a length of 2." - ), - list( - param = "ndraws", - bad = list("a", c(1, 2, 3)), - regex = "`ndraws` must be a numeric vector that has a length of 1." + regex = "`y_lim` must be a numeric vector that has a length of 2." ) ) @@ -109,36 +104,91 @@ test_that("parameter validation works", { # ---- run model and check of it works ---- -# model_coa_standard$generated_quantities +# standard_gaussian$generated_quantities -# bayesplot::ppc_dens_overlay(y = as.vector(Y), yrep = model_coa_standard$generated_quantities) +# bayesplot::ppc_dens_overlay(y = as.vector(Y), y_rep = standard_gaussian$generated_quantities) # rstan::traceplot(fit$model, pars = c("alpha0", "alpha1", # "sigma", "lp__")) -test_that("test COA_standard model results to make sure its consisitent", { - mean_p0 <- model_coa_standard$summary[1] +test_that("test COA_standard gaussian model results to make sure its consistent", { + mean_p0 <- standard_gaussian$summary[1] expected_mean_p0 <- 0.2658 expect_equal(mean_p0, expected_mean_p0, tolerance = 0.05) }) -test_that("check to see if model_coa_standard classes", { - expect_type(model_coa_standard, "list") - expect_s4_class(model_coa_standard$model, "stanfit") - expect_s3_class(model_coa_standard$coas, "data.frame") - expect_s3_class(model_coa_standard$all_estimates, "data.frame") - expect_type(model_coa_standard$summary, "double") - expect_true(is.matrix(model_coa_standard$summary)) - expect_true(is.matrix(model_coa_standard$generated_quantities$yrep)) - expect_type(model_coa_standard$generated_quantities, "list") - expect_true(is.numeric(model_coa_standard$time)) +test_that("check standard_gaussian classes", { + expect_type(standard_gaussian, "list") + expect_s4_class(standard_gaussian$model, "stanfit") + expect_s3_class(standard_gaussian$coas, "data.frame") + expect_s3_class(standard_gaussian$all_estimates, "data.frame") + expect_type(standard_gaussian$summary, "double") + expect_true(is.matrix(standard_gaussian$summary)) + expect_true(is.matrix(standard_gaussian$generated_quantities$y_rep)) + expect_type(standard_gaussian$generated_quantities, "list") + expect_true(is.numeric(standard_gaussian$time)) +}) + + +test_that("check to see if coa returns proper info", { + expect_true("coas" %in% names(standard_gaussian)) + expect_equal(nrow(standard_gaussian$coas), model_param_ex$tsteps) + expect_equal( + colnames(standard_gaussian$coas), + c( + "time", + "x", + "y", + "x_lower", + "x_upper", + "y_lower", + "y_upper" + ) + ) + + for (col in colnames(standard_gaussian$coas)) { + expect_type(standard_gaussian$coas[[col]], "double") + expect_true(all(is.finite(standard_gaussian$coas[[col]]))) + } +}) + +test_that("check to see model converged and has a good rhat", { + rhat <- standard_gaussian$summary[, "Rhat"] + expect_true(all(rhat > 0.95 & rhat < 1.05)) +}) + + +# ----- check if gq retruns the correct length ------ + +test_that("check to see if gq is the correct length", { + expected <- 11 + expect_true(nrow(standard_gaussian$generated_quantities$y_rep) %in% expected) +}) + + +#### LOGISTIC #### +test_that("test COA_standard logistic model results to make sure its consistent", { + mean_p0 <- standard_logistic$summary[1] + expected_mean_p0 <- 0.5849 + expect_equal(mean_p0, expected_mean_p0, tolerance = 0.05) }) +test_that("check standard_logistic classes", { + expect_type(standard_logistic, "list") + expect_s4_class(standard_logistic$model, "stanfit") + expect_s3_class(standard_logistic$coas, "data.frame") + expect_s3_class(standard_logistic$all_estimates, "data.frame") + expect_type(standard_logistic$summary, "double") + expect_true(is.matrix(standard_logistic$summary)) + expect_true(is.matrix(standard_logistic$generated_quantities$y_rep)) + expect_type(standard_logistic$generated_quantities, "list") + expect_true(is.numeric(standard_logistic$time)) +}) test_that("check to see if coa returns proper info", { - expect_true("coas" %in% names(model_coa_standard)) - expect_equal(nrow(model_coa_standard$coas), model_param_ex$tsteps) + expect_true("coas" %in% names(standard_logistic)) + expect_equal(nrow(standard_logistic$coas), model_param_ex$tsteps) expect_equal( - colnames(model_coa_standard$coas), + colnames(standard_logistic$coas), c( "time", "x", @@ -150,14 +200,14 @@ test_that("check to see if coa returns proper info", { ) ) - for (col in colnames(model_coa_standard$coas)) { - expect_type(model_coa_standard$coas[[col]], "double") - expect_true(all(is.finite(model_coa_standard$coas[[col]]))) + for (col in colnames(standard_logistic$coas)) { + expect_type(standard_logistic$coas[[col]], "double") + expect_true(all(is.finite(standard_logistic$coas[[col]]))) } }) test_that("check to see model converged and has a good rhat", { - rhat <- model_coa_standard$summary[, "Rhat"] + rhat <- standard_logistic$summary[, "Rhat"] expect_true(all(rhat > 0.95 & rhat < 1.05)) }) @@ -166,5 +216,5 @@ test_that("check to see model converged and has a good rhat", { test_that("check to see if gq is the correct length", { expected <- 11 - expect_true(nrow(model_coa_standard$generated_quantities$yrep) %in% expected) + expect_true(nrow(standard_logistic$generated_quantities$y_rep) %in% expected) }) diff --git a/tests/testthat/test-COA_Tag_Integrated.R b/tests/testthat/test-COA_Tag_Integrated.R index 55c6308..eaf8faf 100644 --- a/tests/testthat/test-COA_Tag_Integrated.R +++ b/tests/testthat/test-COA_Tag_Integrated.R @@ -1,20 +1,21 @@ -# ----- Model checked from setup-test-env is object model_coa_tag_int----- +# ----- Model checked from setup-test-env is object tag_int_gaussian----- # ---- test each argument if it errors appropriately ----- coa_args <- list( - nind = model_param_ex$nind, - nrec = model_param_ex$nrec, - ntime = model_param_ex$tsteps, - ntrans = model_param_ex$ntrans, - ntest = nsentinel, - y = Y, - test = testY, - recX = rlocs$east, - recY = rlocs$north, - xlim = example_extent$xlim, - ylim = example_extent$ylim, - testX = array(testloc$east, dim = c(nsentinel)), - testY = array(testloc$north, dim = c(nsentinel)), + n_ind = model_param_ex$nind, + n_rec = model_param_ex$nrec, + n_time = model_param_ex$tsteps, + n_trans = model_param_ex$ntrans, + n_trans_test = model_param_ex$ntrans, + n_test = n_sentinel, + det = Y, + det_test = testY, + rec_x = rlocs$east, + rec_y = rlocs$north, + x_lim = example_extent$xlim, + y_lim = example_extent$ylim, + test_x = array(testloc$east, dim = c(n_sentinel)), + test_y = array(testloc$north, dim = c(n_sentinel)), chains = 2, warmup = 1000, iter = 2000, @@ -27,69 +28,74 @@ call_coa_tagint <- function(overrides) { # ----- create param tables params_table <- list( list( - param = "nind", + param = "n_ind", bad = list("bc", NA, c(1, 2)), - regex = "`nind` must be a numeric vector that has a length of 1." + regex = "`n_ind` must be a numeric vector that has a length of 1." ), list( - param = "nrec", + param = "n_rec", bad = list("bc", NA, c(1, 2)), - regex = "`nrec` must be a numeric vector that has a length of 1." + regex = "`n_rec` must be a numeric vector that has a length of 1." ), list( - param = "ntime", + param = "n_time", bad = list("bc", NA, c(1, 2)), - regex = "`ntime` must be a numeric vector that has a length of 1." + regex = "`n_time` must be a numeric vector that has a length of 1." ), list( - param = "ntrans", + param = "n_trans", bad = list(c(model_param_ex$ntrans, model_param_ex$ntrans), "1"), - regex = "`ntrans` must be a numeric vector that has a length of 1." + regex = "`n_trans` must be a numeric vector that has a length of 1." ), list( - param = "ntest", + param = "n_trans_test", + bad = list(c(model_param_ex$ntrans, model_param_ex$ntrans), "1"), + regex = "`n_trans_test` must be a numeric vector that has a length of 1." + ), + list( + param = "n_test", bad = list(c(3, 6, 3), "1"), - regex = "`ntest` must be a numeric vector that has a length of 1." + regex = "`n_test` must be a numeric vector that has a length of 1." ), list( - param = "y", + param = "det", bad = list(c(1, 2, 3), "bc"), - regex = "`y` must be a 3-dimensional numeric array." + regex = "`det` must be a 3-dimensional numeric array." ), list( - param = "test", + param = "det_test", bad = list(c(1, 2, 3), "bc"), - regex = "`test` must be a 3-dimensional numeric array." + regex = "`det_test` must be a 3-dimensional numeric array." ), list( - param = "recX", + param = "rec_x", bad = list("bc", NA), - regex = "`recX` must be a numeric vector that has a length of 1." + regex = "`rec_x` must be a numeric vector that has a length of 1." ), list( - param = "recY", + param = "rec_y", bad = list("bc", NA), - regex = "`recY` must be a numeric vector that has a length of 1." + regex = "`rec_y` must be a numeric vector that has a length of 1." ), list( - param = "xlim", + param = "x_lim", bad = list("bc", c(1, 2, 3)), - regex = "`xlim` must be a numeric vector that has a length of 2." + regex = "`x_lim` must be a numeric vector that has a length of 2." ), list( - param = "ylim", + param = "y_lim", bad = list("bc", c(1, 2, 3)), - regex = "`ylim` must be a numeric vector that has a length of 2." + regex = "`y_lim` must be a numeric vector that has a length of 2." ), list( - param = "testX", + param = "test_x", bad = list("bc", NA), - regex = "`testX` must be a numeric array with length equal to 1 \\(the number of test tags\\)\\." + regex = "`test_x` must be a numeric array with length equal to 1 \\(the number of test tags\\)\\." ), list( - param = "testY", + param = "test_y", bad = list("bc", NA), - regex = "`testY` must be a numeric array with length equal to 1 \\(the number of test tags\\)\\." + regex = "`test_y` must be a numeric array with length equal to 1 \\(the number of test tags\\)\\." ) ) @@ -125,36 +131,36 @@ test_that("parameter validation works", { # ---- run model and check of it works ---- -# rstan::traceplot(model_coa_tag_int$model, pars = c("alpha0", "alpha1", +# rstan::traceplot(tag_int_gaussian$model, pars = c("alpha0", "alpha1", # "sigma", "lp__")) -test_that("test COA_TagInt model results to make sure its consisitent", { - mean_p0 <- model_coa_tag_int$summary[1] +test_that("test COA_TagInt model results to make sure its consistent", { + mean_p0 <- tag_int_gaussian$summary[1] expected_mean_p0 <- 0.5008 expect_equal(mean_p0, expected_mean_p0, tolerance = 0.05) }) -test_that("check to see if model_coa_tag_int classes", { - expect_type(model_coa_tag_int, "list") - expect_s4_class(model_coa_tag_int$model, "stanfit") - expect_s3_class(model_coa_tag_int$coas, "data.frame") - expect_s3_class(model_coa_tag_int$all_estimates, "data.frame") - expect_type(model_coa_tag_int$summary, "double") - expect_true(is.matrix(model_coa_tag_int$summary)) - expect_type(model_coa_tag_int$generated_quantities, "list") - expect_true(is.matrix(model_coa_tag_int$generated_quantities$yrep)) - expect_true(is.matrix(model_coa_tag_int$generated_quantities$testrep)) - expect_true(is.numeric(model_coa_tag_int$time)) +test_that("check tag_int_gaussian classes", { + expect_type(tag_int_gaussian, "list") + expect_s4_class(tag_int_gaussian$model, "stanfit") + expect_s3_class(tag_int_gaussian$coas, "data.frame") + expect_s3_class(tag_int_gaussian$all_estimates, "data.frame") + expect_type(tag_int_gaussian$summary, "double") + expect_true(is.matrix(tag_int_gaussian$summary)) + expect_type(tag_int_gaussian$generated_quantities, "list") + expect_true(is.matrix(tag_int_gaussian$generated_quantities$y_rep)) + expect_true(is.matrix(tag_int_gaussian$generated_quantities$test_rep)) + expect_true(is.numeric(tag_int_gaussian$time)) }) test_that("check to see if coa returns proper info", { - expect_true("coas" %in% names(model_coa_tag_int)) - expect_equal(nrow(model_coa_tag_int$coas), model_param_ex$tsteps) + expect_true("coas" %in% names(tag_int_gaussian)) + expect_equal(nrow(tag_int_gaussian$coas), model_param_ex$tsteps) expect_equal( - colnames(model_coa_tag_int$coas), + colnames(tag_int_gaussian$coas), c( "time", "x", @@ -166,14 +172,14 @@ test_that("check to see if coa returns proper info", { ) ) - for (col in colnames(model_coa_tag_int$coas)) { - expect_type(model_coa_tag_int$coas[[col]], "double") - expect_true(all(is.finite(model_coa_tag_int$coas[[col]]))) + for (col in colnames(tag_int_gaussian$coas)) { + expect_type(tag_int_gaussian$coas[[col]], "double") + expect_true(all(is.finite(tag_int_gaussian$coas[[col]]))) } }) test_that("check to see model converged and has a good rhat", { - rhat <- model_coa_tag_int$summary[, "Rhat"] + rhat <- tag_int_gaussian$summary[, "Rhat"] expect_true(all(rhat > 0.95 & rhat < 1.05)) }) @@ -181,8 +187,63 @@ test_that("check to see model converged and has a good rhat", { test_that("check to see if gq is the correct length", { expected <- 11 - expect_true(nrow(model_coa_tag_int$generated_quantities$yrep) %in% expected) + expect_true(nrow(tag_int_gaussian$generated_quantities$y_rep) %in% expected) expect_true( - nrow(model_coa_tag_int$generated_quantities$testrep) %in% expected + nrow(tag_int_gaussian$generated_quantities$test_rep) %in% expected + ) +}) + + +#### LOGISTIC #### +test_that("test COA_standard logistic model results to make sure its consistent", { + mean_p0 <- tag_int_logistic$summary[1] + expected_mean_p0 <- 0.4899 + expect_equal(mean_p0, expected_mean_p0, tolerance = 0.05) +}) + +test_that("check tag_int_logistic classes", { + expect_type(tag_int_logistic, "list") + expect_s4_class(tag_int_logistic$model, "stanfit") + expect_s3_class(tag_int_logistic$coas, "data.frame") + expect_s3_class(tag_int_logistic$all_estimates, "data.frame") + expect_type(tag_int_logistic$summary, "double") + expect_true(is.matrix(tag_int_logistic$summary)) + expect_true(is.matrix(tag_int_logistic$generated_quantities$y_rep)) + expect_type(tag_int_logistic$generated_quantities, "list") + expect_true(is.numeric(tag_int_logistic$time)) +}) + +test_that("check to see if coa returns proper info", { + expect_true("coas" %in% names(tag_int_logistic)) + expect_equal(nrow(tag_int_logistic$coas), model_param_ex$tsteps) + expect_equal( + colnames(tag_int_logistic$coas), + c( + "time", + "x", + "y", + "x_lower", + "x_upper", + "y_lower", + "y_upper" + ) ) + + for (col in colnames(tag_int_logistic$coas)) { + expect_type(tag_int_logistic$coas[[col]], "double") + expect_true(all(is.finite(tag_int_logistic$coas[[col]]))) + } +}) + +test_that("check to see model converged and has a good rhat", { + rhat <- tag_int_logistic$summary[, "Rhat"] + expect_true(all(rhat > 0.95 & rhat < 1.05)) +}) + + +# ----- check if gq retruns the correct length ------ + +test_that("check to see if gq is the correct length", { + expected <- 11 + expect_true(nrow(tag_int_logistic$generated_quantities$y_rep) %in% expected) }) diff --git a/tests/testthat/test-COA_TimeVarying.R b/tests/testthat/test-COA_TimeVarying.R index 25109af..8212286 100644 --- a/tests/testthat/test-COA_TimeVarying.R +++ b/tests/testthat/test-COA_TimeVarying.R @@ -1,19 +1,19 @@ -# ----- Model checked from setup-test-env is object model_coa_time_vary ----- +# ----- Model checked from setup-test-env is object time_vary_gaussian ----- # ---- test each argument if it errors appropriately ----- -# ---- Check if nind errors ----- +# ---- Check if n_ind errors ----- # Base arguments for COA_Standard coa_args <- list( - nind = model_param_ex$nind, - nrec = model_param_ex$nrec, - ntime = model_param_ex$tsteps, - ntrans = model_param_ex$ntrans, - y = Y, - recX = rlocs$east, - recY = rlocs$north, - xlim = example_extent$xlim, - ylim = example_extent$ylim, + n_ind = model_param_ex$nind, + n_rec = model_param_ex$nrec, + n_time = model_param_ex$tsteps, + n_trans = model_param_ex$ntrans, + det = Y, + rec_x = rlocs$east, + rec_y = rlocs$north, + x_lim = example_extent$xlim, + y_lim = example_extent$ylim, chains = 2, warmup = 1000, iter = 2000, @@ -28,49 +28,49 @@ call_coa_timevarying <- function(overrides) { params_table <- list( list( - param = "nind", + param = "n_ind", bad = list("a", NA, c(1, 2)), - regex = "`nind` must be a numeric vector that has a length of 1." + regex = "`n_ind` must be a numeric vector that has a length of 1." ), list( - param = "nrec", + param = "n_rec", bad = list("a", NA, c(1, 2)), - regex = "`nrec` must be a numeric vector that has a length of 1." + regex = "`n_rec` must be a numeric vector that has a length of 1." ), list( - param = "ntime", + param = "n_time", bad = list("a", NA, c(1, 2)), - regex = "`ntime` must be a numeric vector that has a length of 1." + regex = "`n_time` must be a numeric vector that has a length of 1." ), list( - param = "ntrans", + param = "n_trans", bad = list(c(model_param_ex$ntrans, model_param_ex$ntrans), "1"), - regex = "`ntrans` must be a numeric vector that has a length of 1." + regex = "`n_trans` must be a numeric vector that has a length of 1." ), list( - param = "y", + param = "det", bad = list(c(1, 2, 3), "a"), - regex = "`y` must be a 3-dimensional numeric array." + regex = "`det` must be a 3-dimensional numeric array." ), list( - param = "recX", + param = "rec_x", bad = list("a", NA), - regex = "`recX` must be a numeric vector that has a length of 1." + regex = "`rec_x` must be a numeric vector that has a length of 1." ), list( - param = "recY", + param = "rec_y", bad = list("a", NA), - regex = "`recY` must be a numeric vector that has a length of 1." + regex = "`rec_y` must be a numeric vector that has a length of 1." ), list( - param = "xlim", + param = "x_lim", bad = list("a", c(1, 2, 3)), - regex = "`xlim` must be a numeric vector that has a length of 2." + regex = "`x_lim` must be a numeric vector that has a length of 2." ), list( - param = "ylim", + param = "y_lim", bad = list("a", c(1, 2, 3)), - regex = "`ylim` must be a numeric vector that has a length of 2." + regex = "`y_lim` must be a numeric vector that has a length of 2." ) ) @@ -104,39 +104,39 @@ test_that("parameter validation works", { # ---- run model and check of it works ---- -# summary(model_coa_time_vary) +# summary(time_vary_gaussian) -# bayesplot::ppc_dens_overlay(y = as.vector(Y), yrep = model_coa_time_vary$generated_quantities) -# rstan::traceplot(model_coa_time_vary$model, pars = c( +# bayesplot::ppc_dens_overlay(y = as.vector(Y), yrep = time_vary_gaussian$generated_quantities) +# rstan::traceplot(time_vary_gaussian$model, pars = c( # # "alpha0", # "alpha1", # "sigma", "lp__")) -# model_coa_time_vary$coas -test_that("test COA_TimeVarying model results to make sure its consisitent", { - mean_p0 <- model_coa_time_vary$summary[1] +# time_vary_gaussian$coas +test_that("test COA_TimeVarying model results to make sure its consistent", { + mean_p0 <- time_vary_gaussian$summary[1] expected_mean_p0 <- 0.4928 expect_equal(mean_p0, expected_mean_p0, tolerance = 0.07) }) -test_that("check to see if model_coa_time_vary classes", { - expect_type(model_coa_time_vary, "list") - expect_s4_class(model_coa_time_vary$model, "stanfit") - expect_s3_class(model_coa_time_vary$coas, "data.frame") - expect_s3_class(model_coa_time_vary$all_estimates, "data.frame") - expect_type(model_coa_time_vary$summary, "double") - expect_true(is.matrix(model_coa_time_vary$summary)) - expect_true(is.matrix(model_coa_time_vary$generated_quantities$yrep)) - expect_type(model_coa_time_vary$generated_quantities, "list") - expect_true(is.numeric(model_coa_time_vary$time)) +test_that("check time_vary_gaussian classes", { + expect_type(time_vary_gaussian, "list") + expect_s4_class(time_vary_gaussian$model, "stanfit") + expect_s3_class(time_vary_gaussian$coas, "data.frame") + expect_s3_class(time_vary_gaussian$all_estimates, "data.frame") + expect_type(time_vary_gaussian$summary, "double") + expect_true(is.matrix(time_vary_gaussian$summary)) + expect_true(is.matrix(time_vary_gaussian$generated_quantities$y_rep)) + expect_type(time_vary_gaussian$generated_quantities, "list") + expect_true(is.numeric(time_vary_gaussian$time)) }) test_that("check to see if coa returns proper info", { - expect_true("coas" %in% names(model_coa_time_vary)) - expect_equal(nrow(model_coa_time_vary$coas), model_param_ex$tsteps) + expect_true("coas" %in% names(time_vary_gaussian)) + expect_equal(nrow(time_vary_gaussian$coas), model_param_ex$tsteps) expect_equal( - colnames(model_coa_time_vary$coas), + colnames(time_vary_gaussian$coas), c( "time", "x", @@ -148,14 +148,14 @@ test_that("check to see if coa returns proper info", { ) ) - for (col in colnames(model_coa_time_vary$coas)) { - expect_type(model_coa_time_vary$coas[[col]], "double") - expect_true(all(is.finite(model_coa_time_vary$coas[[col]]))) + for (col in colnames(time_vary_gaussian$coas)) { + expect_type(time_vary_gaussian$coas[[col]], "double") + expect_true(all(is.finite(time_vary_gaussian$coas[[col]]))) } }) test_that("check to see model converged and has a good rhat", { - rhat <- model_coa_time_vary$summary[, "Rhat"] + rhat <- time_vary_gaussian$summary[, "Rhat"] expect_true(all(rhat > 0.95 & rhat < 1.05)) }) @@ -163,5 +163,61 @@ test_that("check to see model converged and has a good rhat", { test_that("check to see if gq is the correct length", { expected <- 11 - expect_true(nrow(model_coa_time_vary$generated_quantities$yrep) %in% expected) + expect_true(nrow(time_vary_gaussian$generated_quantities$y_rep) %in% expected) +}) + + +#### LOGISTIC #### +test_that("test COA_TimeVarying model results to make sure its consistent", { + mean_p0 <- time_vary_logistic$summary[1] + expected_mean_p0 <- 0.49815 + expect_equal(mean_p0, expected_mean_p0, tolerance = 0.07) +}) + + +test_that("check time_vary_logistic classes", { + expect_type(time_vary_logistic, "list") + expect_s4_class(time_vary_logistic$model, "stanfit") + expect_s3_class(time_vary_logistic$coas, "data.frame") + expect_s3_class(time_vary_logistic$all_estimates, "data.frame") + expect_type(time_vary_logistic$summary, "double") + expect_true(is.matrix(time_vary_logistic$summary)) + expect_true(is.matrix(time_vary_logistic$generated_quantities$y_rep)) + expect_type(time_vary_logistic$generated_quantities, "list") + expect_true(is.numeric(time_vary_logistic$time)) +}) + + +test_that("check to see if coa returns proper info", { + expect_true("coas" %in% names(time_vary_logistic)) + expect_equal(nrow(time_vary_logistic$coas), model_param_ex$tsteps) + expect_equal( + colnames(time_vary_logistic$coas), + c( + "time", + "x", + "y", + "x_lower", + "x_upper", + "y_lower", + "y_upper" + ) + ) + + for (col in colnames(time_vary_logistic$coas)) { + expect_type(time_vary_logistic$coas[[col]], "double") + expect_true(all(is.finite(time_vary_logistic$coas[[col]]))) + } +}) + +test_that("check to see model converged and has a good rhat", { + rhat <- time_vary_logistic$summary[, "Rhat"] + expect_true(all(rhat > 0.95 & rhat < 1.05)) +}) + +# ----- check if gq retruns the correct length ------ + +test_that("check to see if gq is the correct length", { + expected <- 11 + expect_true(nrow(time_vary_logistic$generated_quantities$yrep) %in% expected) }) diff --git a/tests/testthat/test-generated_quantities.R b/tests/testthat/test-generated_quantities.R index e79ed23..d401095 100644 --- a/tests/testthat/test-generated_quantities.R +++ b/tests/testthat/test-generated_quantities.R @@ -1,17 +1,19 @@ # ----- make all models into a list ----- all_models <- list( - model_coa_standard, - model_coa_time_vary, - model_coa_tag_int + standard_gaussian, + standard_logistic, + time_vary_gaussian, + time_vary_logistic, + tag_int_gaussian, + tag_int_logistic ) # ---- do the same for the data ----- -all_data <- list( - standata, - standata, - standata_1 +all_data <- rep( + list(standata, standata_test_tag), + times = c(4, 2) ) # set the number of draws to test -ndraws_test <- 5 +n_draws_test <- 5 # ----- create function to loop over for errors ----- @@ -20,9 +22,9 @@ call_generated_quantities <- function(overrides) { } # ----- gq_arguments to check ------ gq_args <- list( - model = model_coa_standard$model, + model = standard_gaussian$model, standata = standata, - ndraws = ndraws_test + n_draws = n_draws_test ) # ----- check if params error properly ------ @@ -33,9 +35,9 @@ params_table <- list( regex = "`model` must be a Stan object \\(from rstan or cmdstanr\\)\\." ), list( - param = "ndraws", + param = "n_draws", bad = list("a", NA, c(1, 2)), - regex = "`ndraws` must be a numeric vector that has a length of 1." + regex = "`n_draws` must be a numeric vector that has a length of 1." ) ) @@ -70,27 +72,27 @@ test_that("parameter validation works", { # create empty list to dump all gc to check -yreps <- list() +y_reps <- list() # ----- loop over generated quantities ----- for (i in seq_along(all_models)) { # Call your function - yreps[[i]] <- generated_quantities( + y_reps[[i]] <- generated_quantities( model = all_models[[i]]$model, standata = all_data[[i]], - ndraws = ndraws_test + n_draws = n_draws_test ) } -# length of each object returned bs = basic_structure -bs_returned <- c(1, 1, 2) +# length of each object returned bs = basic_structure +bs_returned <- rep(c(1, 2), times = c(4, 2)) # length of yrep give the 5 draws -length_yrep <- ndraws_test +length_y_rep <- n_draws_test # ----- checks the structure of the structure of generated_quantites ------- test_that("generated_quantities returns correct structure", { - for (s in seq_along(yreps)) { - bs <- yreps[[s]] + for (s in seq_along(y_reps)) { + bs <- y_reps[[s]] expect_type(bs, "list") expect_length(bs, bs_returned[s]) @@ -99,7 +101,7 @@ test_that("generated_quantities returns correct structure", { post_draws <- bs[[n]] expect_type(post_draws, "list") - expect_length(post_draws, length_yrep) + expect_length(post_draws, length_y_rep) for (h in seq_along(post_draws)) { one_draw <- post_draws[[h]] @@ -120,8 +122,8 @@ test_that("generated_quantities returns correct structure", { # do not test actual values as these will change test_that("generated_quantities returns correct integer ", { - for (s in seq_along(yreps)) { - bs <- yreps[[s]] + for (s in seq_along(y_reps)) { + bs <- y_reps[[s]] for (n in seq_along(bs)) { post_draws <- bs[[n]] diff --git a/tests/testthat/test-transform_gq.R b/tests/testthat/test-transform_gq.R index 483e906..a2e33d6 100644 --- a/tests/testthat/test-transform_gq.R +++ b/tests/testthat/test-transform_gq.R @@ -1,39 +1,41 @@ # ----- create example gq to transform ------- # ----- make all models into a list ----- all_models <- list( - model_coa_standard, - model_coa_time_vary, - model_coa_tag_int + standard_gaussian, + standard_logistic, + time_vary_gaussian, + time_vary_logistic, + tag_int_gaussian, + tag_int_logistic ) # ---- do the same for the data ----- -all_data <- list( - standata, - standata, - standata_1 +all_data <- rep( + list(standata, standata_testtag), + times = c(4, 2) ) # set the number of draws to test -ndraws_test <- 5 -yreps <- list() +n_draws_test <- 5 +y_reps <- list() # ----- loop over generated quantities ----- for (i in seq_along(all_models)) { # Call your function - yreps[[i]] <- generated_quantities( + y_reps[[i]] <- generated_quantities( model = all_models[[i]]$model, standata = all_data[[i]], - ndraws = ndraws_test + n_draws = n_draws_test ) } # test if the returned object matches the correct format -bs_returned <- c(1, 1, 2) -bs_names <- c(rep("yrep", 3), "testrep") +bs_returned <- rep(c(1, 2), times = c(4, 2)) +bs_names <- c(rep("y_rep", 6), "test_rep") test_that("check transformation of gq to matrix", { - for (i in seq_along(yreps)) { - tran_gq <- transform_gq(yreps[[i]]) + for (i in seq_along(y_reps)) { + tran_gq <- transform_gq(y_reps[[i]]) expect_type(tran_gq, "list") expect_length(tran_gq, bs_returned[i]) @@ -49,12 +51,12 @@ test_that("check transformation of gq to matrix", { }) test_that("check row and column names of gq in matrix", { - for (i in seq_along(yreps)) { - tran_gq <- transform_gq(yreps[[i]]) + for (i in seq_along(y_reps)) { + tran_gq <- transform_gq(y_reps[[i]]) for (n in seq_along(tran_gq)) { post_draws <- tran_gq[[n]] # check row names - expect_true(all(grepl("^(yrep|testrep)_[0-9]+$", rownames(post_draws)))) + expect_true(all(grepl("^(y_rep|test_rep)_[0-9]+$", rownames(post_draws)))) # check column names expect_true(all(grepl( @@ -67,13 +69,3 @@ test_that("check row and column names of gq in matrix", { } } }) - -# for (i in 1:n_draws) { -# y_rep_mat[i, ] <- as.vector(draws$yrep[i, , , ]) -# } -# # make sure there's no NA and make sure obs vfallls within a range -# for (i in 1:n_draws) { -# expect_false(unique(is.na( y_rep_mat[i, ]))) -# expect_true(all(y_rep_mat[i, ] >= 0 & y_rep_mat[i, ] <= 25)) -# } -# }