9.4 Multivariate Normal Distribution

In this section, we introduce how to work with multivariate normal distribution in R.

First, let’s review the definition of a multivariate normal distribution. Suppose a \(p\)-dimensional random vector \({\bf x} \sim N({\bf \mu}, {\bf \Sigma})\), where \({\bf \mu}\) is the mean vector and \({\bf \Sigma}\) is the covariance matrix. We have \[ f({\bf x}) = \frac{1}{(2\pi)^{p/2}\sqrt{|{\bf \Sigma}|}}\exp\left[-\frac{1}{2}({\bf x}-{\bf\mu})^T{\bf \Sigma}^{-1}({\bf x}-{\bf\mu})\right]. \]

9.4.1 Generate Multivate Normal Distribution

To randomly generate \({\bf x} \sim N({\bf \mu}, {\bf \Sigma})\), you can use the mvrnorm(n = 1, mu, Sigma) function in the MASS package, which is preloaded in the base R. The function mvnorm() takes three arguments.

n: targeted sample size
mu: mean vector
Sigma: covariance matrix

Let’s take a look at an example where the distribution is \(\bf\mu = (1, 2)^T\) and \[{\bf \Sigma} = \left[\begin{matrix} 1&0.5\\ 0.5&1 \end{matrix}\right]\]

library(MASS)
#> 
#> Attaching package: 'MASS'
#> The following object is masked from 'package:dplyr':
#> 
#>     select
n <- 10
mu <- c(1, 2)
Sigma <- matrix(c(1, 0.5, 0.5, 1), 2, 2)
x <- mvrnorm(n, mu, Sigma)
x
#>            [,1]     [,2]
#>  [1,] 0.5317478 1.591908
#>  [2,] 2.4575530 2.868658
#>  [3,] 0.5484207 2.079917
#>  [4,] 1.6036965 1.085302
#>  [5,] 0.7599533 2.066511
#>  [6,] 2.3854120 1.848962
#>  [7,] 1.0867794 1.785803
#>  [8,] 1.2315478 1.703268
#>  [9,] 0.7357121 1.083617
#> [10,] 0.7476221 1.690725

x contains \(n=10\) two-dimensional normal distributed variable, where each row represents one observation.

9.4.2 Sample Mean and Sample Covariance Matrix

In practical applications, we usually have a set of observations and would like to get a sample estimate of the parameters \({\bf \mu}\) and \({\bf \Sigma}\).

For the \({\bf \mu}\), you can directly use the sample mean as the estimate which is unbiased. \[\hat {\bf\mu} = [\hat\mu_j],\] where \[\hat\mu_j = \frac{1}{n}\sum_{i=1}^nx_{ij}.\]

mu_hat <- apply(x, 2, mean)
mu_hat
#> [1] 1.208844 1.780467

To estimate \({\bf\Sigma}\), we will use the sample covariance matrix estimate: \[\hat{\bf \Sigma} = [\hat\Sigma_{jk}],\] where \[\hat\Sigma_{jk} = \frac{1}{n-1}\sum_{i=1}^n(x_{ij}-\hat\mu_j)(x_{ik}-\hat\mu_k).\]

To compute this matrix, we can use the function cov() on the matrix x.

sigma_hat <- cov(x)
sigma_hat
#>           [,1]      [,2]
#> [1,] 0.5164012 0.1483479
#> [2,] 0.1483479 0.2643104
1/(n-1)*sum((x[,1]-mu_hat[1])*(x[,2]-mu_hat[2])) #compute \hat\Sigma_{12} by formula
#> [1] 0.1483479

For cov() function, in addition to computing the sample covariance matrix, you can also use it to compute the covariance among two vectors \(y\) and \(z\). \[cov(y, z) = \frac{1}{n-1}\sum_{i=1}^n(y_{i}-\bar y)(z_i-\bar z).\]

Let’s try to compute the covariance between the first and second column of \(x\).

cov(x[,1], x[,2])
#> [1] 0.1483479
cov(x)[1, 2]
#> [1] 0.1483479

You can see that this is exactly the same as \(\hat\Sigma_{12}\).

9.4.3 Sample Correlation Matrix

Another useful measure is the sample correlation matrix. It is defined as \[\widehat{corr}(x) = [\widehat{corr}(x)_{jk}],\] where \[\widehat{corr}(x)_{jk} = \frac{\hat\Sigma_{jk}} {\sqrt{\hat\Sigma_{jj}\hat\Sigma_{kk}}}. \]

You can use cor() on the matrix x to get its sample correlation matrix.

cor(x)
#>           [,1]      [,2]
#> [1,] 1.0000000 0.4015417
#> [2,] 0.4015417 1.0000000
sigma_hat[1,2] / sqrt(sigma_hat[1, 1] * sigma_hat[2, 2]) #verify the [1,2] element
#> [1] 0.4015417

Similar to the cov() function, you can also us the cor() function to compute the sample correlation between two vectors.

cor(x[,1], x[,2])
#> [1] 0.4015417
cor(x)[1, 2]
#> [1] 0.4015417

You can see that the sample correlation matrix include all the pairwise sample correlations between any of the two columns.