R/posteriorZ.R

In bpgmm: Bayesian Model Selection Approach for Parsimonious Gaussian Mixture Models

#' updatePostZ
#'
#' @param X X
#' @param m m
#' @param n n
#' @param thetaYList thetaYList
#'
#' @examples
#' set.seed(100)
#' n <- 10
#' p <- 2
#' q <- 1
#' K <- 2
#' m <- 1
#' muBar <- c(0, 0)
#' qVec <- c(1, 1)
#' constraint <- c(0, 0, 0)
#' X <- t(
#'   fabMix::simData(
#'     sameLambda = TRUE,
#'     sameSigma = TRUE,
#'     K.true = K,
#'     n = n,
#'     q = q,
#'     p = p,
#'     sINV_values = 1 / ((1:p))
#'   )$data
#' )
#' hparam <- new(
#'   "Hparam",
#'   alpha1 = 0.567755037123148,
#'   alpha2 = 1.1870201935945,
#'   delta = 2,
#'   ggamma = 2,
#'   bbeta = 3.39466184520673
#' )
#' ZOneDim <- sample(seq_len(m), n, replace = TRUE)
#' thetaYList <-
#'   new(
#'     "ThetaYList",
#'     tao = 0.366618687752634,
#'     psy = list(structure(
#'       c(
#'         4.18375613018654,
#'         0, 0, 5.46215996830771
#'       ),
#'       .Dim = c(2L, 2L)
#'     )),
#'     M = list(structure(
#'       c(
#'         3.27412045866392,
#'         -2.40544145363349
#'       ),
#'       .Dim = 1:2
#'     )),
#'     lambda = list(structure(
#'       c(
#'         2.51015961514781,
#'         -0.0741189919182549
#'       ),
#'       .Dim = 2:1
#'     )),
#'     Y = list(structure(
#'       c(
#'         -0.244239011725104,
#'         -0.26876172736886,
#'         0.193431511203083,
#'         0.41624466812811,
#'         -0.54581548068437,
#'         -0.0479517628308146,
#'         -0.633383997203325,
#'         0.856855296613208,
#'         0.792850576988512,
#'         0.268208848994559
#'       ),
#'       .Dim = c(1L, 10L)
#'     ))
#'   )
#' \donttest{
#' updatePostZ(X, m, n, thetaYList)
#' }
updatePostZ <- function(X, m, n, thetaYList) {
  tao <- thetaYList@tao
  psy <- thetaYList@psy
  M <- thetaYList@M
  lambda <- thetaYList@lambda


  pMat <- matrix(NA, m, n)
  ## evaluate density
  dMat <- matrix(NA, m, n)

  for (k in 1:m) {
    for (i in 1:n) {
      dMat[k, i] <- mvtnorm::dmvnorm(X[, i],
        mean = M[[k]], sigma = psy[[k]] + lambda[[k]] %*% t(lambda[[k]]),
        log = T
      )
    }
  }

  for (k in 1:m) {
    dMat[k, ] <- dMat[k, ] + log(tao[k])
  }

  for (i in 1:n) {
    for (k in 1:m) {
      pMat[k, i] <- calculateRatio(dMat[k, i], dMat[, i])
    }
  }

  ZOneDim <- c()
  for (i in 1:n) {
    tempProb <- as.numeric(pMat[, i])
    ZOneDim[i] <- sample(x = 1:m, size = 1, prob = tempProb)
  }
  ZOneDim
}