The multivariate normal and the multivariate
A supposedly simple task in statistics courses and related applications
is to generate random variates from a multivariate
The multivariate normal distribution can be defined in various ways, one is with its stochastic representation
where
In what follows, we assume
Contours of equal density are ellipsoids; all so-called elliptical distributions which admit a density share this property.
A positive definite (semi-definite) matrix
for a lower triangular matrix
The stochastic representation (1), together with the
Cholesky decomposition of
> ## Setup
> mu <- 1:2 # mean vector of X
> Sigma <- matrix(c(4, 2, 2, 3), ncol=2) # covariance matrix of X
> n <- 1000 # sample size
> d <- 2 # dimension
> ## Step 1: Compute the Cholesky factor of Sigma
> L <- t(chol(Sigma)) # t() as chol() returns an upper triangular matrix
> ## Step 2: Generate iid standard normal random variates
> set.seed(271) # set seed for reproducibility
> Z <- matrix(rnorm(n*d), nrow=d, ncol=n) # (d,n)-matrix
> ## Step 3: Reconstruct the stochastic representation
> X <- mu + L %*% Z # (d,n)-matrix of realizations N_d(mu, Sigma)
This idea for sampling
> require(mvtnorm)
> set.seed(271)
> X. <- rmvnorm(n, mean=mu, sigma=Sigma, method="chol") # (n,d)-matrix
> stopifnot(identical(t(X), X.)) # check equality of the random numbers
The default method (method="eigen"
) utilizes the eigendecomposition of
mvrnorm()
of the recommended R package
MASS provides the same
approach; see (Venables and B. D. Ripley 2013) or (Ripley 1987 98). Note however, that
although the internally drawn independent standard normal random
variates are identical, the two algorithms compute different matrices
> require(MASS)
> X.. <- mvrnorm(n, mu=mu, Sigma=Sigma) # (n,d)-matrix
The left-hand side of Figure 1 displays 1000 log-returns of daily stock prices for BMW and Siemens in the decade from 1985-01-02 to 1994-12-30; the data is available in the R package evir (Pfaff 2012). This time period includes various historically “extreme” events such as the stock market crash 1987 known as “Black Monday” (1987-10-19), one of the Monday demonstrations in Leipzig (1989-10-16), and the August Putsch attempt against Mikhail Gorbachev (1991-08-19).
A comparison with simulated data of the same sample size from a fitted
bivariate normal distribution on the right-hand side of
Figure 1 shows that the bivariate normal
distribution is not an adequate model here to account for such (joint)
“extreme” events (in the upper right or lower left corner of the
bivariate distribution); this can also be checked with formal
statistical or graphical goodness-of-fit tests. The bivariate
The multivariate
where
By introducing the additional random factor
As for the multivariate normal distribution, the density (4) has ellipsoidal level sets and thus belongs to the class of elliptical distributions.
As the CRAN Task View “Distributions” reveals, the
R packages mvtnorm and
mnormt (see Azzalini 2012 for the latter) provide functions for drawing random variates from the
multivariate normal and
The left-hand side of Figure 2 shows 2608 simulated
vectors from a bivariate
> require(QRM)
> fit <- fit.mst(X, method = "BFGS") # fit a multivariate t distribution
> mu <- fit$mu # estimated location vector
> Sigma <- as.matrix(fit$Sigma) # estimated scale matrix
> nu <- fit$df # estimated degrees of freedom
In comparison to the sample from the multivariate normal distribution on
the right-hand side of Figure 1, the multivariate
We now turn to the task of generating vectors of random variates from a
multivariate
and try to generate
The obvious generalization of drawing rmvnorm(n, mean=mu, sigma=Sigma)
(with n
being mu
being
Sigma
being rmvt(n, mean=mu, sigma=
Sigma, df=nu)
(with nu
being
> require(mvtnorm)
> n <- 1000
> mu <- 1:2
> Sigma <- matrix(c(4, 2, 2, 3), ncol=2)
> nu <- 3
> set.seed(271)
> X1 <- try(rmvt(n, mean=mu, sigma=Sigma, df=nu)) # error; 'mean' not allowed anymore
In mvtnorm (
To understand why it is dangerous not to throw an error in this
situation, let us look at the mean. Theoretically, the mean of
rmvt()
, a
specified argument mean
was passed to the internal call of rmvnorm()
via the ellipsis “...
”. This implies that one actually generated a
sample from
where
1.0080 | 1.0794 | 1.7757 | 8.5417 | 76.0418 | 751.0417 | |
101 | 11 | 2 | 1.1 | 1.01 | 1.001 |
It follows from (1) and (7) that previous
versions of rmvt()
with specified argument mean
actually sampled
from a random vector
The distribution of such rmvt()
, we can still mimic this
previous behavior of rmvt(n, mean=mu, sigma=Sigma, df=
nu)
by
> set.seed(271)
> ## exactly the random variates drawn by rmvt(n, mean=mu, sigma=Sigma, df=nu)
> ## in versions of mvtnorm before 0.9-9996:
> X12 <- rmvt(n, sigma=Sigma, df=nu, delta=mu, type="Kshirsagar")
> colMeans(X12) # => wrong (sample) mean
[1] 1.5380 2.7955
The result is shown on the right-hand side of
Figure 2. In contrast to many other sampling
functions in R (even including mnormt’s rmt()
),
rmvt()
does not have an argument mean
and previous versions of
rmvt()
thus generated random variates from (8) instead if
this argument was provided.
Remark 1
To properly specify the location vector, R users are often tempted to do the following:
> X2 <- mu + rmvt(n, sigma=Sigma, df=nu)
The problem with this approach is not visible for the human eye here! To make it more pronounced, let us blow up the location vector:
> set.seed(271)
> X2 <- 20*mu + rmvt(n, sigma=Sigma, df=nu)
The left-hand side of Figure 3 reveals the
problem. Indeed, we added the vector 20*mu
to the matrix returned
by rmvt()
. R solves this problem by sufficiently often
repeating the vector elements to obtain a matrix of the same size such
that addition makes sense.
> head(matrix(20*mu, nrow=n, ncol=d))
[,1] [,2]
[1,] 20 20
[2,] 40 40
[3,] 20 20
[4,] 40 40
[5,] 20 20
[6,] 40 40
As we can see (and more experienced R users know this
fact well), matrices are filled and worked on column-wise. So every
second sample has a different mean vector (alternating between (20, 20)
and (40, 40)). We thus sampled a mixture of two
In order to take care of the correct mean, there are several possibilities, some are:
> set.seed(271)
> X21 <- matrix(mu, nrow=n, ncol=d, byrow=TRUE) + rmvt(n, sigma=Sigma, df=nu)
> set.seed(271)
> X22 <- rep(mu, each=n) + rmvt(n, sigma=Sigma, df=nu)
> set.seed(271)
> X23 <- sweep(rmvt(n, sigma=Sigma, df=nu), MARGIN=2, STATS=mu, FUN="+")
> stopifnot(identical(X21, X22),
identical(X21, X23)) # check equality of the random numbers
The last approach is implemented in rmvt()
in terms of the argument
delta
(if the argument type
attains its default "shifted"
):
> set.seed(271)
> X24 <- rmvt(n, sigma=Sigma, df=nu, delta=mu)
> stopifnot(identical(X21, X24))
sigma
After having taken care of the mean vector
> set.seed(271)
> X3 <- rmvt(n, sigma=Sigma, df=nu, delta=mu)
> cov(X3)
[,1] [,2]
[1,] 9.8843 4.9204
[2,] 4.9204 7.6861
As we see, the sample covariance matrix is not close to
For X3
(displayed on the right-hand side of Figure 3)
therefore shows a sample from the correct distribution as specified by
(5); note that the same construction principle has been used
to create the left-hand side of Figure 2.
Finally, let us mention that if
> set.seed(271)
> ## sample correlation matrix of a t3 sample with scale matrix Sigma
> cor(rmvt(1e6, sigma=Sigma, df=3, delta=mu))
[,1] [,2]
[1,] 1.00000 0.57667
[2,] 0.57667 1.00000
> ## correlation matrix corresponding to Sigma
> cov2cor(Sigma)
[,1] [,2]
[1,] 1.00000 0.57735
[2,] 0.57735 1.00000
Remark 2 The user should also be aware of the fact that
rmvt(n, delta=mu, sigma=Sigma, df=0, ...)
is equivalent tormvnorm(n, mean=mu, sigma=Sigma, ...)
. This is counter-intuitive as the multivariate normal distribution arises as, not . This might be problematic for very small degrees of freedom parameters due to truncation to 0. Note that mvtnorm ( 0.9-9996) now also allows to specify df=Inf
for the (same) limiting multivariate normal distribution. The casedf=0
remains for backward compatibility.
The R package mnormt provides the function rmt()
for
sampling from the multivariate rmt(n, mean=mu, S=Sigma, df=nu)
provides the correct answer to task
(5). Note that rmt()
has an argument mean
, but actually
means the location vector (see ?rmt
). Furthermore, there is a subtle
difference between mnormt and mvtnorm. mnormt’s rmnorm()
uses a
Cholesky decomposition of rmvt()
with method="chol"
, the vectors of
random variates generated by rmt(n, mean=mu, S=Sigma, df=nu)
and those
by rmvt(n, sigma=Sigma, df=nu,
delta=mu, method="chol")
are not
identical. The reason for this is that the order in which the normal and
the rmt()
and rmvt()
.
Table 2 collects the various calls of rmvt()
and
corresponding stochastic representations of
Call | Stochastic representation of |
---|---|
X1 <- rmvt(n, mean=mu, sigma=Sigma, df=nu) ; |
gives an error |
X12 <- rmvt(n, sigma=Sigma, df=nu, delta=mu, type="Kshirsagar") ; |
|
X2 <- mu + rmvt(n, sigma=Sigma, df=nu) ; |
mixture of two |
X21 <- matrix(mu, nrow=n, ncol=d, byrow=TRUE) + rmvt(n, sigma=Sigma, df=nu) ; |
as X3 |
X22 <- rep(mu, each=n) + rmvt(n, sigma=Sigma, df=nu) ; |
as X3 |
X23 <- sweep(rmvt(n, sigma=Sigma, df=nu), MARGIN=2, STATS=mu, FUN="+") ; |
as X3 |
X24 <- rmvt(n, sigma=Sigma, df=nu, delta=mu) ; |
as X3 |
X3 <- rmvt(n, sigma=Sigma, df=nu, delta=mu) ; |
To summarize, let
The location vector
rmvt(n, sigma=Sigma, df=nu, delta=mu)
, where n
is the sample
size mu
the location vector Sigma
the scale
matrix nu
the degrees of freedom parameter
The argument sigma
of rmvt()
denotes the scale matrix cov2cor()
).
Acknowledgments
I would like to thank Martin Mächler (ETH Zürich) and Gabriel Doyon (ETH
Zürich) for valuable feedback on an early version of this work.
Furthermore, I would like to thank the Associate Editor Bettina Grün
(Johannes Kepler Universität Linz) and two anonymous reviewers for their
feedback during the review process.
mvtnorm, MASS, evir, mnormt, QRM
Distributions, Econometrics, Environmetrics, ExtremeValue, Finance, MixedModels, NumericalMathematics, Psychometrics, Robust, TeachingStatistics
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Hofert, "On Sampling from the Multivariate t Distribution", The R Journal, 2013
BibTeX citation
@article{RJ-2013-033, author = {Hofert, Marius}, title = {On Sampling from the Multivariate t Distribution}, journal = {The R Journal}, year = {2013}, note = {https://rjournal.github.io/}, volume = {5}, issue = {2}, issn = {2073-4859}, pages = {129-136} }