diff --git a/DESCRIPTION b/DESCRIPTION
index 1f25d4a..5c57910 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -8,8 +8,8 @@ Description: Tool-set to support Bayesian evidence synthesis.  This
     for details on applying this package while Neuenschwander et al. (2010)
     <doi:10.1177/1740774509356002> and Schmidli et al. (2014)
     <doi:10.1111/biom.12242> explain details on the methodology.
-Version: 1.6-7
-Date: 2023-06-26
+Version: 1.7-0
+Date: 2023-07-19
 Authors@R: c(person("Novartis", "Pharma AG", role = "cph")
             ,person("Sebastian", "Weber", email="sebastian.weber@novartis.com", role=c("aut", "cre"))
             ,person("Beat", "Neuenschwander", email="beat.neuenschwander@novartis.com", role="ctb")
@@ -36,7 +36,8 @@ Imports:
     dplyr,
     stats,
     utils,
-    matrixStats
+    matrixStats,
+    abind
 LinkingTo:
     BH (>= 1.72.0),
     Rcpp (>= 0.12.0),
diff --git a/NAMESPACE b/NAMESPACE
index 4b34a26..9eca642 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -18,6 +18,7 @@ S3method(dmix,betaMix)
 S3method(dmix,default)
 S3method(dmix,gammaMix)
 S3method(dmix,gammaPoissonMix)
+S3method(dmix,mvnormMix)
 S3method(dmix,normMix)
 S3method(ess,betaMix)
 S3method(ess,default)
@@ -41,11 +42,13 @@ S3method(oc2S,normMix)
 S3method(plot,EM)
 S3method(plot,gMAP)
 S3method(plot,mix)
+S3method(plot,mvnormMix)
 S3method(pmix,betaBinomialMix)
 S3method(pmix,betaMix)
 S3method(pmix,default)
 S3method(pmix,gammaMix)
 S3method(pmix,gammaPoissonMix)
+S3method(pmix,mvnormMix)
 S3method(pmix,normMix)
 S3method(pos1S,betaMix)
 S3method(pos1S,default)
@@ -62,10 +65,12 @@ S3method(postmix,normMix)
 S3method(preddist,betaMix)
 S3method(preddist,default)
 S3method(preddist,gammaMix)
+S3method(preddist,mvnormMix)
 S3method(preddist,normMix)
 S3method(predict,gMAP)
 S3method(print,EMbmm)
 S3method(print,EMgmm)
+S3method(print,EMmvnmm)
 S3method(print,EMnmm)
 S3method(print,betaBinomialMix)
 S3method(print,betaMix)
@@ -78,23 +83,27 @@ S3method(print,gammaExpMix)
 S3method(print,gammaMix)
 S3method(print,gammaPoissonMix)
 S3method(print,mix)
+S3method(print,mvnormMix)
 S3method(print,normMix)
 S3method(qmix,betaBinomialMix)
 S3method(qmix,betaMix)
 S3method(qmix,default)
 S3method(qmix,gammaMix)
 S3method(qmix,gammaPoissonMix)
+S3method(qmix,mvnormMix)
 S3method(qmix,normMix)
 S3method(rmix,betaBinomialMix)
 S3method(rmix,betaMix)
 S3method(rmix,default)
 S3method(rmix,gammaMix)
 S3method(rmix,gammaPoissonMix)
+S3method(rmix,mvnormMix)
 S3method(rmix,normMix)
 S3method(robustify,betaMix)
 S3method(robustify,default)
 S3method(robustify,gammaMix)
 S3method(robustify,normMix)
+S3method(sigma,mvnormMix)
 S3method(sigma,normMix)
 S3method(summary,betaBinomialMix)
 S3method(summary,betaMix)
@@ -102,10 +111,12 @@ S3method(summary,gMAP)
 S3method(summary,gMAPpred)
 S3method(summary,gammaMix)
 S3method(summary,gammaPoissonMix)
+S3method(summary,mvnormMix)
 S3method(summary,normMix)
 S3method(support,betaMix)
 S3method(support,default)
 S3method(support,gammaMix)
+S3method(support,mvnormMix)
 S3method(support,normMix)
 export("likelihood<-")
 export("sigma<-")
@@ -127,6 +138,7 @@ export(mixbeta)
 export(mixcombine)
 export(mixfit)
 export(mixgamma)
+export(mixmvnorm)
 export(mixnorm)
 export(mn2beta)
 export(mn2gamma)
@@ -151,6 +163,7 @@ export(robustify)
 export(sigma)
 import(Formula)
 import(Rcpp)
+import(abind)
 import(assertthat)
 import(checkmate)
 import(ggplot2)
diff --git a/NEWS.md b/NEWS.md
index bb59454..d39b42e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,16 @@
+# RBesT 1.7-0 - July 19th, 2023
+
+## Enhancements
+
+* Implementation multivariate normal mixture support in a first
+  version. This includes density evaluation, basic summary functions
+  and multivariate normal EM fitting. Support is not yet as complete
+  as for other densites, but will be expanded in upcoming releases.
+* Change the default for the option `constrain_gt1` of the EM for beta
+  mixtures to `TRUE`. This will be default constrain the fitted
+  parameters of each beta mixture component to be greater than unity
+  as required for finite ESS elir calculations.
+
 # RBesT 1.6-7 - June 26th, 2023
 
 ## Bug fixes
diff --git a/R/EM_bmm_ab.R b/R/EM_bmm_ab.R
index 7c021fb..92c2cb2 100644
--- a/R/EM_bmm_ab.R
+++ b/R/EM_bmm_ab.R
@@ -1,16 +1,10 @@
 
 ## EM for Beta Mixture Models (BMM) with Nc components
 
-EM_bmm_ab <- function(x, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,a=0.005,b=0.005), constrain_gt1)
-{
+EM_bmm_ab <- function(x, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,a=0.005,b=0.005), constrain_gt1=TRUE) {
     N <- length(x)
     assert_that(N+Nc >= Ninit)
 
-    if(missing(constrain_gt1)) {
-        message("Since version 1.6-0 of RBesT the EM for beta mixtures constrains a > 1 & b > 1 by default.\nUse constrain_gt1=FALSE for an unconstrained fit. To avoid this message use constrain_gt1=TRUE.\nIn a future version the new default will be used without this message.")
-        constrain_gt1 <- TRUE
-    }
-
     assert_logical(constrain_gt1, any.missing=FALSE, len=1)
 
     ## check data for 0 and 1 values which are problematic, but may be
diff --git a/R/EM_gmm.R b/R/EM_gmm.R
index 52763ec..0473a55 100644
--- a/R/EM_gmm.R
+++ b/R/EM_gmm.R
@@ -1,6 +1,5 @@
 ## EM for GMM with Nc components
-EM_gmm <- function(x, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(weight=0.005,alpha=0.005,beta=0.005))
-{
+EM_gmm <- function(x, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(weight=0.005,alpha=0.005,beta=0.005)) {
     N <- length(x)
     assert_that(N+Nc >= Ninit)
 
diff --git a/R/EM_mnmm.R b/R/EM_mnmm.R
index bdae935..0470fc2 100644
--- a/R/EM_mnmm.R
+++ b/R/EM_mnmm.R
@@ -4,8 +4,7 @@
 ## cluster weights and covariance matrices (taken from the knn
 ## determined clusters)
 
-EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005))
-{
+EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005)) {
     ## in case X is no matrix, interpret it as uni-variate case
     if(!is.matrix(X))
         X <- matrix(X,ncol=1)
@@ -175,8 +174,10 @@ EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol
     ##    tauEst <- sqrt(as.vector(covEst))
     ##}
 
-    mixEst <- list(p=pEst, mean=muEst, sigma=covEst)
+    ##mixEst <- list(p=pEst, mean=muEst, sigma=covEst)
 
+    mixEst <- do.call(mixmvnorm, lapply(1:Nc, function(i) c(pEst[i], muEst[i,,drop=FALSE], matrix(covEst[i,,], Nd, Nd))))
+    
     ## give further details
     attr(mixEst, "df") <- df
     attr(mixEst, "nobs") <- N
@@ -187,25 +188,22 @@ EM_mnmm <- function(X, Nc, mix_init, Ninit=50, verbose=FALSE, Niter.max=500, tol
     attr(mixEst, "tol") <- tol
     attr(mixEst, "traceLli") <- traceLli
     attr(mixEst, "traceMix") <- traceMix
+    attr(mixEst, "x") <- X
+
+    class(mixEst) <- c("EM", "EMmvnmm", "mvnormMix", "mix")
 
-    invisible(mixEst)
+    mixEst
 }
 
 ## uni-variate case
-EM_nmm <- function(X, Nc, mix_init, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005))
-{
-    res <- EM_mnmm(X=X, Nc=Nc, mix_init=mix_init, verbose=verbose, Niter.max=Niter.max, tol=tol, Neps=Neps, eps=eps)
-    convert <- function(est) suppressWarnings(do.call(mixnorm, lapply(1:Nc, function(i) c(est$p[i], est$mean[i,1], sqrt(est$sigma[i,,])))))
-    mixEst <- convert(res)
-    attr(mixEst, "df") <- attr(res, "df")
-    attr(mixEst, "nobs") <- attr(res, "nobs")
-    attr(mixEst, "lli") <- attr(res, "lli")
-    attr(mixEst, "Nc") <- attr(res, "Nc")
-    attr(mixEst, "tol") <- attr(res, "tol")
-    attr(mixEst, "traceLli") <- attr(res, "traceLli")
-    attr(mixEst, "traceMix") <-  lapply(attr(res, "traceMix"), convert)
-    attr(mixEst, "x") <- X
+EM_nmm <- function(X, Nc, mix_init, verbose=FALSE, Niter.max=500, tol, Neps, eps=c(w=0.005,m=0.005,s=0.005)) {
+    if(is.matrix(X)) {
+        assert_matrix(X, any.missing=FALSE, ncols=1)
+    }
+    mixEst <- EM_mnmm(X=X, Nc=Nc, mix_init=mix_init, verbose=verbose, Niter.max=Niter.max, tol=tol, Neps=Neps, eps=eps)
+    rownames(mixEst) <- c("w", "m", "s")
     class(mixEst) <- c("EM", "EMnmm", "normMix", "mix")
+    likelihood(mixEst) <- "normal"
     mixEst
 }
 
@@ -215,6 +213,12 @@ print.EMnmm <- function(x, ...) {
     NextMethod()
 }
 
+#' @export
+print.EMmvnmm <- function(x, ...) {
+    cat("EM for Multivariate Normal Mixture Model\nLog-Likelihood = ", logLik(x), "\n\n",sep="")
+    NextMethod()
+}
+
 ## given a vector of means and a covariance matrix for a multi-variate
 ## normal (and optionally a vector of df), return a vector of
 ## coefficients with a deterministic mapping
diff --git a/R/EM_plot.R b/R/EM_plot.R
index 6231415..e5f4c53 100644
--- a/R/EM_plot.R
+++ b/R/EM_plot.R
@@ -52,6 +52,7 @@
 #' @export
 plot.EM <- function(x, size=1.25, link=c("identity", "logit", "log"), ...) {
     pl <- list()
+    if(inherits(x, "mvnormMix")) stop("Multivariate normal mixtures plotting not supported.")
     pl$mixdist <- plot.mix(x, size=size, ...)
     ## in verbose mode we output EM fit diagnostics
     if(getOption("RBesT.verbose", FALSE)) {
diff --git a/R/RBesT-package.R b/R/RBesT-package.R
index 1b5964a..9b33b57 100644
--- a/R/RBesT-package.R
+++ b/R/RBesT-package.R
@@ -64,4 +64,5 @@
 #' @import ggplot2
 #' @import Formula
 #' @import checkmate
+#' @import abind
 NULL
diff --git a/R/mix.R b/R/mix.R
index 31fb63f..490e0ab 100644
--- a/R/mix.R
+++ b/R/mix.R
@@ -105,11 +105,13 @@ rmix <- function(mix, n) UseMethod("rmix")
     cl <- grep("mix$", class(mix), ignore.case=TRUE, value=TRUE)
     dl <- dlink(mix)
     if(inherits(mix, "normMix")) s <- sigma(mix)
+    if(inherits(mix, "mvnormMix")) s <- sigma(mix)
     mix <- mix[,...,drop=FALSE]
     if(rescale) mix[1,] <- mix[1,] / sum(mix[1,])
     class(mix) <- cl
     dlink(mix) <- dl
     if(inherits(mix, "normMix")) sigma(mix) <- s
+    if(inherits(mix, "mvnormMix")) sigma(mix) <- s
     mix
 }
 #' @export
@@ -126,7 +128,7 @@ rmix <- function(mix, n) UseMethod("rmix")
 ## IMPLEMENTATION DETAILS
 
 #' @export
-dmix.default <- function(mix, x, log=FALSE) "Unknown mixture"
+dmix.default <- function(mix, x, log=FALSE) stop("Unknown mixture")
 
 ## default implementation which only needs the density function;
 ## assumption is that the first argument of the density corresponds to
@@ -165,9 +167,34 @@ dmix.betaBinomialMix <- function(mix, x, log=FALSE) dmix_impl(Curry(dBetaBinomia
 #' @export
 dmix.gammaPoissonMix <- function(mix, x, log=FALSE) dmix_impl(Curry(.dnbinomAB, n=attr(mix, "n")), mix, x, log)
 
+#' @export
+dmix.mvnormMix <- function(mix, x, log=FALSE) {
+    p <- mvnormdim(mix[-1,1])
+    Nc <- ncol(mix)
+    if(is.matrix(x)) {
+        Nx <- nrow(x)
+        assert_matrix(x, any.missing=FALSE, nrows=Nx, ncols=p)
+    } else if(is.vector(x)) {
+        Nx <- 1
+        assert_vector(x, any.missing=FALSE, len=p)
+    } else {
+        stop("x must a vector or a matrix.")
+    }
+    assert_that(is.mixidentity_link(mix))
+    comp_res <- matrix(NA_real_, nrow=Nx, ncol=Nc)
+    for(i in 1:Nc) {
+        S <- mvnormsigma(mix[-1,i])
+        comp_res[,i] <- log(mix[1,i]) + mvtnorm::dmvnorm(x, mix[2:(p+1), i], sigma=S, log=TRUE, checkSymmetry=FALSE)
+    }
+    res <- matrixStats::rowLogSumExps(comp_res)
+    if(!log)
+        res <- exp(res)
+    return(res)
+}
+
 ## DISTRIBUTION FUNCTIONS
 #' @export
-pmix.default <- function(mix, q, lower.tail = TRUE, log.p=FALSE) "Unknown mixture"
+pmix.default <- function(mix, q, lower.tail = TRUE, log.p=FALSE) stop("Unknown mixture")
 
 pmix_impl <- function(dist, mix, q, lower.tail = TRUE, log.p=FALSE) {
     Nc <- ncol(mix)
@@ -217,10 +244,14 @@ pmix.betaBinomialMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) {
 #' @export
 pmix.gammaPoissonMix <- function(mix, q, lower.tail = TRUE, log.p=FALSE) pmix_impl(Curry(.pnbinomAB, n=attr(mix, "n")), mix, q, lower.tail, log.p)
 
+#' @export
+pmix.mvnormMix <- function(mix, q, ...) stop("Multivariate normal mixture cumulative density not supported.")
+
+
 ## QUANTILE FUNCTION
 
 #' @export
-qmix.default <- function(mix, p, lower.tail = TRUE, log.p=FALSE) "Unknown mixture"
+qmix.default <- function(mix, p, lower.tail = TRUE, log.p=FALSE) stop("Unknown mixture")
 
 qmix_impl <- function(quant, mix, p, lower.tail = TRUE, log.p=FALSE) {
     Nc <- ncol(mix)
@@ -286,11 +317,11 @@ qmix.betaBinomialMix  <- function(mix, p, lower.tail = TRUE, log.p=FALSE) {
 ## internal redefinition of negative binomial
 ##.qnbinomAB <- function(p, a, b, lower.tail = TRUE, log.p = FALSE ) qnbinom(p, size=a, prob=b/(b+1), lower.tail = lower.tail, log.p = log.p )
 .qnbinomAB <- function(p, a, b, n=1, lower.tail = TRUE, log.p = FALSE ) qnbinom(p, size=a, prob=b/(b+n), lower.tail = lower.tail, log.p = log.p )
-#' @export
 ##qmix.gammaPoissonMix <- function(mix, p, lower.tail = TRUE, log.p=FALSE) qmix_impl(Curry(.qnbinomAB, n=attr(mix, "n")), mix, p, lower.tail, log.p, discrete=TRUE)
 
 ## switched to numeric implementation as discretization seems to cause
 ## some trouble in the above definitions
+#' @export
 qmix.gammaPoissonMix <- function(mix, p, lower.tail = TRUE, log.p=FALSE) {
     assert_that(is.dlink_identity(attr(mix, "link")))
     ## numerical evaulation
@@ -308,11 +339,13 @@ qmix.gammaPoissonMix <- function(mix, p, lower.tail = TRUE, log.p=FALSE) {
     ind
 }
 
+#' @export
+qmix.mvnormMix <- function(mix, p, ...) stop("Multivariate normal mixture quantiles not supported.")
 
 ### RANDOM NUMBER GENERATION
 
 #' @export
-rmix.default <- function(mix, n) "Unknown mixture"
+rmix.default <- function(mix, n) stop("Unknown mixture")
 
 rmix_impl <- function(rng, mix, n) {
     ind <-  sample.int(ncol(mix), n, replace = TRUE, prob = mix[1,])
@@ -343,6 +376,29 @@ rmix.betaBinomialMix  <- function(mix, n) {
 #' @export
 rmix.gammaPoissonMix  <- function(mix, n) rmix_impl(Curry(.rnbinomAB, n=attr(mix, "n")),  mix, n)
 
+#' @export
+rmix.mvnormMix <- function(mix, n) {
+    ## sample the mixture components
+    ind <-  sample.int(ncol(mix), n, replace = TRUE, prob = mix["w",])
+    ## sort these
+    sidx <- order(ind)
+    ## ensure we can sort into the original random order
+    oidx <- seq(1, n)[sidx]
+    sind <- ind[sidx]
+    ## now sind[oidx] == ind
+    ## count how many times we need to sample which component
+    r <- rle(sind)
+    p <- mvnormdim(mix[-1,1])
+    samp <- do.call(rbind,
+                    mapply(function(comp, cn) {
+                        m <- mix[2:(p+1),comp]
+                        S <- mvnormsigma(mix[-1,comp])
+                        rmvnorm(cn, m, S, checkSymmetry=FALSE)
+                    },
+                    r$values, r$lengths, SIMPLIFY=FALSE))[oidx,,drop=FALSE]
+    attr(samp, "ind") <- ind
+    samp
+}
 
 #' @export
 print.mix <- function(x, digits, ...) {
diff --git a/R/mixbeta.R b/R/mixbeta.R
index d609e2a..b7ffd9f 100644
--- a/R/mixbeta.R
+++ b/R/mixbeta.R
@@ -56,6 +56,7 @@ NULL
 #' @export
 mixbeta <- function(..., param=c("ab", "ms", "mn")) {
     mix <- mixdist3(...)
+    assert_matrix(mix, nrows=3, any.missing=FALSE)
     param <- match.arg(param)
     mix[c(2,3),] <- switch(param,
                            ab=mix[c(2,3),],
diff --git a/R/mixcombine.R b/R/mixcombine.R
index 996a35b..43033c9 100644
--- a/R/mixcombine.R
+++ b/R/mixcombine.R
@@ -47,5 +47,6 @@ mixcombine <- function(..., weight, rescale=TRUE) {
     dlink(mix) <- dl
     likelihood(mix) <- lik
     if("normMix" %in% cl) sigma(mix) <- sigma(comp[[1]])
+    if("mvnormMix" %in% cl) sigma(mix) <- sigma(comp[[1]])
     mix
 }
diff --git a/R/mixdiff.R b/R/mixdiff.R
index cebab4f..ccdff48 100644
--- a/R/mixdiff.R
+++ b/R/mixdiff.R
@@ -90,6 +90,7 @@ mixnormdiff <- function(mix1, mix2) {
 #' @rdname mixdiff
 #' @export
 dmixdiff <- function(mix1, mix2, x) {
+    assert_that(!inherits(mix1, "mvnormMix") & !inherits(mix2, "mvnormMix"), msg="Multivariate normal mixture density not supported.")
     if(inherits(mix1, "normMix") & inherits(mix2, "normMix"))
         return(dmix(mixnormdiff(mix1, mix2), x=x))
     interval <- support(mix1)
@@ -123,6 +124,7 @@ dmixdiff <- function(mix1, mix2, x) {
 #' @rdname mixdiff
 #' @export
 pmixdiff <- function(mix1, mix2, q, lower.tail=TRUE) {
+    assert_that(!inherits(mix1, "mvnormMix") & !inherits(mix2, "mvnormMix"), msg="Multivariate normal mixture density not supported.")
     if(inherits(mix1, "normMix") & inherits(mix2, "normMix"))
         return(pmix(mixnormdiff(mix1, mix2), q=q, lower.tail=lower.tail))
     interval <- support(mix1)
@@ -154,6 +156,7 @@ pmixdiff <- function(mix1, mix2, q, lower.tail=TRUE) {
 #' @rdname mixdiff
 #' @export
 qmixdiff <- function(mix1, mix2, p, lower.tail = TRUE) {
+    assert_that(!inherits(mix1, "mvnormMix") & !inherits(mix2, "mvnormMix"), msg="Multivariate normal mixture density not supported.")
     if(inherits(mix1, "normMix"))
         return(qmix(mixnormdiff(mix1, mix2), p=p, lower.tail=lower.tail))
     interval <- support(mix1)
diff --git a/R/mixdist3.R b/R/mixdist3.R
index 0c10bdc..dd7cda9 100644
--- a/R/mixdist3.R
+++ b/R/mixdist3.R
@@ -4,8 +4,9 @@
 mixdist3 <- function(...) {
     args <- list(...)
     Nc <- length(args)
-    if(!all(sapply(args, length) == 3))
-        stop("All components must have 3 parameters.")
+    l <- sapply(args, length)
+    if(!all(l >= 3) | !all(l == l[1]))
+        stop("All components must have equal number of parameters.")
     res <- do.call(cbind, args)
     if(is.null(names(args)))
         colnames(res) <- paste("comp", seq(Nc), sep="")
diff --git a/R/mixfit.R b/R/mixfit.R
index 9f73b64..b949b7b 100644
--- a/R/mixfit.R
+++ b/R/mixfit.R
@@ -81,19 +81,18 @@
 #'
 
 #' @export
-mixfit <- function(sample, type=c("norm", "beta", "gamma"), thin, ...) UseMethod("mixfit")
+mixfit <- function(sample, type=c("norm", "beta", "gamma", "mvnorm"), thin, ...) UseMethod("mixfit")
 
 #' @describeIn mixfit Performs an EM fit for the given
 #' sample. Thinning is applied only if thin is specified.
 #' @export
-mixfit.default <- function(sample, type=c("norm", "beta", "gamma"), thin, ...)
-{
+mixfit.default <- function(sample, type=c("norm", "beta", "gamma", "mvnorm"), thin, ...) {
     type <- match.arg(type)
-    assert_that(type %in% c("norm", "beta", "gamma"))
-    EM <- switch(type, norm=EM_nmm, beta=EM_bmm_ab, gamma=EM_gmm)
+    assert_that(type %in% c("norm", "beta", "gamma", "mvnorm"))
+    EM <- switch(type, norm=EM_nmm, beta=EM_bmm_ab, gamma=EM_gmm, mvnorm=EM_mnmm)
     if(!missing(thin)) {
         assert_that(thin >= 1)
-        sample <- sample[seq(1,NROW(sample),by=thin),drop=FALSE]
+        sample <- asub(sample, seq(1,NROW(sample),by=thin), dims=1, drop=FALSE)
     }
     EM(sample, ...)
 }
@@ -107,8 +106,7 @@ mixfit.default <- function(sample, type=c("norm", "beta", "gamma"), thin, ...)
 #' normal mixture will set the reference scale to the estimated
 #' sigma in \code{\link{gMAP}} call.
 #' @export
-mixfit.gMAP <- function(sample, type, thin, ...)
-{
+mixfit.gMAP <- function(sample, type, thin, ...) {
     family <- sample$family$family
     ## automatically thin sample as estimated by gMAP function
     if(missing(thin)) {
@@ -129,8 +127,7 @@ mixfit.gMAP <- function(sample, type, thin, ...)
 #' @describeIn mixfit Fits a mixture density for each prediction from
 #' the \code{\link{gMAP}} prediction.
 #' @export
-mixfit.gMAPpred <- function(sample, type, thin, ...)
-{
+mixfit.gMAPpred <- function(sample, type, thin, ...) {
     if(attr(sample, "type") == "response") {
         type <- switch(attr(sample, "family")$family, binomial = "beta", gaussian = "norm", poisson = "gamma", "unknown")
         family <- attr(sample, "family")$family
@@ -157,10 +154,8 @@ mixfit.gMAPpred <- function(sample, type, thin, ...)
 #' expected to have 3 dimensions which are nested as iterations,
 #' chains, and draws.
 #' @export
-mixfit.array <- function(sample, type, thin, ...)
-{
-    Nmcmc <- prod(dim(sample))
-    if(dim(sample)[3] != 1)
+mixfit.array <- function(sample, type, thin, ...) {
+    if(type != "mvnorm" & dim(sample)[3] != 1)
         stop("Only univariate data is supported.")
     mixfit.default(sample, type, thin, ...)
 }
diff --git a/R/mixgamma.R b/R/mixgamma.R
index 611441c..750ee8c 100644
--- a/R/mixgamma.R
+++ b/R/mixgamma.R
@@ -70,6 +70,7 @@ NULL
 #' @export
 mixgamma <- function(..., param=c("ab", "ms", "mn"), likelihood=c("poisson", "exp")) {
     mix <- mixdist3(...)
+    assert_matrix(mix, nrows=3, any.missing=FALSE)
     param <- match.arg(param)
     likelihood <- match.arg(likelihood)
     mix[c(2,3),] <- switch(param,
diff --git a/R/mixmvnorm.R b/R/mixmvnorm.R
new file mode 100644
index 0000000..35873c4
--- /dev/null
+++ b/R/mixmvnorm.R
@@ -0,0 +1,193 @@
+#' @name mixmvnorm
+#'
+#' @title Multivariate Normal Mixture Density
+#'
+#' @description The multivariate normal mixture density and auxiliary
+#'     functions.
+#'
+#' @param ... List of mixture components.
+#' @param sigma Reference covariance.
+#' @param param Determines how the parameters in the list are
+#'     interpreted. See details.
+#' @param object Multivariate normal mixture object.
+#'
+#' @details Each entry in the \code{...} argument list is a numeric
+#'     vector defining one component of the mixture multivariate
+#'     normal distribution. The first entry of the component defining
+#'     vector is the weight of the mixture component followed by the
+#'     vector of means in each dimension and finally a specification
+#'     of the covariance matrix, which depends on the chosen
+#'     parametrization. The covariance matrix is expected to be given
+#'     as numeric vector in a column-major format, which is standard
+#'     conversion applied to matrices by the vector concatenation
+#'     function \code{\link[base:c]{c}}. Please refer to the examples
+#'     section below.
+#'
+#' Each component defining vector can be specified in different ways
+#' as determined by the \code{param} option:
+#' 
+#' \describe{
+#' \item{ms}{Mean vector and covariance matrix \code{s}. Default.}
+#' \item{mn}{Mean vector and number of observations. \code{n} determines the covariance for each component via the relation \eqn{\Sigma/n} with \eqn{\Sigma} being the known reference covariance.}
+#' }
+#'
+#' The reference covariance \eqn{\Sigma} is the known covariance in
+#' the normal-normal model (observation covariance). The function
+#' \code{sigma} can be used to query the reference covariance and may
+#' also be used to assign a new reference covariance, see examples
+#' below. In case \code{sigma} is not specified, the user has to
+#' supply \code{sigma} as argument to functions which require a
+#' reference covariance.
+#'
+#' @family mixdist
+#'
+#' @return Returns a multivariate normal mixture with the specified
+#'     mixture components.
+#'
+#' @examples
+#'
+#' S <- diag(c(1, 2)) %*% matrix(c(1, 0.5, 0.5, 1), 2, 2) %*% diag(c(1, 2)) 
+#' mvnm1 <- mixmvnorm(rob=c(0.2, c(0, 0), diag(c(5, 5))),
+#'                    inf=c(0.8, c(0.5, 1), S/10), sigma=S)
+#'
+#' print(mvnm1)
+#' summary(mvnm1)
+#'
+#' set.seed(657846)
+#' mixSamp1 <- rmix(mvnm1, 500)
+#' colMeans(mixSamp1)
+#'
+#'
+NULL
+
+#' @rdname mixmvnorm
+#' @export
+mixmvnorm <- function(..., sigma, param=c("ms", "mn")) {
+    ## length of first mean vector determines dimension
+    mix <- mixdist3(...)
+    param <- match.arg(param)
+    Nc <- ncol(mix)
+    n <- colnames(mix)
+    if(param == "ms") {
+        ## mean vector & covariance parametrization
+        l <- nrow(mix)
+        p <- (sqrt(1 + 4 * (l - 1)) - 1) / 2
+        assert_integerish(p, lower=1, any.missing=FALSE, len=1)
+        ## in this case we expect c(weight, mean, as.numeric(cov))
+        mix <- do.call(mixdist3,
+                       lapply(1:Nc,
+                              function(co) c(mix[1,co], mvnorm(mix[2:(p+1),co], matrix(mix[(1+p+1):(1+p+p^2),co], p, p)))
+                              ))
+    }
+    if(param == "mn") {
+        ## mean vector & number of observations
+        l <- nrow(mix)
+        p <- l - 2
+        assert_integerish(p, lower=1, any.missing=FALSE, len=1)
+        assert_matrix(sigma, any.missing=FALSE, nrows=p, ncols=p)
+        mix <- do.call(mixdist3, lapply(1:Nc,
+                                        function(co) {
+                                            assert_numeric(mix[l,co], lower=0, finite=TRUE, any.missing=FALSE)
+                                            c(mix[1,co], mvnorm(mix[2:(p+1),co], sigma/mix[l,co] ))
+                                        }
+                                        ))
+    }
+    colnames(mix) <- n
+    p <- mvnormdim(mix[-1,1])
+    rownames(mix) <- c("w", mvnorm_label(mix[-1,1]))
+    if(!missing(sigma)) {
+        assert_matrix(sigma, any.missing=FALSE, nrows=p, ncols=p)
+        attr(mix, "sigma") <- sigma
+    }
+    class(mix) <- c("mvnormMix", "mix")
+    likelihood(mix) <- "mvnormal"
+    mix
+}
+
+#' @keywords internal
+mvnorm <- function(mean, sigma) {
+    ## TODO: far more checks!!! Allow to pass in directly a cholesky factor
+    assert_numeric(mean, finite=TRUE, any.missing=FALSE)
+    p <- length(mean)
+    assert_matrix(sigma, any.missing=FALSE, nrows=p, ncols=p)
+    rho <- cov2cor(sigma)
+    s <- sqrt(diag(sigma))
+    mvn <- c(mean, s, rho[lower.tri(rho)])
+    names(mvn) <- mvnorm_label(mvn)
+    mvn
+}
+
+#' @keywords internal
+mvnormdim <- function(mvn) {
+    p <- (-3 + sqrt(9 + 8 * length(mvn))) / 2
+    assert_integerish(p, lower=1, any.missing=FALSE, len=1)
+    as.integer(p)
+}
+
+#' @keywords internal
+mvnorm_label <- function(mvn) {
+    p <- mvnormdim(mvn)
+    if(p > 1) {
+        Rho_labs_idx <- outer(1:p, 1:p, paste, sep=",")[lower.tri(diag(p))]
+        lab <- c(paste0("m[", 1:p,"]"), paste0("s[", 1:p, "]"), paste0("rho[", Rho_labs_idx, "]"))
+    } else {
+        lab <- c(paste0("m[", 1:p,"]"), paste0("s[", 1:p, "]"))
+    }
+    lab
+}
+
+#' @keywords internal
+mvnormsigma <- function(mvn) {
+    p <- mvnormdim(mvn)
+    n <- length(mvn)
+    s <- mvn[(p+1):(2*p)]
+    Rho <- diag(p)
+    Rho[lower.tri(Rho)] <- mvn[(2*p+1):n]
+    Rho[upper.tri(Rho)] <- t(Rho)[upper.tri(Rho)]
+    diag(s, nrow=p) %*% Rho %*% diag(s, nrow=p)
+}
+
+
+#' @rdname mixmvnorm
+#' @method print mvnormMix
+#' @param x The mixture to print
+#' @export
+print.mvnormMix <- function(x, ...) {
+    cat("Multivariate normal mixture\n")
+    cat(paste0("Outcome dimension: ", mvnormdim(x[-1,1]), "\n"))
+    if(!is.null(sigma(x))) {
+        cat("Reference covariance:\n")
+        print(sigma(x), ...)
+    }
+    NextMethod()
+}
+
+#' @rdname mixmvnorm
+#' @method summary mvnormMix
+#' @export
+summary.mvnormMix <- function(object, ...) {
+    w <- object[1,]
+    p <- mvnormdim(object[-1,1])
+    m <- object[2:(p+1),,drop=FALSE]
+    Nc <- ncol(object)
+    mmix <- rowSums(sweep(m, 2, w, "*"))
+    ## Cov(x,x) = E(x x') - E(x) E(x')
+    ## E(x x') = Sigma + m m' (see matrix cookbook eq 377)
+    if(Nc == 1) {
+        S <- w[1] * mvnormsigma(object[-1,1])
+    } else {
+        S <- -1 * tcrossprod(mmix)
+        for(i in 1:Nc) {
+            S <- S + w[i] * ( mvnormsigma(object[-1,i]) + tcrossprod(unname(m[,i,drop=FALSE]) ))
+        }
+    }
+    list(mean=mmix, cov=S)
+}
+
+#' @rdname mixmvnorm
+#' @method sigma mvnormMix
+#' @export
+#' @export sigma
+sigma.mvnormMix <- function(object, ...) {
+    attr(object, "sigma")
+}
diff --git a/R/mixnorm.R b/R/mixnorm.R
index 66c3ab6..f9dd7b8 100644
--- a/R/mixnorm.R
+++ b/R/mixnorm.R
@@ -50,7 +50,7 @@
 #' summary(nm)
 #' plot(nm)
 #'
-#' set.seed(1)
+#' set.seed(57845)
 #' mixSamp <- rmix(nm, 500)
 #' plot(nm, samp=mixSamp)
 #'
@@ -74,6 +74,7 @@ NULL
 #' @export
 mixnorm <- function(..., sigma, param=c("ms", "mn")) {
     mix <- mixdist3(...)
+    assert_matrix(mix, nrows=3, any.missing=FALSE)
     param <- match.arg(param)
     mix[c(2,3),] <- switch(param,
                            ms=mix[c(2,3),],
diff --git a/R/mixplot.R b/R/mixplot.R
index 314c3fd..5766e0c 100644
--- a/R/mixplot.R
+++ b/R/mixplot.R
@@ -1,4 +1,8 @@
-#' Plot mixture distributions
+#' @name mixplot
+#' 
+#' @title Plot mixture distributions
+#'
+#' @description Plotting for mixture distributions
 #'
 #' @param x mixture distribution
 #' @param prob  defining lower and upper percentile of x-axis. Defaults to the 99\% central probability mass.
@@ -55,6 +59,7 @@
 #' pl <- plot(nm)
 #' pl + ggtitle("Normal 2-Component Mixture")
 #'
+#' @rdname mixplot
 #' @method plot mix
 #' @export
 plot.mix <- function(x, prob=0.99, fun=dmix, log=FALSE, comp=TRUE, size=1.25, ...) {
@@ -106,3 +111,10 @@ plot.mix <- function(x, prob=0.99, fun=dmix, log=FALSE, comp=TRUE, size=1.25, ..
     }
     pl
 }
+
+#' @rdname mixplot
+#' @method plot mvnormMix
+#' @export
+plot.mvnormMix <- function(x, prob=0.99, fun=dmix, log=FALSE, comp=TRUE, size=1.25, ...) {
+    stop("Multivariate normal mixture plotting not supported.")
+}
diff --git a/R/preddist.R b/R/preddist.R
index 71c8bad..34c6eb0 100644
--- a/R/preddist.R
+++ b/R/preddist.R
@@ -78,10 +78,10 @@
 #' @export
 preddist <- function(mix, ...) UseMethod("preddist")
 #' @export
-preddist.default <- function(mix, ...) "Unknown distribution"
+preddist.default <- function(mix, ...) stop("Unknown distribution")
 
 #' @describeIn preddist Obtain the matching predictive distribution
-#' for a beta distribution, the BetaBinomial.
+#'     for a beta distribution, the BetaBinomial.
 #' @export
 preddist.betaMix <- function(mix, n=1, ...) {
     attr(mix, "n") <- n
@@ -90,11 +90,12 @@ preddist.betaMix <- function(mix, n=1, ...) {
 }
 
 #' @describeIn preddist Obtain the matching predictive distribution
-#' for a Normal with constant standard deviation. Note that the
-#' reference scale of the returned Normal mixture is meaningless as the
-#' individual components are updated appropriatley.
+#'     for a Normal with constant standard deviation. Note that the
+#'     reference scale of the returned Normal mixture is meaningless
+#'     as the individual components are updated appropriatley.
 #' @param sigma The fixed reference scale of a normal mixture. If left
-#' unspecified, the default reference scale of the mixture is assumed.
+#'     unspecified, the default reference scale of the mixture is
+#'     assumed.
 #' @export
 preddist.normMix <- function(mix, n=1, sigma, ...) {
     if(missing(sigma)) {
@@ -116,7 +117,7 @@ preddist.normMix <- function(mix, n=1, sigma, ...) {
 }
 
 #' @describeIn preddist Obtain the matching predictive distribution
-#' for a Gamma. Only Poisson likelihoods are supported.
+#'     for a Gamma. Only Poisson likelihoods are supported.
 #' @export
 preddist.gammaMix <- function(mix, n=1, ...) {
     assert_set_equal(likelihood(mix), "poisson")
@@ -127,3 +128,7 @@ preddist.gammaMix <- function(mix, n=1, ...) {
     mix
 }
 
+#' @describeIn preddist Multivariate normal mixtures predictive
+#'     densities are not (yet) supported.
+#' @export
+preddist.mvnormMix <- function(mix, ...) stop("Multivariate normal mixture predictive density not supported.")
diff --git a/R/support.R b/R/support.R
index 7e8b8b2..43ab9d2 100644
--- a/R/support.R
+++ b/R/support.R
@@ -7,10 +7,12 @@
 #' @keywords internal
 support <- function(mix) UseMethod("support")
 #' @export
-support.default <- function(mix) "Unknown mixture"
+support.default <- function(mix) stop("Unknown mixture")
 #' @export
 support.betaMix <- function(mix) mixlink(mix, c(0,1))
 #' @export
 support.gammaMix <- function(mix) mixlink(mix, c(0,Inf))
 #' @export
 support.normMix <- function(mix) mixlink(mix, c(-Inf,Inf))
+#' @export
+support.mvnormMix <- function(mix) matrix(c(-Inf, Inf), nrow=mvnormdim(mix[-1,1]), ncol=2)
diff --git a/R/sysdata.rda b/R/sysdata.rda
index bf74a39..e7f2104 100644
Binary files a/R/sysdata.rda and b/R/sysdata.rda differ
diff --git a/man/mix.Rd b/man/mix.Rd
index f2d5356..8a65722 100644
--- a/man/mix.Rd
+++ b/man/mix.Rd
@@ -123,7 +123,8 @@ Other mixdist:
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
+\code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{plot.mix}()}
+\code{\link{mixplot}}
 }
 \concept{mixdist}
diff --git a/man/mixbeta.Rd b/man/mixbeta.Rd
index a37db11..8f2f6a5 100644
--- a/man/mixbeta.Rd
+++ b/man/mixbeta.Rd
@@ -88,8 +88,9 @@ print(bm4)
 Other mixdist: 
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
+\code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{mix}},
-\code{\link{plot.mix}()}
+\code{\link{mixplot}},
+\code{\link{mix}}
 }
 \concept{mixdist}
diff --git a/man/mixcombine.Rd b/man/mixcombine.Rd
index c7dd900..4925797 100644
--- a/man/mixcombine.Rd
+++ b/man/mixcombine.Rd
@@ -43,8 +43,9 @@ mixcombine(bm, unif, weight=c(9, 1))
 Other mixdist: 
 \code{\link{mixbeta}()},
 \code{\link{mixgamma}()},
+\code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{mix}},
-\code{\link{plot.mix}()}
+\code{\link{mixplot}},
+\code{\link{mix}}
 }
 \concept{mixdist}
diff --git a/man/mixfit.Rd b/man/mixfit.Rd
index 0a418e5..475497c 100644
--- a/man/mixfit.Rd
+++ b/man/mixfit.Rd
@@ -8,9 +8,9 @@
 \alias{mixfit.array}
 \title{Fit of Mixture Densities to Samples}
 \usage{
-mixfit(sample, type = c("norm", "beta", "gamma"), thin, ...)
+mixfit(sample, type = c("norm", "beta", "gamma", "mvnorm"), thin, ...)
 
-\method{mixfit}{default}(sample, type = c("norm", "beta", "gamma"), thin, ...)
+\method{mixfit}{default}(sample, type = c("norm", "beta", "gamma", "mvnorm"), thin, ...)
 
 \method{mixfit}{gMAP}(sample, type, thin, ...)
 
diff --git a/man/mixgamma.Rd b/man/mixgamma.Rd
index 33627f2..e1d97cb 100644
--- a/man/mixgamma.Rd
+++ b/man/mixgamma.Rd
@@ -106,8 +106,9 @@ gfmix <- mixgamma(rob1=c(0.15, mn2gamma(2, 1)), rob2=c(0.15, ms2gamma(2, 5)), in
 Other mixdist: 
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
+\code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
-\code{\link{mix}},
-\code{\link{plot.mix}()}
+\code{\link{mixplot}},
+\code{\link{mix}}
 }
 \concept{mixdist}
diff --git a/man/mixmvnorm.Rd b/man/mixmvnorm.Rd
new file mode 100644
index 0000000..3c352d1
--- /dev/null
+++ b/man/mixmvnorm.Rd
@@ -0,0 +1,91 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mixmvnorm.R
+\name{mixmvnorm}
+\alias{mixmvnorm}
+\alias{print.mvnormMix}
+\alias{summary.mvnormMix}
+\alias{sigma.mvnormMix}
+\title{Multivariate Normal Mixture Density}
+\usage{
+mixmvnorm(..., sigma, param = c("ms", "mn"))
+
+\method{print}{mvnormMix}(x, ...)
+
+\method{summary}{mvnormMix}(object, ...)
+
+\method{sigma}{mvnormMix}(object, ...)
+}
+\arguments{
+\item{...}{List of mixture components.}
+
+\item{sigma}{Reference covariance.}
+
+\item{param}{Determines how the parameters in the list are
+interpreted. See details.}
+
+\item{x}{The mixture to print}
+
+\item{object}{Multivariate normal mixture object.}
+}
+\value{
+Returns a multivariate normal mixture with the specified
+    mixture components.
+}
+\description{
+The multivariate normal mixture density and auxiliary
+    functions.
+}
+\details{
+Each entry in the \code{...} argument list is a numeric
+    vector defining one component of the mixture multivariate
+    normal distribution. The first entry of the component defining
+    vector is the weight of the mixture component followed by the
+    vector of means in each dimension and finally a specification
+    of the covariance matrix, which depends on the chosen
+    parametrization. The covariance matrix is expected to be given
+    as numeric vector in a column-major format, which is standard
+    conversion applied to matrices by the vector concatenation
+    function \code{\link[base:c]{c}}. Please refer to the examples
+    section below.
+
+Each component defining vector can be specified in different ways
+as determined by the \code{param} option:
+
+\describe{
+\item{ms}{Mean vector and covariance matrix \code{s}. Default.}
+\item{mn}{Mean vector and number of observations. \code{n} determines the covariance for each component via the relation \eqn{\Sigma/n} with \eqn{\Sigma} being the known reference covariance.}
+}
+
+The reference covariance \eqn{\Sigma} is the known covariance in
+the normal-normal model (observation covariance). The function
+\code{sigma} can be used to query the reference covariance and may
+also be used to assign a new reference covariance, see examples
+below. In case \code{sigma} is not specified, the user has to
+supply \code{sigma} as argument to functions which require a
+reference covariance.
+}
+\examples{
+
+S <- diag(c(1, 2)) \%*\% matrix(c(1, 0.5, 0.5, 1), 2, 2) \%*\% diag(c(1, 2)) 
+mvnm1 <- mixmvnorm(rob=c(0.2, c(0, 0), diag(c(5, 5))),
+                   inf=c(0.8, c(0.5, 1), S/10), sigma=S)
+
+print(mvnm1)
+summary(mvnm1)
+
+set.seed(657846)
+mixSamp1 <- rmix(mvnm1, 500)
+colMeans(mixSamp1)
+
+
+}
+\seealso{
+Other mixdist: 
+\code{\link{mixbeta}()},
+\code{\link{mixcombine}()},
+\code{\link{mixgamma}()},
+\code{\link{mixnorm}()},
+\code{\link{mixplot}},
+\code{\link{mix}}
+}
+\concept{mixdist}
diff --git a/man/mixnorm.Rd b/man/mixnorm.Rd
index 596aa7f..cac7eb3 100644
--- a/man/mixnorm.Rd
+++ b/man/mixnorm.Rd
@@ -86,7 +86,7 @@ print(nm)
 summary(nm)
 plot(nm)
 
-set.seed(1)
+set.seed(57845)
 mixSamp <- rmix(nm, 500)
 plot(nm, samp=mixSamp)
 
@@ -110,7 +110,8 @@ Other mixdist:
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
-\code{\link{mix}},
-\code{\link{plot.mix}()}
+\code{\link{mixmvnorm}()},
+\code{\link{mixplot}},
+\code{\link{mix}}
 }
 \concept{mixdist}
diff --git a/man/plot.mix.Rd b/man/mixplot.Rd
similarity index 93%
rename from man/plot.mix.Rd
rename to man/mixplot.Rd
index f5f5822..7fe703c 100644
--- a/man/plot.mix.Rd
+++ b/man/mixplot.Rd
@@ -1,10 +1,14 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/mixplot.R
-\name{plot.mix}
+\name{mixplot}
+\alias{mixplot}
 \alias{plot.mix}
+\alias{plot.mvnormMix}
 \title{Plot mixture distributions}
 \usage{
 \method{plot}{mix}(x, prob = 0.99, fun = dmix, log = FALSE, comp = TRUE, size = 1.25, ...)
+
+\method{plot}{mvnormMix}(x, prob = 0.99, fun = dmix, log = FALSE, comp = TRUE, size = 1.25, ...)
 }
 \arguments{
 \item{x}{mixture distribution}
@@ -27,7 +31,7 @@ density in addition to the density.}
 A \code{\link[ggplot2]{ggplot}} object is returned.
 }
 \description{
-Plot mixture distributions
+Plotting for mixture distributions
 }
 \details{
 Plot function for mixture distribution objects. It shows
@@ -91,6 +95,7 @@ Other mixdist:
 \code{\link{mixbeta}()},
 \code{\link{mixcombine}()},
 \code{\link{mixgamma}()},
+\code{\link{mixmvnorm}()},
 \code{\link{mixnorm}()},
 \code{\link{mix}}
 }
diff --git a/man/preddist.Rd b/man/preddist.Rd
index d12d2ed..ac2c6e7 100644
--- a/man/preddist.Rd
+++ b/man/preddist.Rd
@@ -5,6 +5,7 @@
 \alias{preddist.betaMix}
 \alias{preddist.normMix}
 \alias{preddist.gammaMix}
+\alias{preddist.mvnormMix}
 \title{Predictive Distributions for Mixture Distributions}
 \usage{
 preddist(mix, ...)
@@ -14,6 +15,8 @@ preddist(mix, ...)
 \method{preddist}{normMix}(mix, n = 1, sigma, ...)
 
 \method{preddist}{gammaMix}(mix, n = 1, ...)
+
+\method{preddist}{mvnormMix}(mix, ...)
 }
 \arguments{
 \item{mix}{mixture distribution}
@@ -23,7 +26,8 @@ preddist(mix, ...)
 \item{n}{predictive sample size, set by default to 1}
 
 \item{sigma}{The fixed reference scale of a normal mixture. If left
-unspecified, the default reference scale of the mixture is assumed.}
+unspecified, the default reference scale of the mixture is
+assumed.}
 }
 \value{
 The function returns for a normal, beta or gamma mixture
@@ -59,12 +63,15 @@ for a beta distribution, the BetaBinomial.
 
 \item \code{preddist(normMix)}: Obtain the matching predictive distribution
 for a Normal with constant standard deviation. Note that the
-reference scale of the returned Normal mixture is meaningless as the
-individual components are updated appropriatley.
+reference scale of the returned Normal mixture is meaningless
+as the individual components are updated appropriatley.
 
 \item \code{preddist(gammaMix)}: Obtain the matching predictive distribution
 for a Gamma. Only Poisson likelihoods are supported.
 
+\item \code{preddist(mvnormMix)}: Multivariate normal mixtures predictive
+densities are not (yet) supported.
+
 }}
 \section{Supported Conjugate Prior-Likelihood Pairs}{
 
diff --git a/tests/testthat/test-EM.R b/tests/testthat/test-EM.R
index f3d3c1d..58ff507 100644
--- a/tests/testthat/test-EM.R
+++ b/tests/testthat/test-EM.R
@@ -41,6 +41,32 @@ ref$norm_bi <- mixnorm(c(0.5, 0, 0.5),
                        c(0.25, -1, 2),
                        param="mn", sigma=1)
 
+p <- 4
+Rho <- diag(p)
+Rho[lower.tri(Rho)] <- c(0.3, 0.8, -0.2, 0.1, 0.5, -0.4)
+Rho[upper.tri(Rho)] <- t(Rho)[upper.tri(Rho)]
+s <- c(1, 2, 3, 4)
+S <- diag(s, p) %*% Rho %*% diag(s, p)
+zero <- rep(0, p)
+
+ref$mvnorm_single <- mixmvnorm(c(1, zero, 5),
+                               param="mn", sigma=S)
+
+ref$mvnorm_heavy <- mixmvnorm(c(0.5, zero, 0.25),
+                              c(0.5, zero + 1, 5),
+                              param="mn", sigma=S)
+
+ref$mvnorm_bi <- mixmvnorm(c(0.5, zero, 0.5),
+                           c(0.25, zero + 1, 5),
+                           c(0.25, zero - 1, 2),
+                           param="mn", sigma=S)
+
+ref$mvnorm_bi_1D <- mixmvnorm(c(0.5, 0, 0.5),
+                              c(0.25, 1, 5),
+                              c(0.25, -1, 2),
+                              param="mn", sigma=S[1,1,drop=FALSE])
+
+
 ref$beta_single <- mixbeta(c(1, 0.3, 10),
                            param="mn")
 
@@ -76,7 +102,7 @@ EM_test <- function(mixTest, seed, Nsim=1e4, verbose=FALSE, ...) {
     set.seed(seed)
     samp <- rmix(mixTest, Nsim)
     EMmix <- mixfit(samp,
-                    type=switch(class(mixTest)[1], gammaMix="gamma", normMix="norm", betaMix="beta"),
+                    type=switch(class(mixTest)[1], gammaMix="gamma", normMix="norm", betaMix="beta", mvnormMix="mvnorm"),
                     thin=1,
                     eps=2,
                     Nc=ncol(mixTest),
@@ -85,10 +111,28 @@ EM_test <- function(mixTest, seed, Nsim=1e4, verbose=FALSE, ...) {
     expect_true(kl < KLthresh)
 }
 
+EM_mvn_test <- function(mixTest, seed, Nsim=1e4, verbose=FALSE, ...) {
+    set.seed(seed)
+    samp <- rmix(mixTest, Nsim)
+    EMmix <- mixfit(samp,
+                    type="mvnorm",
+                    thin=1,
+                    eps=2,
+                    Nc=ncol(mixTest),
+                    verbose=verbose, ...)
+    expect_equal(summary(mixTest)$mean, summary(EMmix)$mean, tolerance=0.1)
+    expect_equal(summary(mixTest)$cov, summary(EMmix)$cov, tolerance=0.1)
+}
+
 test_that("Normal EM fits single component",     EM_test(ref$norm_single, 3453563, Nsim, verbose))
 test_that("Normal EM fits heavy-tailed mixture", EM_test(ref$norm_heavy,  9275624, Nsim, verbose))
 test_that("Normal EM fits bi-modal mixture",     EM_test(ref$norm_bi,     9345726, Nsim, verbose))
 
+test_that("Multivariate Normal EM fits single component",     EM_mvn_test(ref$mvnorm_single, 3453563, Nsim, verbose))
+test_that("Multivariate Normal EM fits heavy-tailed mixture", EM_mvn_test(ref$mvnorm_heavy,  9275624, Nsim, verbose))
+test_that("Multivariate Normal EM fits bi-modal mixture",     EM_mvn_test(ref$mvnorm_bi,     9345726, Nsim, verbose))
+test_that("Multivariate Normal EM fits bi-modal mixture 1D",     EM_mvn_test(ref$mvnorm_bi_1D,     9345726, Nsim, verbose))
+
 test_that("Gamma EM fits single component",      EM_test(ref$gamma_single, 9345835, Nsim, verbose))
 test_that("Gamma EM fits heavy-tailed mixture",  EM_test(ref$gamma_heavy,  5629389, Nsim, verbose))
 test_that("Gamma EM fits bi-modal mixture",      EM_test(ref$gamma_bi,     9373515, Nsim, verbose))
diff --git a/tests/testthat/test-mixdist.R b/tests/testthat/test-mixdist.R
index 4f0d169..e552795 100644
--- a/tests/testthat/test-mixdist.R
+++ b/tests/testthat/test-mixdist.R
@@ -22,8 +22,6 @@ betaMix <- mixbeta(c(0.8, 11, 4), c(0.2, 1, 1))
 gamma <- mixgamma(c(1, 5, 10), param="mn")
 gammaMix <- mixgamma(rob=c(0.25, 8, 0.5), inf=c(0.75, 8, 10), param="mn")
 
-nm <- mixnorm(rob=c(0.2, 0, 2), inf=c(0.8, 2, 2), sigma=5)
-
 norm <- mixnorm(c(1, 0, sqrt(2)), sigma=1)
 
 normMix <- mixnorm(c(0.2, 0, 2), c(0.8, 2, 2), sigma=1)
@@ -108,3 +106,57 @@ test_that("Gamma mixture CDF function is consistent", consistent_cdf(gammaMix, s
 ## problematic Beta density
 bm <- mixbeta(c(1.0, 298.30333970, 146.75306521))
 test_that("Problematic (1) BetaBinomial CDF function is consistent", consistent_cdf(preddist(bm, n=50), 0:50))
+
+## tests for the multivariate normal mixture density
+p <- 4
+Rho <- diag(p)
+Rho[lower.tri(Rho)] <- c(0.3, 0.8, -0.2, 0.1, 0.5, -0.4)
+Rho[upper.tri(Rho)] <- t(Rho)[upper.tri(Rho)]
+s <- c(1, 2, 3, 4)
+S <- diag(s, p) %*% Rho %*% diag(s, p)
+
+mvn_consistent_dimension <- function(mix, p) {
+    s <- summary(mix)
+    expect_numeric(s$mean, any.missing=FALSE, len=p)
+    expect_matrix(s$cov, any.missing=FALSE, nrows=p, ncols=p)
+}
+
+test_that("Multivariate normal mixture has consistent dimensionality",
+{
+    for(i in 1:(nrow(S)-1)) {
+        p_sub <- 4-i
+        S_sub <- S[-c(1:i), -c(1:i), drop=FALSE]
+        mvn_consistent_dimension(mixmvnorm(c(1, rep(0, p_sub), S_sub), sigma=S_sub), p_sub)
+    }
+})
+
+test_that("Multivariate normal mixture has consistent initialization",
+{
+    p <- nrow(S)
+    mv1 <- mixmvnorm(c(1, rep(0, p), S), sigma=S, param="ms")
+    mv2 <- mixmvnorm(c(1, rep(0, p), 1), sigma=S, param="mn")
+    mv3 <- mixmvnorm(c(1, rep(0, p), 2), sigma=S, param="mn")
+
+    expect_equal(summary(mv1)$cov, S, tolerance=eps_lower)
+    expect_equal(summary(mv2)$cov, S, tolerance=eps_lower)
+    expect_equal(summary(mv3)$cov, S/2, tolerance=eps_lower)
+})
+
+mvn_consistent_summaries <- function(mix, S=Nsamp_equant) {
+    samp <- rmix(mix, S)
+    m <- colMeans(samp)
+    expect_equal(colMeans(samp), summary(mix)$mean, tolerance=eps)
+    expect_equal(unname(cov(samp)), summary(mix)$cov, tolerance=eps)
+}
+
+test_that("Multivariate normal mixture has consistent summaries",
+{
+    p <- nrow(S)
+    mv1 <- mixmvnorm(c(1, rep(0, p), S), sigma=S, param="ms")
+    mv2 <- mixmvnorm(c(1, rep(0, p), 1), sigma=S, param="mn")
+    mv3 <- mixmvnorm(c(0.2, rep(0, p), 2), c(0.8, rep(1, p), 6), sigma=S, param="mn")
+
+    mvn_consistent_summaries(mv1)
+    mvn_consistent_summaries(mv2)
+    mvn_consistent_summaries(mv3)
+})
diff --git a/tools/make-ds.R b/tools/make-ds.R
index 289cbfd..e438293 100644
--- a/tools/make-ds.R
+++ b/tools/make-ds.R
@@ -43,7 +43,7 @@ make_internal_ds <- function() {
     calibration_meta["MD5"]  <- vals["MD5"]
 
     pkg_create_date  <- Sys.time()
-    pkg_sha <- "252a9f8"
+    pkg_sha <- "d5f8521"
 
     if (gsub("\\$", "", pkg_sha) == "Format:%h") {
         pkg_sha <- system("git rev-parse --short HEAD", intern=TRUE)