R/qn0.R

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

#' Mixing probability for creating new cluster
#'
#' @param lambda A scalar defining the parameter for the Truncate Poisson distribution
#'    that controls the number of change points (or its initial values)
#' @param kstar A scalar with the number maximum of change points in all clusters
#' @param Yn A vector or matrix with data sequences for a cluster
#' @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 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 N A scalar representing the number of data sequences
#' @param alpha0 A scalar defining the parameter for the Dirichlet process prior
#'    that controls the number of clusters (or its initial values)
#'
#' @returns A numerical value representing the mixing value term used to compute the probability that the given data sequence should be a singleton cluster
#'
#' @note
#' This function is called within [gibbs_alg()]. It should not be called alone.
#'
#' @seealso [gibbs_alg()]
#'
#' @examples
#' qn0(alpha0 = 1/100, w = 10, N = 5, M = 50, bs = 1000, as = 2, kstar = 2, lambda = 2, Yn = data[,1])
#'
#' @export
#'
qn0 <- function(alpha0, w, N, M, bs, as, kstar, lambda, Yn){

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


    res2 <-c(res2, sum(apply(mk, 1,
                             FUN = function(x){qn0_mk(w=w, m0=m0, bs=bs, as=as, M=M, km=km,
                                                      lambda=lambda, mk=x, Yn=Yn, kstar=kstar)})))
  }
  return(alpha0*sum(res2))
}