R/postK.R

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

#' Marginal probability of K
#'
#' @param kstar A scalar with the number maximum of change points in all clusters
#' @param Y A matrix M x N with the data sequences
#' @param sigma2 A vector with the variances of the data sequences (or its initial values)
#' @param cluster A vector containing the cluster assignments for the data sequences (or its initial values)
#' @param M A scalar representing the number of points available for each data sequence
#' @param w A scalar representing the minimum number of points in each interval between two change points
#' @param lambda A scalar defining the parameter for the Truncate Poisson distribution
#'    that controls the number of change points (or its initial values)
#' @param clusteri A scalar with the index of a cluster
#'
#' @returns A numerical value corresponding to the sampled number of change points, k, for a given cluster
#'
#' @note
#' This function is called within the Gibbs sampler, but it can also de called separately.
#'
#' @seealso [gibbs_alg()]
#'
#' @examples
#' postK(kstar = 2, w = 10, M = 50, Y = data, cluster = c(1,1,2,1,2),
#' sigma2 = apply(data, 2, var), lambda = 2, clusteri = 1)
#'
#' @export
#'
postK <- function(kstar, w, M, Y, cluster, sigma2, lambda, clusteri){

  # filtering the cells in cluster i
  cellsn <- which(cluster == clusteri)

  # selecting the Y values for the cells in cluster i
  Yn <- as.matrix(Y[,cellsn])

  # selecting sigma2 for cells in cluster i
  sigma2n <- sigma2[cellsn]

  Cr <- sum(cluster == clusteri)

  res2 <- c()
  min_H <- numeric(kstar+1)
  min_x <- Inf
  for(k in 0:kstar){
    m0 <- M - 1 -(k+1)*w
    mk <- RcppAlgos::permuteGeneral(0:m0, k + 1,
                                    constraintFun = "sum",
                                    comparisonFun = "==", limitConstraints = m0,
                                    repetition = TRUE)
    Hlist <- apply(mk, 1,
                   FUN = function(x){postK_mk(k = k, m0 = m0, w = w, M = M, Yn = Yn,
                                              sigma2n = sigma2n, cellsn = cellsn, mk = x, Cr = Cr)})

    min_x <- logsumexp(Hlist, min_x)$min_x


  }


  for(k in 0:kstar){
    m0 <- M - 1 -(k+1)*w
    mk <- RcppAlgos::permuteGeneral(0:m0, k + 1,
                                    constraintFun = "sum",
                                    comparisonFun = "==", limitConstraints = m0,
                                    repetition = TRUE)

    Hlist <- apply(mk, 1,
                   FUN = function(mkx){postK_mk(k = k, m0 = m0, w = w, M = M, Yn = Yn,
                                                sigma2n = sigma2n, cellsn = cellsn, mk = mkx, Cr = Cr)})

    res1 <- c()
    res1 <- exp(logsumexp(Hlist, min_x)$x)*exp(apply(mk, 1, function(z){(lfactorial(m0) - sum(lfactorial(z))) - m0*log(k+1)}))

    res2 <- c(res2, sum(res1)*pk(k = k, kstar = kstar, lambda = lambda))

  }

  return(c(0:kstar)[which.max(res2)])

}