R/Eigen.HMM_fit.R

In MatrixHMM: Parsimonious Families of Hidden Markov Models for Matrix-Variate Longitudinal Data

#' Fitting Parsimonious Hidden Markov Models for Matrix-Variate Longitudinal Data
#'
#' Fits parsimonious Hidden Markov Models for matrix-variate longitudinal data using ECM algorithms.
#' The models are based on the matrix-variate normal, matrix-variate t, and matrix-variate contaminated normal distributions.
#' Parallel computing is implemented and highly recommended for faster model fitting.
#'
#' @param Y An array with dimensions \code{p} x \code{r} x \code{num} x \code{t}, where \code{p} is the number of
#'     variables in the rows of each data matrix, \code{r} is the number of variables in the columns of each
#'     data matrix, \code{num} is the number of data observations, and \code{t} is the number of time points.
#' @param init.par A list of initial values for starting the algorithms, as generated by the \code{Eigen.HMM_init()} function.
#' @param tol A numeric value specifying the tolerance level for the ECM algorithms' convergence.
#' @param maxit A numeric value specifying the maximum number of iterations for the ECM algorithms.
#' @param nThreads A positive integer indicating the number of cores to use for parallel processing.
#' @param verbose A logical value indicating whether to display the running output.
#'
#' @return A list containing the following elements:
#' \item{results}{A list of the results from the fitted models.}
#' \item{c.time}{A numeric value providing information on the computational time required to fit all models for each state.}
#' \item{models}{A data frame listing the models that were fitted.}
#' @export
#' @importFrom foreach %dopar%
#' @examples
#' data(simData)
#' Y <- simData$Y
#' init <- Eigen.HMM_init(Y = Y, k = 2, density = "MVT", mod.row = "EEE", mod.col = "EE", nstartR = 10)
#' fit <- Eigen.HMM_fit(Y = Y, init.par = init, nThreads = 1)
Eigen.HMM_fit <- function(Y, init.par = NULL, tol = 0.001, maxit = 500, nThreads = 1, verbose = FALSE) {
  k <- init.par[[2]]
  list.comb <- init.par[[3]]
  pt.mod <- nrow(list.comb)
  density <- init.par[[6]]

  tol2 <- 0.001
  maxit2 <- 100

  comb <- function(x, ...) {
    lapply(
      seq_along(x),
      function(i) c(x[[i]], lapply(list(...), function(y) y[[i]]))
    )
  }
  oper <- vector(mode = "list", length = length(k))
  time <- numeric(length(k))

  Pars.HMM <- function(Y, k, density, init.par = NULL, mod.row = NULL, mod.col = NULL, tol = 0.001, tol2 = 0.001, maxit = 500, maxit2 = 100) {
    ptm <- proc.time()

    # Functions

    dMVnorm <- function(X, M, U, V) {
      tr <- function(x) {
        return(sum(diag(x)))
      }

      num <- dim(X)[3] # sample size
      p <- nrow(X) # rows of X
      r <- ncol(X) # columns of X

      if (is.na(num)) {
        X <- as.matrix(X)
        delta <- tr(solve(U) %*% (X - M) %*% solve(V) %*% t(X - M))
      } else {
        delta <- sapply(1:num, function(i) tr(solve(U) %*% (X[, , i] - M) %*% solve(V) %*% t(X[, , i] - M)))
      }

      pdf <- (2 * pi)^(-(p * r) / 2) * det(U)^(-r / 2) * det(V)^(-p / 2) * exp(-1 / 2 * delta)

      return(pdf)
    }
    dMVT <- function(X, M, U, V, nu) {
      num <- dim(X)[3] # sample size
      p <- nrow(X) # rows of X
      r <- ncol(X) # columns of X

      if (is.na(num)) {
        X <- as.matrix(X)
        delta <- tr(solve(U) %*% (X - M) %*% solve(V) %*% t(X - M))
      } else {
        delta <- sapply(1:num, function(i) tr(solve(U) %*% (X[, , i] - M) %*% solve(V) %*% t(X[, , i] - M)))
      }

      pdfvar <- (1 + delta / nu)^(-0.5 * ((p * r) + nu))
      pdfconst <- (det(U)^(-r / 2) * det(V)^(-p / 2) * gamma(0.5 * ((p * r) + nu))) / ((pi * nu)^(0.5 * (p * r)) * gamma(nu * 0.5))

      PDF <- pdfconst * pdfvar

      return(PDF)
    }
    dMVcont <- function(X, M, U, V, alpha, eta) {
      dMVnorm <- function(X, M, U, V) {
        num <- dim(X)[3] # sample size
        p <- nrow(X) # rows of X
        r <- ncol(X) # columns of X

        if (is.na(num)) {
          X <- as.matrix(X)
          delta <- tr(solve(U) %*% (X - M) %*% solve(V) %*% t(X - M))
        } else {
          delta <- sapply(1:num, function(i) tr(solve(U) %*% (X[, , i] - M) %*% solve(V) %*% t(X[, , i] - M)))
        }

        pdf <- (2 * pi)^(-(p * r) / 2) * det(U)^(-r / 2) * det(V)^(-p / 2) * exp(-1 / 2 * delta)

        return(pdf)
      }

      pdf <- alpha * dMVnorm(X = X, M = M, U = U, V = V) + (1 - alpha) * dMVnorm(X = X, M = M, U = eta * U, V = V)

      return(pdf)
    }

    tr <- function(x) {
      return(sum(diag(x)))
    }
    dec <- function(x) {
      if ((x %% 1) != 0) {
        nchar(strsplit(sub("0+$", "", as.character(x)), ".", fixed = TRUE)[[1]][[2]])
      } else {
        return(0)
      }
    }

    # Dimensions

    num0 <- dim(Y)[3]
    p <- nrow(Y)
    r <- ncol(Y)
    t <- dim(Y)[4]

    num <- num0 * t
    Yresh <- array(Y, dim = c(p, r, num))

    ## Objects related to the two covariance matrices

    WRY <- array(0, dim = c(p, p, k))
    WCY <- array(0, dim = c(r, r, k))

    phiY.k <- phiX.k <- numeric(k)
    deltaUY.k <- gammaUY.k <- array(0, dim = c(p, p, k))
    deltaVY.k <- gammaVY.k <- array(0, dim = c(r, r, k))

    temp.numdeltaR <- ftemp.r <- tempomegaR <- V_EVV.UY.K <- array(0, dim = c(p, p, k)) # for VEI, VEE, VEV, MM object
    temp.phi2 <- numeric(k) # for EVI, EVE, EVV
    tempW_EEV <- tempW_EV <- vector("list", k) # for EEV, VEV, EV
    ftemp.c <- tempomegaC <- array(0, dim = c(r, r, k)) # MM object

    ## Other objects

    post <- dens <- densN.temp <- array(0, c(num, k), dimnames = list(1:(num), paste("comp.", 1:k, sep = "")))
    post2 <- array(NA, c(k, k, num0, t - 1)) # transition probabilities
    w <- logw <- matrix(0, nrow = num, ncol = k)
    Av <- Bv <- matrix(NA, nrow = num, ncol = k)

    # Preliminary definition of convergence criteria for EM/MM algorithms

    check <- 0
    check2 <- 0
    loglik.old <- -Inf
    loglik.new <- NULL
    ll <- NULL
    mark <- 1
    MM.r.old <- -Inf
    MM.c.old <- -Inf
    m.iter <- 0
    m.iter2 <- 0

    ### Algorithm ###

    M <- init.par$M
    sigmaUY <- init.par$sigmaUY
    sigmaVY <- init.par$sigmaVY
    prior <- init.par$prior
    f <- init.par$f
    A <- init.par$A
    B <- init.par$B
    PI <- init.par$PI
    nu <- init.par$nu
    alpha <- init.par$alpha
    eta <- init.par$eta

    ei.UY <- ei.VY <- vector(mode = "list", length = k)

    for (j in 1:k) {
      ei.UY[[j]] <- eigen(sigmaUY[, , j])
      ei.VY[[j]] <- eigen(sigmaVY[, , j])

      phiY.k[j] <- prod(ei.UY[[j]][["values"]])^(1 / p)
      deltaUY.k[, , j] <- diag(ei.UY[[j]][["values"]] / phiY.k[j])
      gammaUY.k[, , j] <- ei.UY[[j]][["vectors"]]

      deltaVY.k[, , j] <- diag(ei.VY[[j]][["values"]])
      gammaVY.k[, , j] <- ei.VY[[j]][["vectors"]]
    }

    if (mod.row %in% c("EVE", "VVE")) {
      gammaU <- rowSums(gammaUY.k, dims = 2) / k
    }
    if (mod.col == "VE") {
      gammaV <- rowSums(gammaVY.k, dims = 2) / k
    }

    ### Estimation ###

    while (check < 1) {
      m.iter <- m.iter + 1

      ### E - STEP ###

      for (j in 1:k) {
        Av[, j] <- as.vector(A[, , j])
        Bv[, j] <- as.vector(B[, , j])
      }

      AvBv <- Av + Bv

      for (i in 1:num) {
        tempmax <- max(AvBv[i, ])
        tempres <- tempmax + log(sum(exp(AvBv[i, ] - tempmax)))

        post[i, ] <- exp(AvBv[i, ] - tempres)
        post[i, ] <- post[i, ] / sum(post[i, ])
      }

      if (density == "MVT") {
        for (j in 1:k) {
          delta <- sapply(1:num, function(i) tr(solve(sigmaUY[, , j]) %*% (Yresh[, , i] - M[, , j]) %*% solve(sigmaVY[, , j]) %*% t(Yresh[, , i] - M[, , j])))

          numer <- (p * r) + nu[j]
          den <- nu[j] + delta

          numer[numer < .Machine$double.xmin] <- .Machine$double.xmin
          den[den < .Machine$double.xmin] <- .Machine$double.xmin

          w[, j] <- numer / den

          numer2 <- digamma(((p * r) + nu[j]) / 2)
          den2 <- log(0.5 * (nu[j] + delta))

          numer2[numer2 < .Machine$double.xmin] <- .Machine$double.xmin
          den2[den2 < .Machine$double.xmin] <- .Machine$double.xmin

          logw[, j] <- numer2 - den2
        }
      }

      if (density == "MVCN") {
        for (j in 1:k) {
          densN.temp[, j] <- alpha[j] * dMVnorm(X = Yresh, M = M[, , j], U = sigmaUY[, , j], V = sigmaVY[, , j])

          densN.temp[, j] <- (densN.temp[, j] <= 10^(-323)) * 10^(-323) + (densN.temp[, j] > 10^(-323)) * densN.temp[, j]

          dens[, j] <- dMVcont(X = Yresh, M = M[, , j], U = sigmaUY[, , j], V = sigmaVY[, , j], alpha = alpha[j], eta = eta[j])

          dens[, j] <- (dens[, j] <= 10^(-323)) * 10^(-323) + (dens[, j] > 10^(-323)) * dens[, j]

          w[, j] <- densN.temp[, j] / dens[, j]
        }
      }

      for (n in 1:num0) {
        for (T in 1:(t - 1)) {
          amax <- max(c(A[n, T, ]))
          bmax <- max(c(B[n, T + 1, ]))

          post2[, , n, T] <- diag(exp(A[n, T, ] - amax), nrow = k, ncol = k) %*% PI %*% (diag(exp(B[n, T + 1, ] - bmax), nrow = k, ncol = k) * diag(f[n, T + 1, ], nrow = k, ncol = k))
          post2[, , n, T] <- post2[, , n, T] / sum(post2[, , n, T]) # posterior probabilities (zz)
        }
      }

      ### M - STEP ###

      prior <- colMeans(matrix(post[1:num0, ], nrow = num0, ncol = k)) # initial prob.
      post2a <- as.matrix(apply(post2, c(1, 2), sum))
      PI <- diag(1 / rowSums(post2a), nrow = k, ncol = k) %*% post2a # Transition probability matrix

      if (density == "MVN") {
        for (j in 1:k) {
          M[, , j] <- rowSums(Yresh * post[, j][slice.index(Yresh, 3)], dims = 2) / sum(post[, j])

          MX <- array(M[, , j], dim = c(p, r, num))
          WRY[, , j] <- tensor::tensor(aperm(tensor::tensor((Yresh - MX) * post[, j][slice.index(Yresh - MX, 3)], solve(sigmaVY[, , j]), 2, 1), c(1, 3, 2)), aperm((Yresh - MX), c(2, 1, 3)), c(2, 3), c(1, 3))
        }
      } else if (density == "MVT") {
        for (j in 1:k) {
          M[, , j] <- rowSums(Yresh * (w[, j] * post[, j])[slice.index(Yresh, 3)], dims = 2) / sum(w[, j] * post[, j])

          MX <- array(M[, , j], dim = c(p, r, num))
          WRY[, , j] <- tensor::tensor(aperm(tensor::tensor((Yresh - MX) * (w[, j] * post[, j])[slice.index(Yresh - MX, 3)], solve(sigmaVY[, , j]), 2, 1), c(1, 3, 2)), aperm((Yresh - MX), c(2, 1, 3)), c(2, 3), c(1, 3))
        }
      } else if (density == "MVCN") {
        for (j in 1:k) {
          M[, , j] <- rowSums(Yresh * ((w[, j] + ((1 - w[, j]) / eta[j])) * post[, j])[slice.index(Yresh, 3)], dims = 2) / sum((w[, j] + ((1 - w[, j]) / eta[j])) * post[, j])

          MX <- array(M[, , j], dim = c(p, r, num))
          WRY[, , j] <- tensor::tensor(aperm(tensor::tensor((Yresh - MX) * (post[, j] * (w[, j] + ((1 - w[, j]) / eta[j])))[slice.index(Yresh - MX, 3)], solve(sigmaVY[, , j]), 2, 1), c(1, 3, 2)), aperm((Yresh - MX), c(2, 1, 3)), c(2, 3), c(1, 3))
        }
      }

      # ROWS COVARIANCE MATRIX Y #

      if (mod.row == "EII") {
        phiY <- tr(rowSums(WRY, dims = 2)) / ((num) * p * r)

        for (j in 1:k) {
          sigmaUY[, , j] <- phiY * diag(1, p, p)
        }
      }

      if (mod.row == "VII") {
        for (j in 1:k) {
          phiY.k[j] <- tr(WRY[, , j]) / (p * r * sum(post[, j]))
          sigmaUY[, , j] <- phiY.k[j] * diag(1, p, p)
        }
      }

      if (mod.row == "EEI") {
        deltaU <- diag(diag(rowSums(WRY, dims = 2)), p, p) / (det(diag(diag(rowSums(WRY, dims = 2)), p, p)))^(1 / p)

        phiY <- (det(diag(diag(rowSums(WRY, dims = 2)), p, p)))^(1 / p) / ((num) * r)

        for (j in 1:k) {
          sigmaUY[, , j] <- phiY * deltaU
        }
      }

      if (mod.row == "VEI") {
        for (j in 1:k) {
          temp.numdeltaR[, , j] <- (1 / phiY.k[j]) * WRY[, , j]
        }

        deltaU <- diag(diag(rowSums(temp.numdeltaR, dims = 2)), p, p) / (det(diag(diag(rowSums(temp.numdeltaR, dims = 2)), p, p)))^(1 / p)

        for (j in 1:k) {
          phiY.k[j] <- (tr(solve(deltaU) %*% WRY[, , j])) / (p * r * sum(post[, j]))

          sigmaUY[, , j] <- phiY.k[j] * deltaU
        }
      }

      if (mod.row == "EVI") {
        for (j in 1:k) {
          deltaUY.k[, , j] <- diag(diag(WRY[, , j]), p, p) / (det(diag(diag(WRY[, , j]), p, p)))^(1 / p)

          temp.phi2[j] <- det(diag(diag(WRY[, , j]), p, p))^(1 / p)
        }

        phiY <- sum(temp.phi2) / ((num) * r)

        for (j in 1:k) {
          sigmaUY[, , j] <- phiY * deltaUY.k[, , j]
        }
      }

      if (mod.row == "VVI") {
        for (j in 1:k) {
          deltaUY.k[, , j] <- diag(diag(WRY[, , j]), p, p) / (det(diag(diag(WRY[, , j]), p, p)))^(1 / p)

          phiY.k[j] <- det(diag(diag(WRY[, , j]), p, p))^(1 / p) / (r * sum(post[, j]))

          sigmaUY[, , j] <- phiY.k[j] * deltaUY.k[, , j]
        }
      }

      if (mod.row == "EEE") {
        for (j in 1:k) {
          sigmaUY[, , j] <- rowSums(WRY, dims = 2) / (num * r)
        }
      }

      if (mod.row == "VEE") {
        for (j in 1:k) {
          temp.numdeltaR[, , j] <- (1 / phiY.k[j]) * WRY[, , j]
        }

        deltaU <- rowSums(temp.numdeltaR, dims = 2) / ((det(rowSums(temp.numdeltaR, dims = 2)))^(1 / p))

        for (j in 1:k) {
          phiY.k[j] <- tr(solve(deltaU) %*% WRY[, , j]) / (p * r * sum(post[, j]))

          sigmaUY[, , j] <- phiY.k[j] * deltaU
        }
      }

      if (mod.row == "EVE") {
        while (check2 < 1) {
          m.iter2 <- m.iter2 + 1

          for (j in 1:k) {
            ftemp.r[, , j] <- tcrossprod(solve(deltaUY.k[, , j]), gammaU) %*% WRY[, , j] - max(eigen(WRY[, , j])$values) * tcrossprod(solve(deltaUY.k[, , j]), gammaU)
          }

          f <- rowSums(ftemp.r, dims = 2)

          MM.r.new <- tr(f %*% gammaU)

          if ((abs(MM.r.new - MM.r.old)) < tol2 | m.iter2 == maxit2) {
            check2 <- 1
            res.svd <- svd(f)
            gammaU <- tcrossprod(res.svd$v, res.svd$u)
          } else {
            res.svd <- svd(f)
            gammaU <- tcrossprod(res.svd$v, res.svd$u)
          }

          MM.r.old <- MM.r.new
        }

        m.iter2 <- 0
        check2 <- 0
        MM.r.old <- -Inf

        for (j in 1:k) {
          deltaUY.k[, , j] <- diag(diag(crossprod(gammaU, WRY[, , j]) %*% gammaU), p, p) / (det(diag(diag(crossprod(gammaU, WRY[, , j]) %*% gammaU), p, p)))^(1 / p)

          temp.phi2[j] <- tr(gammaU %*% tcrossprod(solve(deltaUY.k[, , j]), gammaU) %*% WRY[, , j])
        }

        phiY <- sum(temp.phi2) / ((num) * p * r)

        for (j in 1:k) {
          sigmaUY[, , j] <- phiY * gammaU %*% tcrossprod(deltaUY.k[, , j], gammaU)
        }
      }

      if (mod.row == "VVE") {
        while (check2 < 1) {
          m.iter2 <- m.iter2 + 1

          for (j in 1:k) {
            ftemp.r[, , j] <- tcrossprod(solve(deltaUY.k[, , j]), gammaU) %*% WRY[, , j] - max(eigen(WRY[, , j])$values) * tcrossprod(solve(deltaUY.k[, , j]), gammaU)
          }

          f <- rowSums(ftemp.r, dims = 2)

          MM.r.new <- tr(f %*% gammaU)

          if ((abs(MM.r.new - MM.r.old)) < tol2 | m.iter2 == maxit2) {
            check2 <- 1
            res.svd <- svd(f)
            gammaU <- tcrossprod(res.svd$v, res.svd$u)
          } else {
            res.svd <- svd(f)
            gammaU <- tcrossprod(res.svd$v, res.svd$u)
          }

          MM.r.old <- MM.r.new
        }

        m.iter2 <- 0
        check2 <- 0
        MM.r.old <- -Inf

        for (j in 1:k) {
          deltaUY.k[, , j] <- diag(diag(crossprod(gammaU, WRY[, , j]) %*% gammaU), p, p) / (det(diag(diag(crossprod(gammaU, WRY[, , j]) %*% gammaU), p, p)))^(1 / p)
          phiY.k[j] <- (det(diag(diag(crossprod(gammaU, WRY[, , j]) %*% gammaU), p, p))^(1 / p)) / (r * sum(post[, j]))
          sigmaUY[, , j] <- phiY.k[j] * gammaU %*% tcrossprod(deltaUY.k[, , j], gammaU)
        }
      }

      if (mod.row == "EEV") {
        for (j in 1:k) {
          tempW_EEV[[j]] <- eigen(WRY[, , j])

          gammaUY.k[, , j] <- tempW_EEV[[j]][["vectors"]]

          tempomegaR[, , j] <- diag(tempW_EEV[[j]][["values"]], p, p)
        }

        deltaU <- rowSums(tempomegaR, dims = 2) / ((det(rowSums(tempomegaR, dims = 2)))^(1 / p))

        phiY <- ((det(rowSums(tempomegaR, dims = 2)))^(1 / p)) / ((num) * r)

        for (j in 1:k) {
          sigmaUY[, , j] <- phiY * gammaUY.k[, , j] %*% tcrossprod(deltaU, gammaUY.k[, , j])
        }
      }

      if (mod.row == "VEV") {
        for (j in 1:k) {
          tempW_EEV[[j]] <- eigen(WRY[, , j])

          gammaUY.k[, , j] <- tempW_EEV[[j]][["vectors"]]

          tempomegaR[, , j] <- diag(tempW_EEV[[j]][["values"]], p, p)

          temp.numdeltaR[, , j] <- (1 / phiY.k[j]) * tempomegaR[, , j]
        }

        deltaU <- rowSums(temp.numdeltaR, dims = 2) / ((det(rowSums(temp.numdeltaR, dims = 2)))^(1 / p))

        for (j in 1:k) {
          phiY.k[j] <- tr(tempomegaR[, , j] %*% solve(deltaU)) / (p * r * sum(post[, j]))

          sigmaUY[, , j] <- phiY.k[j] * gammaUY.k[, , j] %*% tcrossprod(deltaU, gammaUY.k[, , j])
        }
      }

      if (mod.row == "EVV") {
        for (j in 1:k) {
          V_EVV.UY.K[, , j] <- WRY[, , j] / ((det(WRY[, , j]))^(1 / p))

          temp.phi2[j] <- det(WRY[, , j])^(1 / p)
        }

        phiY <- sum(temp.phi2) / ((num) * r)

        for (j in 1:k) {
          sigmaUY[, , j] <- phiY * V_EVV.UY.K[, , j]
        }
      }

      if (mod.row == "VVV") {
        for (j in 1:k) {
          sigmaUY[, , j] <- WRY[, , j] / (r * sum(post[, j]))
        }
      }

      # COLUMNS COVARIANCE MATRIX Y #

      if (density == "MVN") {
        for (j in 1:k) {
          MX <- array(M[, , j], dim = c(p, r, num))
          WCY[, , j] <- tensor::tensor(aperm(tensor::tensor(aperm(Yresh - MX, c(2, 1, 3)) * post[, j][slice.index(aperm(Yresh - MX, c(2, 1, 3)), 3)], solve(sigmaUY[, , j]), 2, 1), c(1, 3, 2)), (Yresh - MX), c(2, 3), c(1, 3))
        }
      } else if (density == "MVT") {
        for (j in 1:k) {
          MX <- array(M[, , j], dim = c(p, r, num))
          WCY[, , j] <- tensor::tensor(aperm(tensor::tensor(aperm(Yresh - MX, c(2, 1, 3)) * (w[, j] * post[, j])[slice.index(aperm(Yresh - MX, c(2, 1, 3)), 3)], solve(sigmaUY[, , j]), 2, 1), c(1, 3, 2)), (Yresh - MX), c(2, 3), c(1, 3))
        }
      } else if (density == "MVCN") {
        for (j in 1:k) {
          MX <- array(M[, , j], dim = c(p, r, num))
          WCY[, , j] <- tensor::tensor(aperm(tensor::tensor(aperm(Yresh - MX, c(2, 1, 3)) * (post[, j] * (w[, j] + ((1 - w[, j]) / eta[j])))[slice.index(aperm(Yresh - MX, c(2, 1, 3)), 3)], solve(sigmaUY[, , j]), 2, 1), c(1, 3, 2)), (Yresh - MX), c(2, 3), c(1, 3))
        }
      }

      if (mod.col == "II") {
        for (j in 1:k) {
          sigmaVY[, , j] <- diag(1, r, r)
        }
      }

      if (mod.col == "EI") {
        deltaVY <- diag(diag(rowSums(WCY, dims = 2)), r, r) / (det(diag(diag(rowSums(WCY, dims = 2)), r, r)))^(1 / r)

        for (j in 1:k) {
          sigmaVY[, , j] <- deltaVY
        }
      }

      if (mod.col == "VI") {
        for (j in 1:k) {
          sigmaVY[, , j] <- diag(diag(WCY[, , j]), r, r) / (det(diag(diag(WCY[, , j]), r, r)))^(1 / r)
        }
      }

      if (mod.col == "EE") {
        for (j in 1:k) {
          sigmaVY[, , j] <- rowSums(WCY, dims = 2) / ((det(rowSums(WCY, dims = 2)))^(1 / r))
        }
      }

      if (mod.col == "VE") {
        while (check2 < 1) {
          m.iter2 <- m.iter2 + 1

          for (j in 1:k) {
            ftemp.c[, , j] <- tcrossprod(solve(deltaVY.k[, , j]), gammaV) %*% WCY[, , j] - max(eigen(WCY[, , j])$values) * tcrossprod(solve(deltaVY.k[, , j]), gammaV)
          }

          f.C <- rowSums(ftemp.c, dims = 2)

          MM.c.new <- tr(f.C %*% gammaV)

          if ((abs(MM.c.new - MM.c.old)) < tol2 | m.iter2 == maxit2) {
            check2 <- 1
            res.svd.C <- svd(f.C)
            gammaV <- tcrossprod(res.svd.C$v, res.svd.C$u)
          } else {
            res.svd.C <- svd(f.C)
            gammaV <- tcrossprod(res.svd.C$v, res.svd.C$u)
          }

          MM.c.old <- MM.c.new
        }

        m.iter2 <- 0
        check2 <- 0
        MM.c.old <- -Inf

        for (j in 1:k) {
          deltaVY.k[, , j] <- diag(diag(crossprod(gammaV, WCY[, , j]) %*% gammaV), r, r) / (det(diag(diag(crossprod(gammaV, WCY[, , j]) %*% gammaV), r, r)))^(1 / r)
        }

        for (j in 1:k) {
          sigmaVY[, , j] <- gammaV %*% tcrossprod(deltaVY.k[, , j], gammaV)
        }
      }

      if (mod.col == "EV") {
        for (j in 1:k) {
          tempW_EV[[j]] <- eigen(WCY[, , j])

          gammaVY.k[, , j] <- tempW_EV[[j]][["vectors"]]

          tempomegaC[, , j] <- diag(tempW_EV[[j]][["values"]], r, r)
        }

        deltaVY <- rowSums(tempomegaC, dims = 2) / ((det(rowSums(tempomegaC, dims = 2)))^(1 / r))

        for (j in 1:k) {
          sigmaVY[, , j] <- gammaVY.k[, , j] %*% tcrossprod(deltaVY, gammaVY.k[, , j])
        }
      }

      if (mod.col == "VV") {
        for (j in 1:k) {
          sigmaVY[, , j] <- WCY[, , j] / ((det(WCY[, , j]))^(1 / r))
        }
      }

      # Tail parameters #

      if (density == "MVT") {
        for (j in 1:k) {
          nu[j] <- stats::optimize(function(v) abs(log(v / 2) + 1 - digamma(v / 2) + (1 / sum(post[, j])) * sum(post[, j] * (logw[, j] - w[, j]))), c(2.01, 100), tol = 0.0000001)$minimum
        }
      }

      if (density == "MVCN") {
        for (j in 1:k) {
          alpha[j] <- max(0.50, sum(post[, j] * w[, j]) / sum(post[, j]))

          delta <- sapply(1:num, function(i) tr(solve(sigmaUY[, , j]) %*% (Yresh[, , i] - M[, , j]) %*% solve(sigmaVY[, , j]) %*% t(Yresh[, , i] - M[, , j])))

          eta.temp <- sum(post[, j] * (1 - w[, j]) * delta) / (sum(post[, j] * (1 - w[, j])) * p * r)

          eta[j] <- max(1.0001, eta.temp)
        }
      }

      ### Density and loglik calculation ###

      if (density == "MVN") {
        for (j in 1:k) {
          dens[, j] <- dMVnorm(X = Yresh, M = M[, , j], U = sigmaUY[, , j], V = sigmaVY[, , j])
        }
      } else if (density == "MVT") {
        for (j in 1:k) {
          dens[, j] <- dMVT(X = Yresh, M = M[, , j], U = sigmaUY[, , j], V = sigmaVY[, , j], nu = nu[j])
        }
      } else if (density == "MVCN") {
        for (j in 1:k) {
          dens[, j] <- dMVcont(X = Yresh, M = M[, , j], U = sigmaUY[, , j], V = sigmaVY[, , j], alpha = alpha[j], eta = eta[j])
        }
      }

      f <- array(dens, c(num0, t, k))

      foo <- matrix(rep(prior, each = num0), ncol = k) * f[, 1, ]
      sumfoo <- rowSums(foo)
      lscale <- log(sumfoo)
      foo <- foo / sumfoo
      A[, 1, ] <- lscale + log(foo) # lalpha [ ,1]

      for (T in 2:t) {
        foo <- foo %*% PI * f[, T, ]
        sumfoo <- rowSums(foo)
        lscale <- lscale + log(sumfoo)
        foo <- foo / sumfoo
        A[, T, ] <- log(foo) + lscale
      }

      loglik.new <- sum(lscale) # log-likelihood

      B[, t, ] <- 0

      foo <- matrix(rep(rep(1 / k, k), each = num0), ncol = k)
      lscale <- log(k)

      for (T in (t - 1):1) {
        foo <- t(PI %*% t(foo * f[, T + 1, ]))
        B[, T, ] <- log(foo) + lscale
        sumfoo <- rowSums(foo)
        foo <- foo / sumfoo
        lscale <- lscale + log(sumfoo)
      }

      ll <- c(ll, loglik.new)

      # stopping rule

      if ((round(loglik.new, dec(tol)) - round(loglik.old, dec(tol))) < tol) {
        check <- 1
      }

      if (m.iter > 1) {
        if (round(loglik.new, dec(tol)) < round(loglik.old, dec(tol))) {
          mark <- 1
        } else {
          mark <- 0
        }
      }

      if (m.iter == maxit) {
        check <- 1
      }

      loglik.old <- loglik.new
    }

    #### Output ####

    # --------------------- #
    # Classification Matrix #
    # --------------------- #

    group <- apply(post, 1, which.max)
    groupT <- array(group, c(num0, t), dimnames = list(1:num0, paste("time", 1:t, sep = " ")))

    # -------------------- #
    # Information criteria #
    # -------------------- #

    # Number of parameters

    M.par <- (p * r) * k

    if (mod.row == "EII") {
      rowparY <- 1
    }
    if (mod.row == "VII") {
      rowparY <- k
    }
    if (mod.row == "EEI") {
      rowparY <- p
    }
    if (mod.row == "VEI") {
      rowparY <- k + (p - 1)
    }
    if (mod.row == "EVI") {
      rowparY <- 1 + k * (p - 1)
    }
    if (mod.row == "VVI") {
      rowparY <- k * p
    }
    if (mod.row == "EEE") {
      rowparY <- p * (p + 1) / 2
    }
    if (mod.row == "VEE") {
      rowparY <- k - 1 + p * (p + 1) / 2
    }
    if (mod.row == "EVE") {
      rowparY <- 1 + k * (p - 1) + p * (p - 1) / 2
    }
    if (mod.row == "VVE") {
      rowparY <- k * p + p * (p - 1) / 2
    }
    if (mod.row == "EEV") {
      rowparY <- p + k * p * (p - 1) / 2
    }
    if (mod.row == "VEV") {
      rowparY <- k + (p - 1) + (k * p * (p - 1) / 2)
    }
    if (mod.row == "EVV") {
      rowparY <- 1 + k * (p * ((p + 1) / 2) - 1)
    }
    if (mod.row == "VVV") {
      rowparY <- k * p * (p + 1) / 2
    }

    if (mod.col == "II") {
      colparY <- 0
    }
    if (mod.col == "EI") {
      colparY <- r - 1
    }
    if (mod.col == "VI") {
      colparY <- k * (r - 1)
    }
    if (mod.col == "EE") {
      colparY <- r * ((r + 1) / 2) - 1
    }
    if (mod.col == "VE") {
      colparY <- k * (r - 1) + r * (r - 1) / 2
    }
    if (mod.col == "EV") {
      colparY <- (r - 1) + k * r * (r - 1) / 2
    }
    if (mod.col == "VV") {
      colparY <- k * (r * ((r + 1) / 2) - 1)
    }

    weights <- k - 1
    pipar <- k * (k - 1)

    if (density == "MVN") {
      nupar <- 0
    } else if (density == "MVT") {
      nupar <- k
    } else if (density == "MVCN") {
      nupar <- 2 * k
    }

    npar <- M.par + rowparY + colparY + weights + pipar + nupar

    name <- c(mod.row, mod.col)

    # to be minimized

    BIC <- -2 * loglik.new + npar * log(num)

    ptm2 <- proc.time() - ptm
    time <- ptm2[3]

    if (density == "MVCN") {
      return(list(
        dens = density, name = name, prior = prior, M = M, U = sigmaUY, V = sigmaVY, alpha = alpha, eta = eta,
        PI = round(PI, digits = 3), loglik = loglik.new, mark = mark, check = check, npar = npar, iter = m.iter, time = time,
        BIC = BIC, class = groupT, pgood = array(w, dim = c(num0, t, k))
      ))
    } else if (density == "MVN") {
      return(list(
        dens = density, name = name, prior = prior, M = M, U = sigmaUY, V = sigmaVY,
        PI = round(PI, digits = 3), loglik = loglik.new, mark = mark, check = check, npar = npar, iter = m.iter, time = time,
        BIC = BIC, class = groupT
      ))
    } else if (density == "MVT") {
      return(list(
        dens = density, name = name, prior = prior, M = M, U = sigmaUY, V = sigmaVY, nu = nu,
        PI = round(PI, digits = 3), loglik = loglik.new, mark = mark, check = check, npar = npar, iter = m.iter, time = time,
        BIC = BIC, class = groupT, w = array(w, dim = c(num0, t, k))
      ))
    }
  }

  for (g in 1:length(k)) {
    ptm <- proc.time()
    if (verbose == TRUE) {
      print(paste(paste("Fitting Parsimonious", density), "HMMs with k =", k[g]))
    }
    cluster <- snow::makeCluster(nThreads, type = "SOCK")
    doSNOW::registerDoSNOW(cluster)

    pb <- progress::progress_bar$new(
      format = "(:spin) [:bar] :percent [Elapsed time: :elapsedfull || Estimated time remaining: :eta]",
      total = pt.mod,
      complete = "=", # Completion bar character
      incomplete = "-", # Incomplete bar character
      current = ">", # Current bar character
      width = 100
    )
    progress <- function(n) {
      pb$tick()
    }
    opts <- list(progress = progress)
    l <- 0

    oper[[g]] <- foreach::foreach(l = 1:pt.mod, .combine = "comb", .multicombine = TRUE, .init = list(list()), .options.snow = opts) %dopar% {
      res <- tryCatch(Pars.HMM(
        Y = Y, k = k[g], density = density, init.par = init.par[[1]][[g]][[1]][[init.par[[5]][l]]],
        mod.row = as.character(list.comb[l, 1]), mod.col = as.character(list.comb[l, 2]),
        tol = tol, tol2 = tol2, maxit = maxit, maxit2 = maxit2
      ), error = function(e) {
        NA
      })

      list(res)
    }

    snow::stopCluster(cluster)
    foreach::registerDoSEQ()

    ptm2 <- proc.time() - ptm
    time[g] <- ptm2[3]
  }

  return(list(results = oper, c.time = time, models = list.comb))
} # fitting function