R/postmk.R

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

#' Marginal probability of m1,m2,m3,...,mk+1
#'
#' @param clusteri A scalar with the index of a cluster
#' @param Y A matrix M x N with the data sequences
#' @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 K A vector containing the number of change points for each cluster (or its initial values)
#' @param sigma2 A vector with the variances of the data sequences (or its initial values)
#'
#' @returns A numerical vector of size k + 1 with the sampled number of observations (or bin size, mk) between each change point for a given cluster
#'
#' @note
#' This function is called within the Gibbs sampler, but it can also be called separately.
#'
#' @examples
#' data(data)
#' postmk(w = 10, M = 50, Y = data, K = c(1, 1), cluster = c(2,1,1,1,1), sigma2 = apply(data, 2, var),
#' clusteri = 1)
#'
#' @export
#'
postmk <- function(w, M, Y, K, cluster, sigma2, 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]

  # number of change points
  k <- K[clusteri]

  Cr <- sum(cluster == clusteri)

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

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

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

  xn <- sample(1:nrow(mk), 1, prob = res1/sum(res1))
  mkn <- mk[xn,]

  return(mkn)

}