R/run_gibbs.R

In BayesCPclust: A Bayesian Approach for Clustering Constant-Wise Change-Point Data

#' Runs the Gibbs sampler algorithm using using initial values for the parameters
#'
#' @param par.values A list with lists with parameters for each cluster.
#'    The first argument in each list is the number of change points,
#'    then the positions for the change points,
#'    where T_1 = 1, T_last = M + 1, and for each interval between change points
#'    you need to specify a value for the constant level. If running the
#'    Gibbs sampler for a dataset with unknown number of change points,
#'    we suggest setting the number of change points for each cluster to be zero.
#'    Check example in README file.
#' @param d A scalar representing the number of clusters.
#' @param M A scalar representing the number of points available for each observation
#' @param w A scalar representing the minimum number of points in each interval between two change points
#' @param N A scalar representing the number of observations
#' @param as The hyperparameter value for the shape parameter in the inverse-gamma prior for the variance
#'    component
#' @param bs The hyperparameter value for the scale parameter in the inverse-gamma prior for the variance
#'    component
#' @param al The hyperparameter value for the shape parameter in the gamma prior for lambda
#' @param bl The hyperparameter value for the scale parameter in the gamma prior for lambda
#' @param a The hyperparameter value for the shape parameter in the gamma prior for alpha0
#' @param b The hyperparameter value for the scale parameter in the gamma prior for alpha0
#' @param alpha0 A scalar defining the parameter for the Dirichlet process prior
#'    that controls the number of clusters (or its initial values)
#' @param lambda A scalar defining the parameter for the Truncate Poisson distribution
#'    that controls the number of change points (or its initial values)
#' @param maxIter A scalar for the number of iteration to run in the Gibbs sampler
#' @param data a matrix of size M x N with data sequences in the columns
#' @param cluster a vector with cluster assignments for each data sequence
#' @param sigma2 a vector with variance components for each data sequence
#'
#' @returns A list with estimates for each iteration of the Gibbs sampler for each parameter
#'
#' @examples
#' \donttest{
#' d = 2 # two clusters
#' N = 5 # 5 data sequences
#' M = 50 # 50 observations for each data sequence
#' maxIter = 10 # number of Gibbs sampler iterations
#'
#' data(data)
#' # initial values for each paramter and each cluster
#' par.values <- list(K = c(0, 0), Tl = list(50, 50), alpha = list(5, 10))
#' #cluster assignment for each data sequence
#' cluster <- kmeans(t(data), 2)$cluster
#' # variance for each data sequence
#' sigma2 <- apply(data, 2, var)
#' res <- run_gibbs(M, N, w = 10, d, as = 2, bs = 100, al = 2, bl = 1000, a = 2,
#'  b = 1000, alpha0 = 1/100, lambda = 2, maxIter = 10, par.values, data,
#'  cluster, sigma2)
#'}
#' @export
#'
run_gibbs <- function(M, N, w, d, as = 2, bs = 100, al = 2, bl = 1000, a = 2, b = 1000, alpha0 = 1/100, lambda = 2, maxIter = 10000, par.values, data, cluster, sigma2){

  # maximum number of change-points
  kstar <- ((M - 1)/w) - 1

  res <- gibbs_alg(alpha0 = alpha0, N = N, w = w, M = M,
                   K = par.values$K,
                   Tl = par.values$Tl,
                   cluster = cluster,
                   alpha = par.values$alpha,
                   sigma2 = sigma2,
                   bs = bs, as = as, al = al, bl = bl, a = a, b = b,
                   kstar = kstar,lambda = lambda,
                   Y = data, d = d, maxIter = maxIter)


  return(res)
}