R/RcppExports.R

In BayesChange: Bayesian Methods for Change Points Analysis

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Simulate epidemiological data
#'
#' @param S0 number of individuals in the population.
#' @param I0 number of infected individuals at time 0.
#' @param max_time maximum observed time.
#' @param beta_vec vector with the infection rate for each discrete time.
#' @param xi_0 the recovery rate of the population, must be in \eqn{(0,1)}.
#' @param user_seed seed for random distribution generation.
#' @return Function \code{sim_epi_data} returns a vector with the simulated infection times.
#'
#' @examples
#'
#' betas <- c(rep(0.45, 25),rep(0.14,25))
#'
#' inf_times <- as.numeric()
#'
#' inf_times <- sim_epi_data(10000, 10, 50, betas, 1/8)
#'
#' @export
sim_epi_data <- function(S0, I0, max_time, beta_vec, xi_0, user_seed = 1234L) {
    .Call(`_BayesChange_sim_epi_data`, S0, I0, max_time, beta_vec, xi_0, user_seed)
}

#' @name detect_cp_uni
#' @export detect_cp_uni
#'
#' @title Detect Change Points on an univariate time series.
#' @description Detect Change Points on an univariate time series.
#'
#' @param data vector of observations.
#' @param n_iterations number of MCMC iteration.
#' @param q probability of performing a split at each iterations.
#' @param a,b,c parameters of the Normal-Gamma prior for \eqn{\mu} and \eqn{\lambda}.
#' @param prior_var_phi parameters for the correlation coefficient in the likelihood.
#' @param par_theta_c,par_theta_d parameters of the shifted Gamma prior for \eqn{\theta}.
#' @param print_progress If TRUE (default) print the progress bar.
#' @param user_seed seed for random distribution generation.
#' @return Function \code{detect_cp_uni} returns a list containing the following components: \itemize{
#' \item{\code{$orders}} a matrix where each row corresponds to the output order of the corresponding iteration.
#' \item{\code{time}} computational time in seconds.
#' \item{\code{$sigma_MCMC}} traceplot for \eqn{\sigma}.
#' \item{\code{$sigma_MCMC_01}} a \eqn{0/1} vector, the \eqn{n}-th element is equal to \eqn{1} if the proposed \eqn{\sigma} was accepted, \eqn{0} otherwise.
#' \item{\code{$theta_MCMC}} traceplot for \eqn{\theta}.
#' }
#'
#' @examples
#'
#' data_vec <- as.numeric(c(rnorm(50,0,0.1), rnorm(50,1,0.25)))
#'
#' out <- detect_cp_uni(data = data_vec,
#'                             n_iterations = 2500,
#'                             q = 0.25)
#'
#'
detect_cp_uni <- function(data, n_iterations, q, a = 1, b = 1, c = 0.1, prior_var_phi = 0.1, par_theta_c = 1, par_theta_d = 1, print_progress = TRUE, user_seed = 1234L) {
    .Call(`_BayesChange_detect_cp_uni`, data, n_iterations, q, a, b, c, prior_var_phi, par_theta_c, par_theta_d, print_progress, user_seed)
}

#' @name detect_cp_multi
#' @export detect_cp_multi
#'
#' @title Detect Change Points on multivariate time series
#' @description Detect Change Points on multivariate time series
#'
#' @param data a matrix where each row is a component of the time series and the columns correpospond to the times.
#' @param n_iterations number of MCMC iterations.
#' @param q probability of performing a split at each iteration.
#' @param k_0,nu_0,S_0,m_0 parameters for the Normal-Inverse-Wishart prior for \eqn{(\mu,\lambda)}.
#' @param par_theta_c,par_theta_d parameters for the shifted Gamma prior for \eqn{\theta}.
#' @param prior_var_phi parameters for the correlation coefficient in the likelihood.
#' @param print_progress If TRUE (default) print the progress bar.
#' @param user_seed seed for random distribution generation.
#' @return Function \code{detect_cp_multi} returns a list containing the following components: \itemize{
#' \item{\code{$orders}} a matrix where each row corresponds to the output order of the corresponding iteration.
#' \item{\code{time}} computational time in seconds.
#' \item{\code{$phi_MCMC}} traceplot for \eqn{\gamma}.
#' \item{\code{$phi_MCMC_01}} a \eqn{0/1} vector, the \eqn{n}-th element is equal to \eqn{1} if the proposed \eqn{\phi} was accepted, \eqn{0} otherwise.
#' \item{\code{$sigma_MCMC}} traceplot for \eqn{\sigma}.
#' \item{\code{$sigma_MCMC_01}} a \eqn{0/1} vector, the \eqn{n}-th element is equal to \eqn{1} if the proposed \eqn{\sigma} was accepted, \eqn{0} otherwise.
#' \item{\code{$theta_MCMC}} traceplot for \eqn{\theta}.
#' }
#'
#' @examples
#'
#' data_mat <- matrix(NA, nrow = 3, ncol = 100)
#'
#' data_mat[1,] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#' data_mat[2,] <- as.numeric(c(rnorm(50,0,0.125), rnorm(50,1,0.225)))
#' data_mat[3,] <- as.numeric(c(rnorm(50,0,0.175), rnorm(50,1,0.280)))
#'
#' out <- detect_cp_multi(data = data_mat,
#'                               n_iterations = 2500,
#'                               q = 0.25,k_0 = 0.25, nu_0 = 4, S_0 = diag(1,3,3), m_0 = rep(0,3),
#'                               par_theta_c = 2, par_theta_d = 0.2, prior_var_phi = 0.1)
#'
#'
detect_cp_multi <- function(data, n_iterations, q, k_0, nu_0, S_0, m_0, par_theta_c = 1, par_theta_d = 1, prior_var_phi = 0.1, print_progress = TRUE, user_seed = 1234L) {
    .Call(`_BayesChange_detect_cp_multi`, data, n_iterations, q, k_0, nu_0, S_0, m_0, par_theta_c, par_theta_d, prior_var_phi, print_progress, user_seed)
}

#' @name detect_cp_epi
#' @export detect_cp_epi
#'
#' @title Detect Change Points on a survival function
#' @description Detect Change Points on a survival function
#'
#' @param data a matrix where each row is a component of the time series and the columns correspond to the times.
#' @param n_iterations number of MCMC iterations.
#' @param q probability of performing a split at each iteration.
#' @param M number of Monte Carlo iterations when computing the likelihood of the survival function.
#' @param xi recovery rate fixed constant for each population at each time.
#' @param a0,b0 parameters for the computation of the integrated likelihood of the survival functions.
#' @param I0_var variance for the Metropolis-Hastings estimation of the proportion of infected at time 0.
#' @param print_progress If TRUE (default) print the progress bar.
#'
#' @param user_seed seed for random distribution generation.
#' @return Function \code{detect_cp_epi} returns a list containing the following components: \itemize{
#' \item{\code{$orders}} a matrix where each row corresponds to the output order of the corresponding iteration.
#' \item{\code{time}} computational time in seconds.
#' \item{\code{$I0_MCMC}} traceplot for \eqn{I_0}.
#' \item{\code{$I0_MCMC_01}} a \eqn{0/1} vector, the \eqn{n}-th element is equal to \eqn{1} if the proposed \eqn{I_0} was accepted, \eqn{0} otherwise.
#' }
#'
#' @examples
#' \donttest{
#' data_mat <- matrix(NA, nrow = 1, ncol = 100)
#'
#' betas <- c(rep(0.45, 25),rep(0.14,75))
#'
#' inf_times <- sim_epi_data(10000, 10, 100, betas, 1/8)
#'
#' inf_times_vec <- rep(0,100)
#' names(inf_times_vec) <- as.character(1:100)
#'
#' for(j in 1:100){
#'  if(as.character(j) %in% names(table(floor(inf_times)))){
#'  inf_times_vec[j] =
#'  table(floor(inf_times))[which(names(table(floor(inf_times))) == j)]}
#' }
#'
#' data_mat[1,] <- inf_times_vec
#'
#' out <- detect_cp_epi(data = data_mat, n_iterations = 250, q = 0.5,
#'                      xi = 1/8, a0 = 40, b0 = 10, M = 250)
#'
#'}
detect_cp_epi <- function(data, n_iterations, q, M, xi, a0, b0, I0_var = 0.01, print_progress = TRUE, user_seed = 1234L) {
    .Call(`_BayesChange_detect_cp_epi`, data, n_iterations, q, M, xi, a0, b0, I0_var, print_progress, user_seed)
}

#' Clustering Epidemiological survival functions with common changes in time
#'
#' @param data a matrix where each entry is the number of infected for a population (row) at a specific discrete time (column).
#' @param n_iterations Second value
#' @param M number of Monte Carlo iterations when computing the likelihood of the survival function.
#' @param B number of orders for the normalisation constant.
#' @param L number of split-merge steps for the proposal step.
#' @param xi recovery rate fixed constant for each population at each time.
#' @param alpha_SM \eqn{\alpha} parameter for the main split-merge algorithm.
#' @param q probability of performing a split when updating the single order for the proposal procedure.
#' @param a0,b0 parameters for the computation of the integrated likelihood of the survival functions.
#' @param I0_var variance for the Metropolis-Hastings estimation of the proportion of infected at time 0.
#' @param avg_blk average number of change points for the random generated orders.
#' @param print_progress If TRUE (default) print the progress bar.
#' @param user_seed seed for random distribution generation.
#' @return Function \code{clust_cp_epi} returns a list containing the following components: \itemize{
#' \item{\code{$clust}} a matrix where each row corresponds to the output cluster of the corresponding iteration.
#' \item{\code{$orders}} a multidimensional matrix where each slice is a matrix with the orders associated to the output cluster of that iteration.
#' \item{\code{time}} computational time in seconds.
#' \item{\code{$llik}} a matrix containing the log-likelihood of each population at each iteration.
#' \item{\code{$rho}} traceplot for the proportion of infected individuals at time 0.
#' }
#'
#'@examples
#'\donttest{
#' data_mat <- matrix(NA, nrow = 5, ncol = 50)
#'
#' betas <- list(c(rep(0.45, 25),rep(0.14,25)),
#'               c(rep(0.55, 25),rep(0.11,25)),
#'               c(rep(0.50, 25),rep(0.12,25)),
#'               c(rep(0.52, 10),rep(0.15,40)),
#'               c(rep(0.53, 10),rep(0.13,40)))
#'
#'  inf_times <- list()
#'
#'  for(i in 1:5){
#'
#'    inf_times[[i]] <- sim_epi_data(10000, 10, 50, betas[[i]], 1/8)
#'
#'    vec <- rep(0,50)
#'    names(vec) <- as.character(1:50)
#'
#'    for(j in 1:50){
#'      if(as.character(j) %in% names(table(floor(inf_times[[i]])))){
#'        vec[j] = table(floor(inf_times[[i]]))[which(names(table(floor(inf_times[[i]]))) == j)]
#'      }
#'    }
#'    data_mat[i,] <- vec
#'  }
#'
#'  out <- clust_cp_epi(data = data_mat, n_iterations = 3000, M = 250, B = 1000, L = 1)
#'
#'}
#' @export
clust_cp_epi <- function(data, n_iterations, M, B, L, xi = 1/8, alpha_SM = 1, q = 0.1, a0 = 4, b0 = 10, I0_var = 0.01, avg_blk = 0.003, print_progress = TRUE, user_seed = 1234L) {
    .Call(`_BayesChange_clust_cp_epi`, data, n_iterations, M, B, L, xi, alpha_SM, q, a0, b0, I0_var, avg_blk, print_progress, user_seed)
}

#' Clustering univariate times series with common changes in time
#'
#' @param data a matrix where each row is an observation and each column corresponds to a discrete time.
#' @param n_iterations number of MCMC iterations.
#' @param B number of orders for the normalisation constant.
#' @param L number of split-merge steps for the proposal step.
#' @param phi,a,b,c parameters of the integrated likelihood.
#' @param q probability of a split in the split-merge proposal and acceleration step.
#' @param alpha_SM \eqn{\alpha} for the main split-merge algorithm.
#' @param print_progress If TRUE (default) print the progress bar.
#' @param user_seed seed for random distribution generation.
#' @return Function \code{clust_cp_uni} returns a list containing the following components: \itemize{
#' \item{\code{$clust}} a matrix where each row corresponds to the output cluster of the corresponding iteration.
#' \item{\code{$orders}} a multidimensional array where each slice is a matrix and represent an iteration. The row of each matrix correspond the order associated to the corresponding cluster.
#' \item{\code{$time}} computational time in seconds.
#' \item{\code{$norm_vec}} a vector containing the normalisation constant computed at the beginning of the algorithm.
#' }
#'
#' @examples
#'\donttest{
#' data_mat <- matrix(NA, nrow = 5, ncol = 100)
#'
#' data_mat[1,] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#' data_mat[2,] <- as.numeric(c(rnorm(50,0,0.125), rnorm(50,1,0.225)))
#' data_mat[3,] <- as.numeric(c(rnorm(50,0,0.175), rnorm(50,1,0.280)))
#' data_mat[4,] <- as.numeric(c(rnorm(25,0,0.135), rnorm(75,1,0.225)))
#' data_mat[5,] <- as.numeric(c(rnorm(25,0,0.155), rnorm(75,1,0.280)))
#'
#' out <- clust_cp_uni(data = data_mat, n_iterations = 5000, B = 1000, L = 1, phi = 0.5)
#'}
#' @export
clust_cp_uni <- function(data, n_iterations, B, L, phi, a = 1, b = 1, c = 1, q = 0.5, alpha_SM = 0.1, print_progress = TRUE, user_seed = 1234L) {
    .Call(`_BayesChange_clust_cp_uni`, data, n_iterations, B, L, phi, a, b, c, q, alpha_SM, print_progress, user_seed)
}

#' Clustering multivariate times series with common changes in time
#'
#' @param data a multidimensional matrix where each element is a matrix whose rows are the observations and columns the dimensions.
#' @param n_iterations number of MCMC iterations.
#' @param B number of orders for the normalisation constant.
#' @param L number of split-merge steps for the proposal step.
#' @param phi,k_0,nu_0,S_0,m_0 parameters of the integrated likelihood.
#' @param q probability of a split in the split-merge proposal and acceleration step.
#' @param alpha_SM \eqn{\alpha} for the main split-merge algorithm.
#' @param print_progress If TRUE (default) print the progress bar.
#' @param user_seed seed for random distribution generation.
#' @return Function \code{clust_cp_multi} returns a list containing the following components: \itemize{
#' \item{\code{$clust}} a matrix where each row corresponds to the output cluster of the corresponding iteration.
#' \item{\code{$orders}} a multidimensional array where each slice is a matrix and represent an iteration. The row of each matrix correspond the order associated to the corresponding cluster.
#' \item{\code{time}} computational time in seconds.
#' \item{\code{$norm_vec}} a vector containing the normalisation constant computed at the beginning of the algorithm.
#' }
#'
#' @examples
#'\donttest{
#' data_array <- array(data = NA, dim = c(3,100,5))
#'
#' data_array[1,,1] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#' data_array[2,,1] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#' data_array[3,,1] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#'
#' data_array[1,,2] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#' data_array[2,,2] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#' data_array[3,,2] <- as.numeric(c(rnorm(50,0,0.100), rnorm(50,1,0.250)))
#'
#' data_array[1,,3] <- as.numeric(c(rnorm(50,0,0.175), rnorm(50,1,0.280)))
#' data_array[2,,3] <- as.numeric(c(rnorm(50,0,0.175), rnorm(50,1,0.280)))
#' data_array[3,,3] <- as.numeric(c(rnorm(50,0,0.175), rnorm(50,1,0.280)))
#'
#' data_array[1,,4] <- as.numeric(c(rnorm(25,0,0.135), rnorm(75,1,0.225)))
#' data_array[2,,4] <- as.numeric(c(rnorm(25,0,0.135), rnorm(75,1,0.225)))
#' data_array[3,,4] <- as.numeric(c(rnorm(25,0,0.135), rnorm(75,1,0.225)))
#'
#' data_array[1,,5] <- as.numeric(c(rnorm(25,0,0.155), rnorm(75,1,0.280)))
#' data_array[2,,5] <- as.numeric(c(rnorm(25,0,0.155), rnorm(75,1,0.280)))
#' data_array[3,,5] <- as.numeric(c(rnorm(25,0,0.155), rnorm(75,1,0.280)))
#'
#' out <- clust_cp_multi(data = data_array, n_iterations = 3000, B = 1000, L = 1,
#'                         phi = 0.1, k_0 = 0.25, nu_0 = 5, S_0 = diag(0.1,3,3), m_0 = rep(0,3))
#'}
#' @export
clust_cp_multi <- function(data, n_iterations, B, L, phi, k_0, nu_0, S_0, m_0, q = 0.5, alpha_SM = 0.1, print_progress = TRUE, user_seed = 1234L) {
    .Call(`_BayesChange_clust_cp_multi`, data, n_iterations, B, L, phi, k_0, nu_0, S_0, m_0, q, alpha_SM, print_progress, user_seed)
}