Chapter 6 Multivariate Normal Distributions
To date, we have studied a number of standard probability distributions: uniform, Bernoulli, binomial, geometric, negative binomial, Poisson, exponential, gamma, normal. These are all random variables that govern exactly one quantity. However throughout this course we have studied joint probability distribution functions: that is the study of multiple quantities simultaneously. This begs the question, are there any standard probability distributions in higher dimensions?
6.1 Matrices: symmetry, eigenvalues and postive definiteness
Throughout this chapter, we will use the language of matrices freely. As a point of notation, we will represent all matrices with variable names in a bold font, for example \(\mathbf{A}\). A vector will be thought of as a matrix in this sense: it is a matrix with either \(1\) row or \(1\) column.
Further when we use the notation \(\mathbf{A} = (a_{ij})\), this labels the entry of \(\mathbf{A}\) in the \(i^{th}\) row and \(j^{th}\) column by \(a_{i,j}\).
We consolidate our knowledge on the definition of symmetric matrices, the determinant of an \(n \times n\) matrix, the definitions of eigenvalues and eigenvectors, and what it means for matrices to be positive definite.
Let \(\mathbf{A}\) be a \(n \times n\) matrix.Then \(\mathbf{A}\) is said to be symmetric if \(\mathbf{A} = \mathbf{A}^{T}\).
Definition 6.1.1 is equivalent to \(\mathbf{A}=(a_{ij})\) satisfying \(a_{ij}=a_{ji}\) for all \(i,j\) with \(1\leq i,j \leq n\).
Are the following matrices symmetric? \[1. \quad \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}; \hspace{40pt} 2. \quad \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix};\] \[3. \quad \begin{pmatrix} 0 & 3 & 4 \\ -3 & 0 & 7 \\ -4 & -7 & 0 \end{pmatrix}; \hspace{40pt} 4. \quad \begin{pmatrix} 9 & 2 & 3 \\ 2 & -1 & -8 \\ 3 & -8 & 0 \end{pmatrix}.\]
The answers are
Yes;
Yes;
No;
Yes.
Let \(\mathbf{A}\) be a \(n \times n\) matrix. A scalar \(\lambda\) is called an eigenvalue of \(\mathbf{A}\) if there is a non-zero vector \(\mathbf{v}\) such that \[\mathbf{Av} = \lambda \mathbf{v}.\] The vector \(\mathbf{v}\) corresponding to \(\lambda\) is called an eigenvector of \(\mathbf{A}\).
It was shown that an eigenvalue \(\lambda\) of an \(n \times n\) matrix \(\mathbf{A}\) satisfies the characteristic polynomial: \[\det(\mathbf{A} - \lambda \mathbf{I_n})=0.\] Recall that the determinant of any \(n \times n\) matrix \(\mathbf{M}\) is calculated by fixing some value \(j\) with \(1\leq j \leq n\), and finding \[\det (\mathbf{M}) = \sum\limits_{i=1}^{n}(-1)^{i+j}a_{ij}M_{ij},\] where \(M_{ij}\) is the minor of \(\mathbf{M}\), that is, the determinant of the \((n-1) \times (n-1)\) matrix obtained by deleting the \(i^{th}\) row and \(j^{th}\) column of \(\mathbf{M}\).
What are the eigenvalues and eigenvectors of \(\mathbf{A} = \begin{pmatrix} 5 & 2 & 0 \\ 2 & 5 & 0 \\ -3 & 4 & 6 \end{pmatrix}\)?
Note that \(\mathbf{A}\) is a \(3\times 3\) matrix, so in the notation of Definition 6.1.3, \(n=3\). Calculate
\[\begin{align*}
A - \lambda I_3 &= \begin{pmatrix} 5 & 2 & 0 \\ 2 & 5 & 0 \\ -3 & 4 & 6 \end{pmatrix} - \lambda \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \\[9pt]
&= \begin{pmatrix} 5-\lambda & 2 & 0 \\ 2 & 5-\lambda & 0 \\ -3 & 4 & 6-\lambda \end{pmatrix}.
\end{align*}\]
So \[\begin{align*} \det (A - \lambda I_3) &= 0 \\[9pt] \det \begin{pmatrix} 5-\lambda & 2 & 0 \\ 2 & 5-\lambda & 0 \\ -3 & 4 & 6-\lambda \end{pmatrix} &=0 \\[9pt] (5-\lambda) \big[ (5-\lambda)(6-\lambda) - 0 \big] - 2 \big[ 2(6-\lambda) -0 \big] &=0 \\[9pt] (5-\lambda) (5-\lambda)(6-\lambda) - 4(6-\lambda) &= 0 \\[9pt] (6-\lambda) \big[ (5-\lambda)(5-\lambda) - 4\big] &=0 \\[9pt] (6-\lambda) (\lambda^2 -10 \lambda +21) &=0 \\[9pt] (6-\lambda) (\lambda-3)(\lambda-7) &=0 \end{align*}\]
So \(\lambda_1=3, \lambda_2 =6, \lambda_3 =7\) are the eigenvalues.
For \(\lambda_1=3\), we are look for a eigenvector \(v_1 = \begin{pmatrix} x \\ y \\ z \end{pmatrix}\) satisfying \[\begin{align*} &(A - 3I_3) v_1 = \mathbf{0} \\[9pt] \iff \hspace{10pt} &\begin{pmatrix} 2 & 2 & 0 \\ 2 & 2 & 0 \\ -3 & 4 & 3 \end{pmatrix}\begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \\[9pt] \iff \hspace{10pt} &\begin{cases} 2x+2y=0 \hspace{100pt} (1)\\ 2x+2y=0 \\ -3x + 4y +3z=0 \hspace{70pt} (2) \end{cases} \end{align*}\]
Equation (1) rearranges to \(y=-x\). Substituting this into equation (2), we obtain \[\begin{align*} -3x-4x+3z &= 0 \\[9pt] 3z &= 7x \end{align*}\]
So clearly \(x=3, y=-3\) and \(z=7\) is a solution. So we have an eigenvector \[v_1 = \begin{pmatrix} 3 \\ -3 \\ 7 \end{pmatrix}.\]
For \(\lambda_2 = 6\), we are looking for eigenvector \(v_2 = \begin{pmatrix} x \\ y \\ z \end{pmatrix}\) satisfying \[\begin{align*} &(A - 6I_3) v_1 = \mathbf{0} \\[9pt] \iff \hspace{10pt} &\begin{pmatrix} -1 & 2 & 0 \\ 2 & -1 & 0 \\ -3 & 4 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \\[9pt] \iff \hspace{10pt} &\begin{cases} -x+2y=0 \\ 2x-y=0 \\ -3x + 4y=0 \end{cases} \end{align*}\]
Similarly to the previous eigenvalue, we solve these three equations simultaneously giving \[v_2 = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}.\]
For \(\lambda_3 = 7\), we are looking for eigenvector \(v_3 = \begin{pmatrix} x \\ y \\ z \end{pmatrix}\) satisfying \[\begin{align*} &(A - 7I_3) v_3 = \mathbf{0} \\[9pt] \iff \hspace{10pt} &\begin{pmatrix} -2 & 2 & 0 \\ 2 & -2 & 0 \\ -3 & 4 & -1 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \\[9pt] \iff \hspace{10pt} &\begin{cases} -2x+2y=0 \\ 2x-2y=0 \\ -3x + 4y-z=0 \end{cases} \end{align*}\]
Similarly to the first eigenvalue, we solve these three equations simultaneously giving \[v_3 = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}.\]
Therefore \(\mathbf{A}\) has eigenvalues \(3,6,7\) with corresponding eigenvectors given respectively by \(\begin{pmatrix} 3 \\ -3 \\ 7 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}, \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}\).A symmetric \(n \times n\) matrix \(\mathbf{A}\) is positive definite if the quantity \(\mathbf{z}^{T} \mathbf{A z}\) is positive for all column vectors \(\mathbf{z}\) of length \(n\).
Note that \(\mathbf{A z}\) is the transformation of \(z\) under \(A\). It follows that \(\mathbf{z}^{T} \mathbf{A z}\) is the dot product of \(\mathbf{z}\) and \(\mathbf{A z}\). This being positive is equivalent to the angle between \(\mathbf{z}\) and \(\mathbf{A z}\) being less than \(\frac{\pi}{2}\). Therefore a positive definite matrix \(\mathbf{A}\) is informally a matrix that transforms all vectors \(\mathbf{z}\) in such a way that they still point in the same general direction.
Determine whether the matrix \(\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\) is positive definite.
For any vector \(\mathbf{z} = \begin{pmatrix} x \\ y \end{pmatrix}\), we have \[\mathbf{z}^T \mathbf{Az} = \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = x^2 + y^2.\] Note for all real numbers \(x,y\), we have \(x^2 \geq 0\) and \(y^2 \geq 0\). It follows that \(\mathbf{z}^T \mathbf{Az} \geq 0\) for all \(\mathbf{z}\).
Determine whether the matrix \(\begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix}\) is positive definite.
Set \(\mathbf{z} = \begin{pmatrix} -1 \\ 1 \end{pmatrix}\). We have \[\mathbf{z}^T \mathbf{Az} = \begin{pmatrix} -1 & 1 \end{pmatrix} \begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix} \begin{pmatrix} -1 \\ 1 \end{pmatrix} = \begin{pmatrix} -1 & 1 \end{pmatrix} \begin{pmatrix} 1 \\ -1 \end{pmatrix} = -2 <0.\] Therefore \(\begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix}\) is not positive definite.
Matrix calculations can be performed in R. To create a matrix variable you can use the matrix function. As input the function takes a vector, below this is `c(2, 4, 3, 1, 5, 7)’, representing the entries, and a parameter nrows representing the number of rows the entries should be sorted into.
R also allows for algebraic manipulations of matrices such as matrix addition, matrix multiplication and scalar multiplication.
Finally R allows for easy calculation of eigenvalues and eigenvectors of a square matrix.
D = matrix( c(13, -4, 2, -4, 11, -2, 2, -2, 8), nrow=3)
ev <- eigen(D)
values <- ev$values
vectors <- ev$vectors
values
vectorsR also allows for algebraic manipulations of matrices
6.2 Definition of Multivariate Normal Distribution
The normal distribution that we saw in Section 1.3 generalises to higher dimensions. Let \(n \in \mathbb{Z}_{\geq 1}\) represent the dimension we are interested in.
Informally a multivariate normal distribution can be defined as follows: Consider \(n\) standard normal distributions \(Z_1, Z_2, \ldots, Z_n\), that is, \(Z_i \sim N(0,1)\) for all \(i\) with \(1 \leq i \leq n\). From this, Form a vector of random variables \(\mathbf{Z} = \begin{pmatrix} Z_1 & Z_2 & \cdots & Z_n \end{pmatrix}^{T}\). Then for any \(n \times n\) non-singular matrix \(\mathbf{A}\) and any vector \(\mathbf{\mu} = \begin{pmatrix} \mu_1 & \mu_2 & \cdots & \mu_n \end{pmatrix}^T\), define \(\mathbf{X} = \mathbf{AZ}+\mathbf{\mu}\). Then the vector of random variables \(\mathbf{X} = \begin{pmatrix} X_1 & X_2 & \cdots & X_n \end{pmatrix}^T\) is distributed by a multivariate normal distribution.
This is saying that every component \(X_i\) is a linear combination of some fixed group of standard normal distributions, and a translation.
Formally we define the multivariate normal distribution by the joint probability density function of \(X_1, X_2, \ldots, X_n\).
Let \(\mathbf{\mu} = \begin{pmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_n \end{pmatrix}\) and \(\mathbf{\Sigma} = (\sigma_{ij})\) be an \(n \times n\) real, symmetric, positive definite matrix all of whose eigenvalues are positive.
A vector of random variables \(\mathbf{X} = \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_n \end{pmatrix}\) is said to have an \(\mathit{n}\)-dimensional normal distribution with parameters \(\mathbf{\mu}\) and \(\mathbf{\Sigma}\), denoted \(N_n(\mathbf{\mu},\mathbf{\Sigma})\), if the joint PDF of \(\mathbf{X}\) is given by \[ f_{\mathbf{X}}(\mathbf{x}) = \frac{1}{\sqrt{(2\pi)^n\cdot \det(\mathbf{\Sigma})}} \exp \left\{ -\frac{1}{2} (\mathbf{x}-\mathbf{\mu})^T \mathbf{\Sigma}^{-1} (\mathbf{x}-\mathbf{\mu}) \right\}.\]
The PDF outlined in the formal Definition 6.2.1 can be derived from the informal construction of \(\mathbf{X}\) at the beginning of the section by taking \(\mathbf{\Sigma} = \mathbf{AA}^T\).
Note for \(n=1\), Definition 6.2.1 simplifies to give the PDF of the one-dimensional normal distribution.
The multivariate normal distribution has the following important properties:
The expectation of \(X_i\) is given by the \(i^{th}\) entry of \(\mathbf{\mu}\), that is, \(E[X_i]=\mu_i\) for all \(i\). Succinctly, the vector \(\mathbf{\mu}\) is the expectation vector of \(X\);
The covariance \(\text{Cov}(X_i,X_j)\), reducing to \(\text{Var}(X_i)\) when \(i=j\), is given by the element in the \(i^{th}\) row and \(j^{th}\) column of \(\mathbf{\Sigma}\), that is, \(\sigma_{ij} = \text{Cov}(X_i,X_j)\) for all \(i\) and \(j\). Succinctly, the matrix \(\Sigma\) is the variance-covariance matrix of \(\mathbf{X}\). It follows that if \(X_1,X_2,\dots,X_n\) are independent then \(\sigma_{ij}=0\) for all \(i \neq j\);
These two points indicate why \(\mathbf{\mu}\) and \(\mathbf{\Sigma}\) are the two input parameters for the multivariate normal distribution.
6.3 Two-Dimensional Normal Distribution
In this section we study the \(2\)-dimensional normal distribution, often referred to as the bivariate normal distribution.
When \(n=2\), or another suitably small value, we can explicitly label the elements of \(\mathbf{\mu}\) and \(\mathbf{\Sigma}\) as \(\begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix}\) and \(\begin{pmatrix} \sigma_1^2 & \sigma_1 \sigma_2 \rho \\ \sigma_1 \sigma_2 \rho & \sigma_2^2 \end{pmatrix}\) respectively. Here \(\sigma_1\) is the standard deviation of \(X_1\), \(\sigma_2\) is the standard deviation of \(X_2\) and \(\rho\) is the correlation coefficient \(\rho(X_1,X_2)\).
Why are the upper right and lower left entries of \(\mathbf{\Sigma}\) equal to \(\sigma_1 \sigma_2 \rho\)?
Show that \(\Sigma\) is non-singular if and only if \(\lvert \rho \rvert <1\) and \(\sigma_1\sigma_2 \neq 0\).
Substitution of these values into Definition 6.2.1 gives \[\begin{align*} f_{X_1,X_2}(x_1,x_2) = \frac{1}{2 \pi \sigma_1 \sigma_2 \sqrt{1-\rho^2}} exp \left\{ - \frac{1}{2(1-\rho^2)} \left[ \left( \frac{x_1-\mu_1}{\sigma_1} \right)^2 \\[5pt] \hspace{150pt}-2\rho \left( \frac{x_1-\mu_1}{\sigma_1} \right)\left( \frac{x_2-\mu_2}{\sigma_2} \right) + \left( \frac{x_2-\mu_2}{\sigma_2} \right)^2 \right] \right\} \end{align*}\]
The following code creates a contour plot of the PDF above with \(\mathbf{\mu} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\) and \(\mathbf{\Sigma} = \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}\).
library(mnormt)
#create bivariate normal distribution
x <- seq(-3, 3, 0.1)
y <- seq(-3, 3, 0.1)
mu <- c(0, 0)
sigma <- matrix(c(2, -1, -1, 2), nrow=2)
f <- function(x, y) dmnorm(cbind(x, y), mu, sigma)
z <- outer(x, y, f)
#create contour plot
contour(x, y, z)A sample of size \(200\) can be observed from a bivariate normal distribution using the following R code.
library(plotly)
library(mixtools)
# Set up means, standard deviations, correlations
mu1 <- 1; mu2 <- 2
sig1 <- 1; sig2 <- 1; rho <- 0;
# Construct mean vector and (co)variance matrix
mu <- c(mu1,mu2)
Sigma <- matrix(c(sig1^2,rho*sig1*sig2,rho*sig1*sig2,sig2^2), ncol=2)
#plot all realisations
X=rmvnorm(n=200, mu, Sigma)
plot(X)In the above code, by setting mu as a vector of length \(n\), and \(Sigma\) as a \(n \times n\) matrix before inputting them into rmvnorm allows one to simulate an \(n\)-dimensional normal distribution.
6.4 Properties of Multivariate Normal Distribution
The following theorem states that the transformation of a multivariate normal distribution is itself a multivariate normal distribution.
Let \(\mathbf{X} \sim N_n(\mathbf{\mu},\mathbf{\Sigma})\) and consider some \(p \times n\) matrix \(\mathbf{D}\). Define a new random vector \(\mathbf{Z} = \mathbf{D}\mathbf{X}\). Then \(\mathbf{Z} \sim N_p (\mathbf{D} \mathbf{\mu}, \mathbf{D} \mathbf{\Sigma} \mathbf{D}^T)\);
As a corollary of this theorem, we can easily calculate the marginal distributions of a multivariate normal distribution.
The marginal distribution of each component \(X_i\) is a one-dimensional normal distribution.
This is a direct consequence of Theorem 6.4.1 by setting \(\mathbf{D}=(0,\dots,0,1,0,\dots,0)\) where the \(1\) is in the \(i^{th}\) position.
Suppose \(\mathbf{X}=(X_1,X_2,X_3)^T \sim N_3(\mathbf{0},\mathbf{\Sigma})\), where \[\mathbf{\Sigma} = \begin{pmatrix} 2 & 1 & 0 \\ 1 & 4 & 0 \\ 0 & 0 & 5 \end{pmatrix}.\] Find the distribution of \(Z=X_1+X_2\).
Suppose \(\mathbf{X}=(X_1,X_2,X_3)^T \sim N_3(\mathbf{0},\mathbf{\Sigma})\), where \[\mathbf{\Sigma} = \begin{pmatrix} 2 & 1 & 0 \\ 1 & 4 & 0 \\ 0 & 0 & 5 \end{pmatrix}.\] Determine the constant \(c\) such that \(Z_1 = 2X_1 + cX_2\) and \(Z_2 = 2X_1 + cX_3\) are independent.
Writing \(\mathbf{Z} = \begin{pmatrix} Z_1 \\ Z_2 \end{pmatrix} = \begin{pmatrix} 2X_1 + cX_2 \\ 2X_1 + cX_3 \end{pmatrix}\) in the form \(\mathbf{DX}\) requires
\[\mathbf{D} = \begin{pmatrix} 2 & c & 0 \\ 2 & 0 & c \end{pmatrix}.\]
By the properties of a multivariate normal distribution,
\[\mathbf{Z} \sim N_2(\mathbf{D0}, \mathbf{D \Sigma D}^T).\]
Calculate \(\mathbf{D0} = \mathbf{0}\) and
\[\begin{align*} \mathbf{D \Sigma D}^T &= \begin{pmatrix} 2 & c & 0 \\ 2 & 0 & c \end{pmatrix} \begin{pmatrix} 2 & 1 & 0 \\ 1 & 4 & 0 \\ 0 & 0 & 5 \end{pmatrix} \begin{pmatrix} 2 & 2 \\ c & 0 \\ 0 & c \end{pmatrix} \\[5pt] &= \begin{pmatrix} 4+c & 2+4c & 0 \\ 4 & 2 & 5c \end{pmatrix} \begin{pmatrix} 2 & 2 \\ c & 0 \\ 0 & c \end{pmatrix} \\[5pt] &= \begin{pmatrix} 8+4c+4c^2 & 8+2c \\ 8+2c & 8+5c^2 \end{pmatrix}. \end{align*}\]
For \(Z_1\) and \(Z_2\) to be independent necessarily \(\text{Cov}(Z_1,Z_2) = 8+2c = 0\). Therefore \(c=-4\).