R/qnj.R

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

#' Mixing probability for getting assigned to an existing cluster
#'
#' @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 N A scalar representing the number of data sequences
#' @param K A vector containing the number of change points for each cluster (or its initial values)
#' @param Tl A list containing a vector for each cluster determining the change-point positions in each cluster
#'    (or its initial values)
#' @param cluster A vector containing the cluster assignments for the data sequences (or its initial values)
#' @param alpha A list containing a vector for each cluster determining the constant level values
#'    for each interval between change points in each cluster (or its initial values)
#'
#' @returns A vector of same size as the vector `cluster` corresponding to the mixing term value used to compute the probability that the given data sequence `Yn` should be part of each existing cluster
#'
#' @note
#' This function is called within the Gibbs sampler. It should not be called alone.
#'
#' @seealso [gibbs_alg()]
#'
#' @examples
#' qnj(N = 5, M = 50, as = 2, bs = 1000, Yn = data[,1], alpha = c(10, 10),
#'  cluster = c(1,1,2,1,2), Tl = c(50,50), K = c(0,0))
#'
#' @export
#'
qnj <- function(N, M, as, bs, Yn, alpha, cluster, Tl, K){
  res3 <- c()
  for(clusteri in 1:length(alpha)){
    alphan <- alpha[[clusteri]]
    kj <- K[clusteri]
    Tlj <- Tl[[clusteri]]
    if(kj == 0){
      Xn <- rep(1, Tlj)

      Nj <- sum(cluster == clusteri)

      Hj <- lgamma(M/2 + as) - (M/2 + as)*log(((t((Yn - Xn*alphan))%*%(Yn - Xn*alphan))/2) + (1/bs))
      -(M/2)*log(2*pi) - lgamma(as) - as*log(bs)
    } else{
      Xn <- apply(diag(kj+1), 2, function(x){rep(x, diff(Tlj))})

      Nj <- sum(cluster == clusteri)

      Hj <- lgamma(M/2 + as) - (M/2 + as)*log(((t(Yn - Xn%*%alphan)%*%(Yn - Xn%*%alphan))/2) + (1/bs))
      -(M/2)*log(2*pi) - lgamma(as) - as*log(bs)
    }

    res3 <- c(res3, Nj*exp(Hj))

  }
  return(res3)
}