# R/utils.R In bpgmm: Bayesian Model Selection Approach for Parsimonious Gaussian Mixture Models

#### Documented in calculateRatiocalculateVarListchangeConstraintFormatclearCurrentThetaYlistcombineClusterParagetIndThetaYgetmodegetRemovedIndThetaYgetZmatlikelihoodlistToStrVecsumerizeZ

```#' Tool for vector to matrix
#'
#' @param ZOneDim a vector.
#' @param m the number of cluster.
#' @param n sample size.
#' @examples
#' m <- 20
#' n <- 500
#' ZOneDim <- sample(seq_len(m), n, replace = TRUE)
#' #'
#' \donttest{
#' getZmat(ZOneDim, m, n)
#' }
getZmat <- function(ZOneDim, m, n) {
Zmat <- matrix(NA, m, n)
for (i in 1:m) {
Zmat[i, ] <- as.numeric(ZOneDim == i)
}
Zmat
}

#' Log scale ratio calculation
#'
#' @param deno denominator.
#' @param nume numerator.
#' @return result of ratio
#'
#' @examples
#' deno <- log(1)
#' nume <- log(2)
#' #'
#' \donttest{
#' calculateRatio(deno, nume)
#' }
calculateRatio <- function(deno, nume) {
## deno nume both in log scale
maxNume <- max(nume)
transDeno <- deno - maxNume
transNume <- nume - maxNume
res <- exp(transDeno) / (sum(exp(transNume)))
res
}

#' Convert list of string to vector of string
#'
#' @param stringList list of string
#' @return vector of string
#'
#' @examples
#' stringList <- list("abc")
#' #'
#' \donttest{
#' listToStrVec(stringList)
#' }
listToStrVec <- function(stringList) {
for (i in 1:length(stringList)) {
stringList[[i]] <-
paste0("(", paste0(stringList[[i]], collapse = ""), ")")
}
unlist(stringList)
}

#' likelihood
#'
#' @param thetaYList thetaYList
#' @param ZOneDim ZOneDim
#' @param qqVec qqVec
#' @param muBar muBar
#' @param X X
#'
#' @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{
#' likelihood(thetaYList, ZOneDim, qVec, muBar, X)
#' }
likelihood <- function(thetaYList, ZOneDim, qqVec, muBar, X) {
m <- length(qqVec)
n <- ncol(X)
xvec <- c()
## 2.1: X
xVal <- 0
for (i in 1:n) {
## for subject i
k <- ZOneDim[i]

# if(is.vector(X[, ZOneDim == k])){
#   meanx = X[, ZOneDim == k]
# }else{
#   meanx = apply(X[, ZOneDim == k], 1, mean)
# }

meanx <- thetaYList@lambda[[k]] %*% thetaYList@Y[[k]][, i] + c(thetaYList@M[[k]])
varx <- thetaYList@psy[[k]]
xvec[i] <- (mvtnorm::dmvnorm(x = X[, i], mean = meanx, sigma = varx, log = T))
xVal <- xVal + mvtnorm::dmvnorm(x = X[, i], mean = meanx, sigma = varx, log = T)
}

xVal
# xvec
}

#' sumerizeZ
#'
#' @param Zlist  Zlist
#' @param index  index
#' @examples
#' Zlist <- list(c(1, 2, 3), c(3, 2, 1), c(2, 2, 2))
#' #'
#' \donttest{
#' sumerizeZ(Zlist)
#' }
sumerizeZ <- function(Zlist, index = 1:length(Zlist)) {
sampleSize <- length(Zlist[[1]])
res <- c()
for (i in 1:sampleSize) {
temp <- c()
for (j in index) {
temp[j] <- Zlist[[j]][i]
}
res[i] <- getmode(temp)
}
return(res)
}

#' getmode
#'
#' @param v  v
#' @examples
#' #'
#' \donttest{
#' getmode(c(1, 1, 2, 3))
#' }
getmode <- function(v) {
v <- v[!is.na(v)]
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}

#' calculateVarList
#'
#' @param psyList psyList
#' @param lambdaList lambdaList
#' @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{
#' calculateVarList(thetaYList@psy, thetaYList@lambda)
#' }
calculateVarList <- function(psyList, lambdaList) {
m <- length(psyList)
varList <- list()
for (i in 1:m) {
tempCov <- psyList[[i]] + lambdaList[[i]] %*% t(lambdaList[[i]])
varList[[i]] <- tempCov
}
return(varList)
}

#' clearCurrentThetaYlist
#'
#' @param thetaYList thetaYList
#' @param clusInd clusInd
#' @param mMax mMax
#' @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)
#'     ))
#'   )
#' clusInd <- rep(1, m)
#' mMax <- 1
#' #'
#' \donttest{
#' clearCurrentThetaYlist(thetaYList, clusInd, mMax)
#' }
clearCurrentThetaYlist <- function(thetaYList, clusInd, mMax) {
resThetaYList <- thetaYList

for (i in 1:mMax) {
if (clusInd[i] == 1) {
resThetaYList@tao[i] <- thetaYList@tao[i]
resThetaYList@psy[[i]] <- thetaYList@psy[[i]]
resThetaYList@M[[i]] <- thetaYList@M[[i]]
resThetaYList@lambda[[i]] <- thetaYList@lambda[[i]]
resThetaYList@Y[[i]] <- thetaYList@Y[[i]]
} else {
resThetaYList@tao[i] <- 0
resThetaYList@psy[[i]] <- NA
resThetaYList@M[[i]] <- NA
resThetaYList@lambda[[i]] <- NA
resThetaYList@Y[[i]] <- NA
}
}
return(resThetaYList)
}

#' combineClusterPara
#'
#' @param oldList oldList
#' @param newList newList
#' @param ind ind
#'
#' @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)
#' newList <- oldList <-
#'   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{
#' combineClusterPara(oldList, newList, 1)
#' }
combineClusterPara <- function(oldList, newList, ind) {
resList <- oldList
resList@tao <- oldList@tao * (1 - newList@tao)
resList@tao[ind] <- newList@tao

resList@psy[ind] <- newList@psy
resList@M[ind] <- newList@M
resList@lambda[ind] <- newList@lambda
resList@Y[ind] <- newList@Y
return(resList)
}

#' getIndThetaY
#'
#' @param thetaYList  thetaYList
#' @param Ind Ind
#' @examples
#' set.seed(100)
#' n <- 10
#' p <- 2
#' q <- 1
#' K <- 2
#' m <- 2
#' 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{
#' getIndThetaY(thetaYList, 1)
#' }
getIndThetaY <- function(thetaYList, Ind) {
new("ThetaYList",
tao = thetaYList@tao[Ind],
psy = thetaYList@psy[Ind],
M = thetaYList@M[Ind],
lambda = thetaYList@lambda[Ind],
Y = thetaYList@Y[Ind]
)
}

#' getRemovedIndThetaY
#'
#' @param thetaYList thetaYList
#' @param Ind Ind
#'
#' @examples
#' set.seed(100)
#' n <- 10
#' p <- 2
#' q <- 1
#' K <- 2
#' m <- 2
#' 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 = c(0.90162050961987, 0.0983794903801295),
#'   psy = list(structure(c(3.68472841602225, 0, 0, 8.34691978354054),
#'   .Dim = c(2L, 2L)), structure(c(0.785011896130842, 0, 0, 1.19022383323437),
#'   .Dim = c(2L, 2L))), M = list(structure(c(
#'     2.96424305287004,
#'     1.08454861414306
#'   ), .Dim = 1:2), structure(c(
#'     -0.232625450433964,
#'     0.984505960868685
#'   ), .Dim = 1:2)), lambda = list(structure(c(
#'     -0.964026624054337,
#'     0.89378616732449
#'   ), .Dim = 2:1), structure(c(
#'     0.533334148228635,
#'     -1.80033696090263
#'   ), .Dim = 2:1)), Y = list(structure(c(
#'     -0.15346475266988,
#'     1.6584112693271, 0.409294936277862, -1.46628591247549, -0.532753243163142,
#'     -0.332143130316749, 0.307558110800446, -0.525374243612587, 0.527667526535661,
#'     0.748193650431916
#'   ), .Dim = c(1L, 10L)), structure(c(
#'     0.571325118638535,
#'     0.542462985882966, 0.559971315637159, -1.73905343105432, -0.583549598471542,
#'     1.71264245945391, -0.327119395945831, 1.02464651767821, -1.11462280255215,
#'     0.81095592501554
#'   ), .Dim = c(1L, 10L))))
#' Ind <- 1
#' #'
#' \donttest{
#' getRemovedIndThetaY(thetaYList, Ind)
#' }
getRemovedIndThetaY <- function(thetaYList, Ind) {
res <- thetaYList

res@tao[Ind] <- NA
res@tao <- res@tao / sum(res@tao, na.rm = T)
res@psy[[Ind]] <- NA
res@M[[Ind]] <- NA
res@lambda[[Ind]] <- NA
res@Y[[Ind]] <- NA

res
}

#' changeConstraintFormat
#'
#' @param strNum strNum
#' @examples
#' #'
#' \donttest{
#' changeConstraintFormat(c(0, 0, 0))
#' }
changeConstraintFormat <- function(strNum) {
res <- ""
res <- gsub(pattern = "0", replacement = "U", x = strNum)
res <- gsub(pattern = "1", replacement = "C", x = res)
res <- gsub(pattern = "\\(", replacement = "", x = res)
res <- gsub(pattern = "\\)", replacement = "", x = res)
return(res)
}
```

## Try the bpgmm package in your browser

Any scripts or data that you put into this service are public.

bpgmm documentation built on July 2, 2020, 2:51 a.m.