Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit

STAT 3690 Lecture 06

2022

Multivariate normal (MVN) distribution (J&W Sec 4.2)

●  Standard normal random vector

φZ (之) = (2π)-p/2 exp(-之T 之/2),    之 = [z1 , . . . , zp ]T  e 皿p

●  (General) normal random vector

-  Def: The distribution of X is MVN iff there exists q e 勿+ , u e 皿q , A e 皿qxp  and Z ~ MVNp (0I) such that X = AZ + u

*  Limit the discussion to non-degenerate cases, i.e., rk(A) = q

*  X ~ MVNq (u, Σ), i.e.,

fx () =  exp(-( - u)T Σ- 1 ( - u)/2},     e q

'  Σ = var(X) = AAT  > 0 

●  Exercise:

1.  Σ = AAT  > 0 兮 rk(A) = q (Hint: SVD of A);

2.  Σ > 0 三 there exists a p x p positive definite matrix, say Σ 1/2, such that Σ = Σ 1/2Σ 1/2  and Σ- 1  = Σ- 1/2Σ- 1/2  (Hint: spectral decomposition of Σ). 

options(digits  =  4)

(Sigma  =  matrix(c ( 10 ,  2  ,2 ,  5),  ncol  =  2 ,  nrow  =  2))

(spectral  =  eigen(Sigma))

(SigmaRoot  =  spectral$vectors  %*%  diag(spectral$values ˆ .5)  %*%  t(spectral$vectors))        (SigmaRootInv  =  spectral$vectors  %*%  diag(spectral$values ˆ- .5)  %*%  t(spectral$vectors))

#  Check  properties  of  root  of  Sigma

(SigmaRoot  %*%  SigmaRoot  -  Sigma)

(solve(SigmaRoot)  -  SigmaRootInv)

(SigmaRootInv  %*%  SigmaRootInv  -  solve(Sigma))

#  SVD  <=>  spectral  decomposition  if  Sigma  is  positive   (semi-)definite

svd(Sigma)

eigen(Sigma)

●  Useful properties of MVN

-  X ~ MVNp (u, Σ) 兮 Z = Σ- 1/2(X - u) ~ MVNp (0I). So, we have a stochastic representation of arbitrary X ~ MVNp (u, Σ): X = Σ 1/2Z + u, where Z ~ MVNp (0I).

-  X ~ MVN iff, for all a e 皿p , aT X has a (univariate) normal distribution.

-  If X ~ MVNp (u, Σ), then AX + b ~ MVNq (Au + b, AΣAT ) for A e 皿qxp  and rk(A) = q .

●  Exercise: Generate six iid samples following bivariate normal MVN2 (u, Σ) with u = [3, 6]T ,    Σ = –  2(10)   5(2) ! .

options(digits  =  4)

set.seed(1)

Mu  =  matrix(c (3 ,  6),  ncol  =  1 ,  nrow  =  2)

Sigma  =  matrix(c ( 10 ,  2  ,2 ,  5),  ncol  =  2 ,  nrow  =  2)

n  =  1000

# Method  1:  via  rnorm()

spectral  =  eigen(Sigma)

SigmaRoot  =  spectral$vectors  %*%  diag(spectral$values ˆ .5)  %*%  t(spectral$vectors) A1  =  matrix(0 ,  nrow  =  n,  ncol  =  length (Mu))

for  (i  in  1:n)  {

A1[i,  ]  =  t(SigmaRoot  %*%  matrix(rnorm (2),  nrow  =  2 ,  ncol  =  1)  +  Mu)

}

#  Method  2:  via  MASS::mvrnorm()

A2  =  MASS::mvrnorm(n,  Mu,  Sigma)