R/mixturemcmc.R

In simonsays1980/finmix: An R package for Bayesian estimation of finite mixture distributions

## Copyright (C) 2013 Lars Simon Zehnder
#
# This file is part of finmix.
#
# finmix is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# finmix is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with finmix. If not, see <http://www.gnu.org/licenses/>.

#' Performs MCMC sampling for finite mixture models
#' 
#' @description 
#' Calling [mixturemcmc()] performs MCMC sampling on the observations stored 
#' in the `fdata` object for the finite mixture model defined in the `model` 
#' object. MCMC sampling is performed with a Gibbs sampler for all finite 
#' mixture models using a prior that must be defined in the `prior` object. 
#' There are possibilities to control the MCMC sampling by hyperparameters 
#' stored in the `mcmc` object. 
#' 
#' @details 
#' ## Performance
#' This function is the central part of the `finmix` package. For MCMC sampling 
#' it relies on high-performance C++ code using the `Rcpp` and `RcppArmadillo` 
#' packages. More specifically, these packages simplify the usage of external 
#' C++ code on the objects in `R` memory (enabled by `R`'s `C` interface). 
#' Execution of MCMC sampling with the default of 10,000 iterations and a 
#' burn-in of 1,000 iterations should finish in a few seconds. 
#' 
#' ## Algorithms
#' The algorithms used here are for the most part specified in the excellent 
#' book \emph{Finite Mixture and Markov Switching Models} by 
#' Sylvia Fr\"uwirth-Schnatter. These algorithms rely on Gibbs sampling by 
#' alternating between sampling the component and weight parameters of the 
#' finite mixture model and the indicators of the data. Thereby, a so-called 
#' random permutation is performed at each iteration of the algorithm, i.e. the 
#' indicators `S` and the component and weight parameters are permuted by their 
#' index. As explained by Fr\"uwirth-Schnatter (2006, Section 3.5.5) label 
#' switching in estimation of finite mixture distributions has to be addressed 
#' explicitly when Bayesian estimation is used. While in maximum likelihood 
#' estimation this is of no concern because only one of the equivalent modes of 
#' likelihood function needs to be found, Bayesian estimation needs to explore 
#' the full mixture posterior distribution and label switching occurs randomly, 
#' but frequently during MCMC sampling. to overcome these issues the sampler is 
#' forced to switch labels in a controlled form by randomly permuting the 
#' labels of the components. This results in a balanced label switching and as 
#' a result the sampler explores the full mixture posterior more thoroughly 
#' leading to more robust estimations. 
#' 
#' ### Starting by sampling the parameters
#' As laid out in the description of the input parameters sampling can start 
#' either by sampling the indicators using starting parameters or by sampling 
#' the parameters using starting indicators. The latter is for example applied, 
#' if indicators are fixed (because they might be known). For starting by 
#' sampling the parameters the slot `@@startpar` in the `mcmc` input argument 
#' must be set to `TRUE` (default) and starting indicators must be present in 
#' slot `@@S` of the `fdata` object. 
#' 
#' @param fdata An `fdata` object storing the observations in slot `@@y` and 
#'   the (starting) indicators in slot `@@S`. If sampling should start by 
#'   sampling the parameters the starting indicators must be defined. 
#' @param model A `model` object specifying the finite mixture model. If it 
#'   should be started by sampling the indicators starting parameters and 
#'   weights must be defined in slots `@@par` and `@@weight` respectively. 
#' @param prior A `prior` object specifying the prior distribution for Bayesian 
#'   estimation. This object must be fully specified regardless, if sampling 
#'   should start with the indicators or parameters. See [priordefine()] for 
#'   choosing automatically a data dependent prior distribution. 
#' @param mcmc An `mcmc` object storing the hyper-parameters for MCMC sampling. 
#'   If slot `@@startpar` is `TRUE` sampling starts by sampling the parameters. 
#'   Henceforth, it needs starting indicators. 
#' @return An object of class [mcmcoutput-class] storing the MCMC 
#'   sampling results. 
#' @export
#' 
#' @examples 
#' # Define a Poisson mixture model with two components.
#' f_model <- model("poisson", par = list(lambda = c(0.3, 1.2)), K = 2)
#' # Simulate data from the mixture model.
#' f_data <- simulate(f_model)
#' # Define the hyper-parameters for MCMC sampling.
#' f_mcmc <- mcmc()
#' # Complete object slots for consistency. 
#' (f_data ~ f_model ~ f_mcmc) %=% mcmcstart(f_data, f_model, f_mcmc)
#' # Define the prior distribution by relying on the data.
#' f_prior <- priordefine(f_data, f_model)
#' # Start MCMC sampling.
#' f_output <- mixturemcmc(f_data, f_model, f_prior, f_mcmc)
#' # Get the sampled model parameters.
#' getPar(f_output)
#' 
#' @seealso 
#' * [fdata-class] for the `fdata` class definition
#' * [model-class] for the `model` class definition
#' * [prior-class] for the `prior` class definition
#' * [prior()] for the `prior` class constructor
#' * [priordefine()] for the advanced class constructor of the `prior` class
#' * [mcmc-class] for the `mcmc` class definition
#' * [mcmc()] for the `mcmc` class constructor
#' * [mcmcstart()] for defining starting parameters and/or indicators
#' 
#' @references 
#' Fr\"uwirth-Schnatter, S. (2006), "Finite Mixture Models and Markov Switching 
#' Models", Springer
"mixturemcmc" <- function(fdata, model, prior, mcmc) {
  ## Check arguments
  mcmc <- .check.args.Mixturemcmc(fdata, model, prior, mcmc, nargs())

  ## Default ordering for MCMC: bycolumn
  setBycolumn(fdata) <- TRUE
  ######################### MCMC SAMPLING #############################
  ## Set the indicators as a default to one for K == 1
  if (model@K == 1) {
    fdata@S <- matrix(1, nrow = fdata@N, ncol = 1)
  }
  dist <- model@dist
  if (dist == "poisson") {
    .do.MCMC.Poisson(fdata, model, prior, mcmc)
  } else if (dist == "cond.poisson") {
    .do.MCMC.CondPoisson(fdata, model, prior, mcmc)
  } else if (dist == "binomial") {
    .do.MCMC.Binomial(fdata, model, prior, mcmc)
  } else if (dist == "exponential") {
    .do.MCMC.Exponential(fdata, model, prior, mcmc)
  } else if (dist == "normal") {
    .do.MCMC.Normal(fdata, model, prior, mcmc)
  } else if (dist == "student") {
    .do.MCMC.Student(fdata, model, prior, mcmc)
  } else if (dist == "normult") {
    .do.MCMC.Normult(fdata, model, prior, mcmc)
  } else if (dist == "studmult") {
    .do.MCMC.Studmult(fdata, model, prior, mcmc)
  }
} ## end mixturemcmc

### Private functions
### These functions are not exported

#' Checks input arguments for MCMC sampling
#' 
#' @description 
#' For internal usage only. This function checks if the input arguments passed 
#' in to [mixturemcmc()] are valid, i.e. `fdata` must contain valid data in 
#' slot `@@y` and in case of starting with sampling the parameters indicators 
#' in slot `@@S`. Furthermore, the data in slot `@@y` must match with the 
#' specified distribution in `@@dist` of the `model` object. 
#' If MCMC sampling should start by sampling the indicators, the `model` object 
#' must contain valid starting parameters in slots `@@par` and `@@weight`. 
#' The `prior` object must contain valid parameters for the prior distribution. 
#' Finally, if a fixed indicator model is used, `@@startpar` in `mcmc` must be 
#' `TRUE` and `@@ranperm` must be `FALSE`.
#' 
#' In addition this function checks for consistency between the four input 
#' objects and modifies the hyper-parameters in the `mcmc` object accordingly 
#' for the user.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the finite mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution. 
#' @param mcmc.obj An `mcmc` object defining the hyper-parameters for MCMC 
#'   sampling.
#' @param n.args An integer specifying how many arguments have been provided 
#'   by the user. As all arguments must be provided values below four throw an 
#'   error. 
#' @return An object of class [mcmc-class]. If any check does not pass an 
#'   error is thrown to let the user know, why MCMC sampling cannot be 
#'   performed with the actual setting.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".check.args.Mixturemcmc" <- function(fdata.obj, model.obj,
                                      prior.obj, mcmc.obj, n.args) {
  ## Check if all arguments are provided
  if (n.args < 4) {
    stop("All arguments must be provided.", call. = FALSE)
  }
  ## Check if 'fdata' object is valid
  if (class(fdata.obj) != "fdata") {
    stop(paste("Unkown argument. Argument 1 must be an ",
      "object of class 'fdata'.",
      sep = ""
    ), call. = FALSE)
  }
  hasY(fdata.obj, verbose = TRUE)
  ## Check if 'model' was provided:
  if (class(model.obj) != "model") {
    stop(paste("Unknown argument. Argument 2 must be an ",
      "object of class 'model'.",
      sep = ""
    ), call. = FALSE)
  }
  ## Check if 'prior' was provided:
  if (class(prior.obj) != "prior") {
    stop(paste("Unknown argument. Argument 3 must be an ",
      "object of class 'prior'.",
      sep = ""
    ), call. = FALSE)
  }
  ## Check if 'mcmc' was provided:
  if (class(mcmc.obj) != "mcmc") {
    stop(paste("Unkown argument. Argument 4 must be an ",
      "object of class 'mcmc'.",
      sep = ""
    ), call. = FALSE)
  }
  ## Check if @startpar in 'mcmc' object and @indicfix in
  ## 'model' object match.
  ## For fixed indicator models indicators are not sampled.
  if (model.obj@indicfix && !mcmc.obj@startpar) {
    mcmc.obj@startpar <- TRUE
  }
  ## Check if @K in 'model' object is one. For a model with
  ## only one component indicators are not sampled.
  if (model.obj@K == 1) {
    mcmc.obj@startpar <- TRUE
  }
  ## If @startpar in 'mcmc.obj' is TRUE, it should be started
  ## by sampling the parameters. In this case starting
  ## indicators must be provided in the 'fdata.obj' object.
  ## If @startpar in 'mcmc.obj' is FALSE it should be started
  ## by sampling the indicators. In this case starting
  ## parameters must be provided in the 'model.obj' object.
  if (model.obj@K > 1) {
    if (mcmc.obj@startpar) {
      if (!hasS(fdata.obj)) {
        stop(paste("For starting with sampling the parameters ",
          "the 'fdata' object must contain starting ",
          "indicator values. See ?mcmcstart for ",
          "generating valid starting values.",
          sep = ""
        ),
        call. = FALSE
        )
      }
    } else {
      if (!hasPar(model.obj)) {
        stop(paste("For starting with sampling the indicators ",
          "the 'model' object must contain starting ",
          "parameter values. See ?mcmcstart for ",
          "generating valid starting values.",
          sep = ""
        ),
        call. = FALSE
        )
      }
      if (!hasWeight(model.obj)) {
        stop(paste("For starting with sampling the indicators ",
          "the 'model' object must contain starting ",
          "weight values. See ?mcmcstart for ",
          "generating valid starting values.",
          sep = ""
        ),
        call. = FALSE
        )
      }
    }
  }
  ## Check if 'fdata' object and 'model' objects match
  ## Call '.check.fdata.model.Mcmcstart()' from 'mcmcstart.R'.
  .check.fdata.model.Mcmcstart(fdata.obj, model.obj)
  ## Check if 'prior' object is valid
  if (!hasPriorPar(prior.obj, model.obj)) {
    stop(paste("Slot @par in 'prior' object is empty. ",
      "For MCMC sampling the prior needs fully ",
      "specified parameters. See ?priordefine for ",
      "generating valid prior parameters.",
      sep = ""
    ),
    call. = FALSE
    )
  }
  if (!model.obj@indicfix && model.obj@K > 1) {
    if (!hasPriorWeight(prior.obj, model.obj)) {
      stop(paste("Slot @weight of 'prior' object is empty. ",
        "For MCMC sampling the prior needs specified ",
        "parameters for the prior of the weights. See ",
        "?priordefine for generating valid prior ",
        "parameters.",
        sep = ""
      ), call. = FALSE)
    }
  }
  ## Check if @indicfix in 'model' object and
  ## @ranperm in 'mcmc' object match.
  ## For a fixed indicator model random permutation
  ## sampling is senseless.
  if (model.obj@indicfix && mcmc.obj@ranperm) {
    mcmc.obj@ranperm <- FALSE
  }
  ## For a model with only one component random permutation
  ## is senseless as well.
  if (model.obj@K == 1 && mcmc.obj@ranperm) {
    mcmc.obj@ranperm <- FALSE
  }
  return(mcmc.obj)
}

#' Checks validity of repetitions for MCMC sampling
#' 
#' @description 
#' For a Binomial model either the `fdata` object or the `model` object must 
#' have specified repetitions in slot `@@T`. This can be either a `matrix` 
#' object of dimension `N x 1` or `1 x 1` (if all repetitions are the 
#' same).
#' 
#' @param data An `fdata` object containing the data.
#' @param model A `model` object specifying the finite mixture model.
#' @return None. If any check does not pass an error is thrown to inform the 
#'   user of the detected inconsistency. 
#' @noRd
#' 
#' @seealso 
#' * [fdata-class] for the `fdata` class definition
#' * [model-class] for the `model` class definition
".valid.Reps.Binomial" <- function(data, model) {
  has.reps <- !all(is.na(data@T))
  N        <- data@N
  if (has.reps) {
    if (data@bycolumn) {
      if (nrow(data@T) != N && nrow(data@T) != 1) {
        stop(paste("Number of repetitions in slot @T of 'data' object ",
          "does not match number of observations in slot @N.",
          sep = ""
        ), call. = FALSE)
      } else if (nrow(data@T) == N) {
        T <- data@T
      } else { ## dimension of T is 1 x 1
        T <- matrix(data@T[1, 1], nrow = N, ncol = 1)
      }
    } else { ## data stored by row
      if (ncol(data@T) != N && ncol(data@T) != 1) {
        stop(paste("Number of repetitions in slot @T of 'data' object ",
          "does not match number of observations slot @N.",
          sep = ""
        ), call. = FALSE)
      } else if (ncol(data@T) == N) {
        T <- t(data@T)
      } else { ## dimension of T is 1 x 1
        T <- matrix(data@T[1, 1], nrow = N, ncol = 1)
      }
    }
  } else { ## then check in model
    has.reps <- !all(is.na(model@T))
    if (has.reps) {
      if (nrow(model@T) != N && nrow(model@T) != 1) {
        stop(paste("Neither 'data' nor 'model' has correctly ",
          "specified repetitions in slot @T for a binomial model.",
          sep = ""
        ), call. = FALSE)
      } else if (nrow(model@T) == N) {
        T <- model@T
      } else { ## dimension of T is 1 x 1
        T <- matrix(model@T[1, 1], nrow = N, ncol = 1)
      }
    } else {
      stop(paste("Neither 'data' object nor 'model' object has ",
        "repetitions in slot @T for a binomial model specified.",
        sep = ""
      ), call. = FALSE)
    }
  }
  ## Check for identifiability ##
  ## Reference: Teicher (1961) ##
  rep.occ <- table(T)
  if (dim(unique(T))[1] == 1) {
    if (T[1, 1] < 2 * model@K - 1) {
      warning(paste("This binomial mixture model is not identifiable. ",
        "For equal repetitions in slot @T it must hold T >= 2K - 1. ",
        "See Teicher (1961) for reference.",
        sep = ""
      ), call. = FALSE)
    }
  } else {
    if (length(rep.occ) != nrow(T)) {
      if (all(dimnames(rep.occ)$T < rep.occ - 1)) {
        warning(paste("This binomial mixture model is not identifiable. ",
          "For varying repetitions 'T_i' in slot @T it must hold T_h ",
          "> r_h - 1, for unique repetitions 'T_h' and their ",
          "respective occurences 'r_h'. See Teicher (1961) ",
          "for reference",
          sep = ""
        ), call. = FALSE)
      } else {
        diff <- diff(sort(unique(T)))
        if (any(diff < rep.occ[1:(length(diff))])) {
          warning(paste("This binomial mixture model is not identifiable. ",
            "For varying repetitions 'T_i' in slot @T it must hold T_h ",
            "- T_(h+1) >= r_h for unique repetitions 'T_h' and ",
            "respective occurrences 'r_h'. See Teicher (1961) ",
            "for reference.",
            sep = ""
          ), call. = FALSE)
        }
      }
    }
  }
}

### MCMC
### For each model the MCMC output has to be prepared

#' Perform MCMC sampling for Poisson mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for Poisson mixture models regarding the specifications in the 
#' passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Poisson" <- function(fdata.obj, model.obj, prior.obj, mcmc.obj) {
  ## Base slots inherited to every derived class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  ranperm <- mcmc.obj@ranperm
  burnin <- mcmc.obj@burnin
  ## Set for MCMC default exposures:
  if (!hasExp(fdata.obj)) {
    fdata.obj@exp <- matrix(1, nrow = N, ncol = 1)
  }
  pars <- list(lambda = array(numeric(), dim = c(M, K)))
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    post.a <- array(numeric(), dim = c(M, K))
    post.b <- array(numeric(), dim = c(M, K))
    post.par <- list(a = post.a, b = post.b)
    posts <- list(par = post.par)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfix(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputfixpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior ##
      hypers <- list(b = array(numeric(), dim = c(M, 1)))
      ## Model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfixhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputfixhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers, post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end hier
    }
    ## end indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## Model with simulated indicators ##
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik,
      mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$S)
    }
    ## Model with simple prior ##
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputbase(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## model with hierarchical prior ##
      hypers <- list(b = array(numeric(), dim = c(M, 1)))
      ## model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      } else {
        ## model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_poisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      }
    } ## end hier
  } ## end no indicfix
}

#' Perform MCMC sampling for Binomial mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for Binomial mixture models regarding the specifications in the 
#' passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Binomial" <- function(fdata.obj, model.obj, prior.obj, mcmc.obj) {
  ## Base slots inherited to every derived 'mcmcoutput' class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  burnin <- mcmc.obj@burnin
  ranperm <- mcmc.obj@ranperm
  pars <- list(p = array(numeric(), dim = c(c(M, K))))
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    post.a <- array(numeric(), dim = c(M, K))
    post.b <- array(numeric(), dim = c(M, K))
    post.par <- list(a = post.a, b = post.b)
    posts <- list(par = post.par)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    if (!mcmc.obj@storepost) {
      mcmcout <- .mcmcoutputfix(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, model = model.obj,
        prior = prior.obj
      )
      .Call("mcmc_binomial_cc", fdata.obj, model.obj, prior.obj,
        mcmc.obj, mcmcout,
        PACKAGE = "finmix"
      )
      return(mcmcout)
    } else {
      ## MCMC output with posterior hyper parameters stored
      mcmcout <- .mcmcoutputfixpost(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, post = posts,
        model = model.obj, prior = prior.obj
      )
      .Call("mcmc_binomial_cc", fdata.obj, model.obj, prior.obj,
        mcmc.obj, mcmcout,
        PACKAGE = "finmix"
      )
      return(mcmcout)
    }
    ## End: indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## Model with simulated indicators
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik, mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$S)
    }
    if (!mcmc.obj@storepost) {
      mcmcout <- .mcmcoutputbase(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, weight = weights,
        entropy = entropies, ST = STm, S = Sm,
        NK = NKm, clust = clustm,
        model = model.obj, prior = prior.obj
      )
      .Call("mcmc_binomial_cc", fdata.obj, model.obj, prior.obj, mcmc.obj,
        mcmcout,
        PACKAGE = "finmix"
      )
      if (mcmc.obj@storeS == 0) {
        mcmcout@S <- as.array(NA)
      }
      return(mcmcout)
    } else {
      ## MCMC output with posterior hyper parameters stored
      mcmcout <- .mcmcoutputpost(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, weight = weights,
        entropy = entropies, ST = STm, S = Sm,
        NK = NKm, clust = clustm, post = posts,
        model = model.obj, prior = prior.obj
      )
      .Call("mcmc_binomial_cc", fdata.obj, model.obj, prior.obj, mcmc.obj,
        mcmcout,
        PACKAGE = "finmix"
      )
      if (mcmc.obj@storeS == 0) {
        mcmcout@S <- as.array(NA)
      }
      return(mcmcout)
    }
  } ## End no indicfix
}

#' Perform MCMC sampling for exponential mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for exponential mixture models regarding the specifications in the 
#' passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Exponential" <- function(fdata.obj, model.obj, prior.obj, mcmc.obj) {
  # Base slots inherited to each derived class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  ranperm <- FALSE
  burnin <- mcmc.obj@burnin
  pars <- list(lambda = array(numeric(), dim = c(M, K)))
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    post.a <- array(numeric(), dim = c(M, K))
    post.b <- array(numeric(), dim = c(M, K))
    post.par <- list(a = post.a, b = post.b)
    posts <- list(par = post.par)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  # Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    # Model with simple prior
    if (!mcmc.obj@storepost) {
      # Model output with NO posterior parameters stored
      mcmcout <- .mcmcoutputfix(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, model = model.obj,
        prior = prior.obj
      )
    } else {
      # Model output with posterior parameters stored
      mcmcout <- .mcmcoutputfixpost(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, post = posts,
        model = model.obj, prior = prior.obj
      )
    }
  } else if (!model.obj@indicfix && K > 1) {
    # Model with simulated indicators
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik, mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      # First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      SM[, 1] <- as.integer(datac$S)
    }
    # Model with simple prior
    # Model output with NO posterior parameters stored
    if (!mcmc.obj@storepost) {
      mcmcout <- .mcmcoutputbase(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, weight = weights,
        entropy = entropies, ST = STm, S = Sm,
        NK = NKm, clust = clustm, model = model.obj,
        prior = prior.obj
      )
    } else {
      # Model output with posterior parameters stored
      mcmcout <- .mcmcoutputpost(
        M = M, burnin = burnin, ranperm = ranperm,
        par = pars, log = logs, weight = weights,
        entropy = entropies, ST = STm, S = Sm,
        NK = NKm, clust = clustm, post = posts,
        model = model.obj, prior = prior.obj
      )
    }
  }
  .Call("mcmc_exponential_cc", fdata.obj, model.obj, prior.obj, mcmc.obj,
    mcmcout,
    PACKAGE = "finmix"
  )
  if (mcmc.obj@storeS == 0) {
    mcmcout@S <- as.array(NA)
  }
  return(mcmcout)
}

#' Perform MCMC sampling for conditional Poisson mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for conditional Poisson mixture models regarding the specifications 
#' in the passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.CondPoisson" <- function(fdata.obj, model.obj, prior.obj, mcmc.obj) {
  if (nrow(fdata.obj@exp) == 1) {
    if (is.na(fdata.obj@exp)) {
      fdata.obj@exp <- matrix(1, nrow = fdata.obj@N, ncol = 1)
    } else {
      fdata.obj@exp <- matrix(fdata.obj@exp, nrow = fdata.obj@N, ncol = 1)
    }
  }
  if (mcmc.obj@ranperm) {
    mcmc.obj@ranperm <- FALSE
  }
  ## base slots inherited to every derived class ##
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  burnin <- mcmc.obj@burnin
  ranperm <- mcmc.obj@ranperm
  pars <- list(
    lambda = array(numeric(), dim = c(M, K)),
    acc = 0.0
  )
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    post.Q <- array(numeric(), dim = c(M, K))
    post.N <- array(numeric(), dim = c(M, K))
    post.par <- list(Q = post.Q, N = post.N)
    posts <- list(par = post.par)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## model with fixed indicators ##
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    ## model with simple prior ##
    if (!prior.obj@hier) {
      ## model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfix(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        ## model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputfixpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## model with hierarchical prior ##
      hypers <- list(b = array(numeric(), dim = c(M, 1)))
      ## model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfixhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        ## model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputfixhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end hier
    }
    ## end indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## model with simulated indicators ##
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik, mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$S)
    }
    ## model with simple prior ##
    if (!prior.obj@hier) {
      ## model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputbase(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      } else {
        ## model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## model with hierarchical prior ##
      hypers <- list(b = array(numeric(), dim = c(M, 1)))
      ## model output with NO posterior parameters stored ##
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      } else {
        ## model output with posterior parameters stored ##
        mcmcout <- .mcmcoutputhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_condpoisson_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(NA)
        }
        return(mcmcout)
      }
    } ## end hier
  } ## end no indicfix
}

#' Perform MCMC sampling for normal mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for normal mixture models regarding the specifications in the 
#' passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Normal" <- function(fdata.obj, model.obj, prior.obj,
                              mcmc.obj) {
  ## Base slots inherited to each derived class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  ranperm <- mcmc.obj@ranperm
  burnin<- mcmc.obj@burnin
  ## Set for MCMC default exposures:
  pars <- list(
    mu = array(numeric(), dim = c(M, K)),
    sigma = array(numeric(), dim = c(M, K))
  )
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    b.post <- array(numeric(), dim = c(M, K))
    B.post <- array(numeric(), dim = c(M, K))
    mu.post <- list(b = b.post, B = B.post)
    c.post <- array(numeric(), dim = c(M, K))
    C.post <- array(numeric(), dim = c(M, K))
    sigma.post <- list(c = c.post, C = C.post)
    par.post <- list(mu = mu.post, sigma = sigma.post)
    posts <- list(par = par.post)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfix(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        mcmcout <- .mcmcoutputfixpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, 1)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfixhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputfixhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers, post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end hier
    }
    ## end indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## Model with simulated indicators
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik,
      mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$s)
    }
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputbase(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored
        mcmcout <- .mcmcoutputpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, 1)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(is.integer(NA))
        }
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          post = posts, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normal_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      }
    } ## end hier
  } ## end no indicfix
}

#' Perform MCMC sampling for Student-t mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for Student-t mixture models regarding the specifications in the 
#' passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Student" <- function(fdata.obj, model.obj, prior.obj,
                               mcmc.obj) {
  ## Base slots inherited to each derived class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  ranperm <- mcmc.obj@ranperm
  burnin <- mcmc.obj@burnin
  ## Set for MCMC default exposures:
  pars <- list(
    mu = array(numeric(), dim = c(M, K)),
    sigma = array(numeric(), dim = c(M, K)),
    df = array(numeric(), dim = c(M, K)),
    acc = array(0.0, dim = c(1, K))
  )
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    b.post <- array(numeric(), dim = c(M, K))
    B.post <- array(numeric(), dim = c(M, K))
    mu.post <- list(b = b.post, B = B.post)
    c.post <- array(numeric(), dim = c(M, K))
    C.post <- array(numeric(), dim = c(M, K))
    sigma.post <- list(c = c.post, C = C.post)
    par.post <- list(mu = mu.post, sigma = sigma.post)
    posts <- list(par = par.post)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfix(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        mcmcout <- .mcmcoutputfixpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, 1)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfixhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputfixhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers, post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        return(mcmcout)
      }
      ## end hier
    }
    ## end indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## Model with simulated indicators
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik,
      mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$s)
    }
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputbase(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored
        mcmcout <- .mcmcoutputpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, 1)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          post = posts, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_student_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        return(mcmcout)
      }
    } ## end hier
  } ## end no indicfix
}

#' Perform MCMC sampling for multivariate normal mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for multivariate normal mixture models regarding the specifications 
#' in the passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Normult" <- function(fdata.obj, model.obj, prior.obj, mcmc.obj) {
  ## Base slots inherited to each derived class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  ranperm <- mcmc.obj@ranperm
  burnin <- mcmc.obj@burnin

  ## Constants simplifying construction
  r <- fdata.obj@r
  s <- r * (r + 1) / 2

  ## Set for MCMC default expousres
  pars <- list(
    mu = array(numeric(), dim = c(M, r, K)),
    sigma = array(numeric(), dim = c(M, s, K)),
    storeinv = mcmc.obj@storeinv
  )
  if (mcmc.obj@storeinv) {
    pars$sigmainv <- array(numeric(), dim = c(M, s, K))
  }
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    b.post <- array(numeric(), dim = c(M, r, K))
    B.post <- array(numeric(), dim = c(M, s, K))
    mu.post <- list(b = b.post, B = B.post)
    c.post <- array(numeric(), dim = c(M, K))
    C.post <- array(numeric(), dim = c(M, s, K))
    sigma.post <- list(c = c.post, C = C.post)
    par.post <- list(mu = mu.post, sigma = sigma.post)
    posts <- list(par = par.post)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfix(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        mcmcout <- .mcmcoutputfixpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, s)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfixhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputfixhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers, post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
      ## end hier
    }
    ## end indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## Model with simulated indicators
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik,
      mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$s)
    }
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputbase(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored
        mcmcout <- .mcmcoutputpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, s)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          post = posts, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_normult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
    } ## end hier
  } ## end no indicfix
}

#' Perform MCMC sampling for multivariate Student-t mixtures
#' 
#' @description 
#' For internal usage only. This function prepares all data containers for MCMC 
#' sampling for multivariate Student-t mixture models regarding the 
#' specifications in the passed-in objects.
#' 
#' @param fdata.obj An `fdata` object containing the data.
#' @param model.obj A `model` object specifying the mixture model.
#' @param prior.obj A `prior` object specifying the prior distribution.
#' @param mcmc.obj An `mcmc` object cotnaining the hyper-parameters for MCMC 
#'   sampling.
#' @param An object of class [mcmcoutput-class] containing the 
#'   results of MCMC sampling.
#' @noRd
#' 
#' @seealso 
#' * [mixturemcmc()] for the calling function
".do.MCMC.Studmult" <- function(fdata.obj, model.obj, prior.obj, mcmc.obj) {
  ## Base slots inherited to each derived class
  K <- model.obj@K
  N <- fdata.obj@N
  M <- mcmc.obj@M
  ranperm <- mcmc.obj@ranperm
  burnin <- mcmc.obj@burnin

  ## Constants simplifying construction
  r <- fdata.obj@r
  s <- r * (r + 1) / 2

  ## Set for MCMC default expousres
  pars <- list(
    mu = array(numeric(), dim = c(M, r, K)),
    sigma = array(numeric(), dim = c(M, s, K)),
    df = array(numeric(), dim = c(M, K)),
    acc = array(0.0, dim = c(1, K)),
    storeinv = mcmc.obj@storeinv
  )
  if (mcmc.obj@storeinv) {
    pars$sigmainv <- array(numeric(), dim = c(M, s, K))
  }
  log.mixlik <- array(numeric(), dim = c(M, 1))
  log.mixprior <- array(numeric(), dim = c(M, 1))
  if (mcmc.obj@storepost) {
    b.post <- array(numeric(), dim = c(M, r, K))
    B.post <- array(numeric(), dim = c(M, s, K))
    mu.post <- list(b = b.post, B = B.post)
    c.post <- array(numeric(), dim = c(M, K))
    C.post <- array(numeric(), dim = c(M, s, K))
    sigma.post <- list(c = c.post, C = C.post)
    par.post <- list(mu = mu.post, sigma = sigma.post)
    posts <- list(par = par.post)
    if (!model.obj@indicfix) {
      posts$weight <- array(numeric(), dim = c(M, K))
    }
  }
  ## Model with fixed indicators
  if (model.obj@indicfix || K == 1) {
    logs <- list(mixlik = log.mixlik, mixprior = log.mixprior)
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model output with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfix(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        mcmcout <- .mcmcoutputfixpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, s)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputfixhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputfixhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          hyper = hypers, post = posts,
          model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
      ## end hier
    }
    ## end indicfix
  } else if (!model.obj@indicfix && K > 1) {
    ## Model with simulated indicators
    log.cdpost <- array(numeric(), dim = c(M, 1))
    logs <- list(
      mixlik = log.mixlik,
      mixprior = log.mixprior,
      cdpost = log.cdpost
    )
    weights <- array(numeric(), dim = c(M, K))
    entropies <- array(numeric(), dim = c(M, 1))
    STm <- array(integer(), dim = c(M, 1))
    Sm <- array(integer(), dim = c(N, mcmc.obj@storeS))
    NKm <- array(integer(), dim = c(M, K))
    clustm <- array(integer(), dim = c(N, 1))
    if (!mcmc.obj@startpar) {
      ## First sample for the indicators
      datac <- dataclass(fdata.obj, model.obj, simS = TRUE)
      Sm[, 1] <- as.integer(datac$s)
    }
    ## Model with simple prior
    if (!prior.obj@hier) {
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputbase(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        ## Model output with posterior parameters stored
        mcmcout <- .mcmcoutputpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, post = posts,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
      ## end no hier
    } else {
      ## Model with hierarchical prior
      hypers <- list(C = array(numeric(), dim = c(M, s)))
      ## Model with NO posterior parameters stored
      if (!mcmc.obj@storepost) {
        mcmcout <- .mcmcoutputhier(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          model = model.obj, prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      } else {
        ## Model with posterior parameters stored
        mcmcout <- .mcmcoutputhierpost(
          M = M, burnin = burnin,
          ranperm = ranperm,
          par = pars, log = logs,
          weight = weights,
          entropy = entropies,
          ST = STm, S = Sm, NK = NKm,
          clust = clustm, hyper = hypers,
          post = posts, model = model.obj,
          prior = prior.obj
        )
        .Call("mcmc_studmult_cc", fdata.obj, model.obj, prior.obj,
          mcmc.obj, mcmcout,
          PACKAGE = "finmix"
        )
        if (mcmc.obj@storeS == 0) {
          mcmcout@S <- as.array(as.integer(NA))
        }
        mcmcout@par$storeinv <- NULL
        return(mcmcout)
      }
    } ## end hier
  } ## end no indicfix
}