R/qn0_mk.R

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

#' Mixing probability for creating new cluster per bin
#'
#' @param km A scalar for the number of changes points in a cluster
#' @param mk A matrix with all possible values to distribute between 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 kstar A scalar with the number maximum of change points in all clusters
#' @param m0 A scalar for the number of positions available to define change-points positions
#' @param Yn A vector with a data sequence
#' @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
#'
#' @returns A numerical value representing the mixing value term used to compute the probability that the given data sequence should be a singleton cluster for a given bin size.
#' @note
#' This function is called within [qn0()]. It should not be called alone.
#'
#' @seealso [qn0()], [gibbs_alg()]
#' @examples
#' data(data)
#' M <- 50; k <- 0; w <- 10;
#' m0 <- M - 1 -(k+1)*w
#' for(k in 0:2){
#' mk <- RcppAlgos::permuteGeneral(0:m0, k + 1,
#' constraintFun = "sum",
#' comparisonFun = "==", limitConstraints = m0,
#' repetition = TRUE)}
#' out <- qn0_mk(w = 10, m0 = m0, bs = 1000, as = 2, M = 50, km = 1,
#'  lambda = 2, mk = mk[1,], Yn = data[,1], kstar = 2)
#' @export
#'
qn0_mk <- function(w, m0, bs, as, M, km, lambda, mk, Yn, kstar){

  Tl <- c()
  if(km == 0){
    Tl <- M
    Xn <- rep(1, Tl)
  } else{
    for(i in 2:(km+1)){
      Tl[1] <- 1
      Tl[i] <- mk[i-1] + Tl[i-1] + w
      #print(i)
    }
    Tl <- append(Tl, M + 1)
    Xn <- apply(diag(km + 1), 2, function(x){rep(x, diff(Tl))})

  }

  Vn <- t(Xn)%*%Xn
  B <- t(Yn)%*%Yn - t(Yn)%*%Xn%*%solve(Vn)%*%t(Xn)%*%Yn

  Hmk <- -((M-(km+1))/2)*log(2*pi) -(1/2)*log(det(Vn)) - as*log(bs) - (((M-(km+1))/2) + as)*log((B/2) + (1/bs)) -
    lgamma(as) + lgamma((((M-(km+1))/2) + as))

  res1 <- exp(Hmk)*exp((lfactorial(m0) - sum(lfactorial(mk))))*((1/(km+1))^m0)*pk(km, kstar, lambda)

  return(res1)
}