Statistical Simulation

Section 5.3 - Section 6.2

Bivariate/Multivariate Normal Distribution

March 30, 2026

Bivariate Normal Distribution (X1 , X2 ) is the bivariate normal distribution probability density
function:  
1 1
f (x1 , x2 ) = p exp − z(x1 , x2 ) ,
2πσ1 σ2 1 − ρ2 2(1 − ρ2 )
where
(x1 − µ1 )2 2ρ(x1 − µ1 )(x2 − µ2 ) (x2 − µ2 )2
z(x1 , x2 ) = − + .
σ12 σ1 σ2 σ22
Then, denote " # " # !
X1 µ1
∼ BN ,Σ ,
X2 µ2
where " #
σ12 ρσ1 σ2
Σ= .
ρσ1 σ2 σ22
Properties:
Cov(X1 , X2 )
1. ρ is the correlation of X1 and X2 , where ρ = .
σ1 σ2
2. The marginal distribution of X1 is N (µ1 , σ12 ), and the marginal distribution of X2 is N (µ2 , σ22 ).

3. The conditional distribution of X1 given X2 is

 
σ1 2 2
X1 |X2 = x2 ∼ N µ1 + ρ (x2 − µ2 ), (1 − ρ )σ1 .
σ2
In addition,  
σ2 2 2
X2 |X1 = x1 ∼ N µ2 + ρ (x1 − µ1 ), (1 − ρ )σ2 .
σ1
4. Recall the simple regression, E(Y |X = x) = β0 + β1 x. The estimators of β0 and β1 are:

β̂0 = ȳ − β̂1 x̄,

and
sy
β̂1 = rxy ,
sx

1
 where rxy is the sample correlation coefficient between x and y, and sx and sy are sample
standard deviations of x and y. From 3,
 
σ1 σ1 σ1
E(X1 |X2 = x2 ) = µ1 + ρ (x2 − µ2 ) = µ1 − ρ µ2 + ρ x2 ,
σ2 σ2 σ2
and V ar(X1 |X2 = x2 ) = (1 − ρ2 )σ12 .

5. If ρ = 0, then X1 and X2 are two independent random variables with X1 ∼ N (µ1 , σ12 ) and
X2 ∼ N (µ2 , σ22 ).

Figure 1: Independent bivariate normal distribution.

6. If ρ 6= 0, the figure of joint distribution of (X1 , X2 ) is as follows.

2
How to draw the bivariate normal distribution with ρ = 0?
Directly use the Box-Muller transformation and polar method.
How to draw the bivariate normal distribution with ρ 6= 0?
σ1
Use E(X1 |X2 = x2 ) = µ1 + ρ (x2 − µ2 ) and V ar(X1 |X2 = x2 ) = (1 − ρ2 )σ12 . Draw samples of
σ2
X2 from N (µ2 , σ22 ) and use the relationship between X1 and X2 .
Step 1. Draw X2∼ N (µ2 , σ22 ). 
σ1 2 2
Step 2. X1 ∼ N µ1 + ρ (X2 − µ2 ), (1 − ρ )σ1 .
σ2

Figure 2: Samples form the bivariate normal distributions.

1 The Polar Method for Generating Normal Random Variables

Y
X and Y are independent standard normal random variables. Let R2 = X 2 + Y 2 and tan Θ = X
.
Then, (R, Θ) is the polar coordinates of (X, Y ). Since X and Y are independent, the joint
density is
1 −(x2 +y2 )/2
f (x, y) = e .
2π

Let
y
d = x2 + y 2 and θ = tan−1 ,
x
and
∂d ∂d
|J| = ∂x ∂y =2
∂θ ∂θ
∂x ∂y

3
 Figure 3: Polar Coordinates.

Then,
1 1 −d/2
f (d, θ) = e ,
2 2π
1 1 −d/2
= e ,
2π 2
= fΘ (θ)fR2 (d), 0 < d < ∞, 0 < θ < 2π.
R2 and Θ are independent, where R2 is from exponential distribution with mean 2 and Θ is
from Uniform (0, 2π).
Box-Muller transformation:
Algorithm:

1. Generate random numbers U1 and U2 .

2. R2 = −2 log U1 and Θ = 2πU2 .

3. Let
p
X = R cos Θ = −2 log U1 cos(2πU2 ) (1)
p
Y = R sin Θ = −2 log U1 sin(2πU2 ) (2)

4
Polar method:
It may cost time on the computations of cos and sin. Let V1 = 2U1 − 1 and V2 = 2U2 − 1, and V1
and V2 are independent Uniform(-1,1). Note that
x V1 y V2
cos θ = = and sin θ = = ,
R r R r
where r2 = V12 + V22 , 0 < r < 1. Note that r2 is from Uniform(0, 1).

Figure 4: Polar Coordinates.

From (2), (3) and r2 is from Uniform(0, 1),

p p V1 p V1
X= −2 log U cos(2πU2 ) = −2 log U = −2 log r2 , (3)
r r
and
p p V2 p V2
Y = −2 log U sin(2πU2 ) = −2 log U = −2 log r2 . (4)
r r
Algorithm:

1. Generate a random numbers U1 and U2 .

2. R2 = −2 log U1 and Θ = 2πU2 .

3. Set V1 = 2U1 − 1, V2 = 2U2 − 1, r2 = V12 + V22 .

4. If r2 ≤ 1, set r r
−2 log r2 −2 log r2
X= V1 and Y = V2 .
r2 r2
Otherwise, return to Step 1.

5
2 Multivariate Normal Distribution

Let Z1 , . . . , Zm be independent and identically distributed normal random variables with mean 0
and variance 1. If for constants aij , i = 1, . . . , n, and j = 1, . . . , m, and mui , i = 1, . . . , n,

X1 = a11 Z1 + a12 Z2 + · · · + a1m Zm + µ1

X2 = a21 Z1 + a22 Z2 + · · · + a2m Zm + µ2
... = ...
Xi = ai1 Z1 + ai2 Z2 + · · · + aim Zm + µi
... = ...
Xn = an1 Z1 + an2 Z2 + · · · + anm Zm + µn .

Then, the vector X1 , . . . , Xn is said to have multivariate normal distribution. For all i,
E(Xi ) = µi , and

Cov(Xi , Xj ) = Cov ( m
P Pm Pm Pm
k=1 aik Zk , r=1 ajr Zr ) = k=1 r=1 Cov(aik Zk , ajr Zr )

= m
P Pm Pm
k=1 r=1 a ik a jr Cov(Z k , Zr ) = k=1 aik ajk ,

where Cov(Zk , Zr ) = 1, if k = r and Cov(Zk , Zr ) = 0, if k 6= r. It means that

 2 2 2 
a11 + a12 + · · · + a1m a11 a21 + a12 a22 + · · · + a1m a2m ··· a11 an1 + a12 an2 + · · · + a1m anm
 a21 a11 + a22 a12 + · · · + a2m a1m a221 + a222 + · · · + a22m ··· a21 an1 + a22 an2 + · · · + a2m anm 
C= .. .. .. .. .
 . . . . 
an1 a11 + an2 a12 + · · · + anm a1m an1 a21 + an2 a22 + · · · + anm a2m ··· an1 + an2 + · · · + a2nm
2 2

The matrix notation is

X0 = AZ0 + µ0 ,

where A = {aij }, X = (X1 , . . . , Xn ) is the multivariate normal vector, Z = (Z1 , . . . , Zm ) is the

row vector of independent standard normal random variables, µ = (µ1 , . . . , µn ) is the vector of
means. Let cij = Cov(Xi , Xj ), and then the covariance matrix is C = AA0 .
Generating a Multivariate Normal Random Vector

1. First, find a matrix A such that C = AA0 .

2. Generate independent standard normals (Z1 , . . . , Zm ).

3. Set X0 = AZ0 + µ0 .

6
 To find such a matrix A, we can use Choleski decomposistion. That is, there is an n × n
lower triangle matrix A such that M = AA0 and M is symmetric and positive define.
Example. Let  
9 4 2
C =  4 8 3 .
 
2 3 7
The matrix equation becomes
     
a11 0 0 a11 a12 a13 9 4 2
 a21 a22 0  ×  0 a22 a23  =  4 8 3  .
     
a31 a32 a33 0 0 a33 2 3 7

Then, 

 a211 =9 a11 = 3

a21 a11 =4 a21 = 4/3





 a31 a11 =2 a31 = 2/3
√

 a21 + a222
2
=8 a22 = 56/3 = 2.4944

a31 a21 + a32 a22 =3 a32 = 0.8463




a231 + a232 + a233

 =7 a33 = 2.4165
Questions:

1. What are the marginal distributions for X1 , X2 , and X3 ?

2. How to generate the multivariate normal distribution?

3. How to check if your multivariate draws are correct?

Exercise in class Let  

3 −2 1
C =  −2 5 3  .
 
1 3 4
Questions:

1. What are the elements in the matrix A?

2. What are the marginal distributions for X1 , X2 , and X3 ?

7
3 Gibbs Sampling

Use the full conditional distribution to draw samples from a multivariate model because

π (xj |x1 , x2 , . . . , xj−1 , xj+1 , . . . , xk ) ∝ π (x1 , x2 , . . . , xj−1 , xj , xj+1 , . . . , xk )

(t) (t) (t) (t)

At the t-step, we have (x1 , x2 , x3 , . . . , xk ).

 
(t+1) (t) (t) (t)
x1 ∼ π x1 |x2 , x3 , . . . , xk
 
(t+1) (t+1) (t) (t)
x2 ∼ π x2 |x1 , x3 , . . . , xk
 
(t+1) (t+1) (t+1) (t)
x3 ∼ π x3 |x1 , x2 , . . . , xk
..
.
 
(t+1) (t+1) (t+1) (t+1)
xk ∼ π xk |x1 , x2 , . . . , xk−1

Example: Use the Gibbs sampling for Bivariate normal distribution

Properties: The conditional distribution of X1 given X2 is

 
σ1 2 2
X1 |X2 = x2 ∼ N µ1 + ρ (x2 − µ2 ), (1 − ρ )σ1 .
σ2
In addition,  
σ2 2 2
X2 |X1 = x1 ∼ N µ2 + ρ (x1 − µ1 ), (1 − ρ )σ2 .
σ1
Algorithm:

8
Homework: Use the Gibb sampling for Multivariate normal distribution
Let (X1 , X2 , X3 )0 follow M N3 ([0, 0, 0]0 , C), where
 
3 −2 1
C =  −2 5 3  .
 
1 3 4

1. What are the marginal distributions for X1 , X2 , and X3 ?

2. What are the marginal distributions for (X1 , X2 )0 , (X1 , X3 )0 , (X2 , X3 )0 ?

3. What are the full conditional distributions:

• f (X1 |X2 , X3 );

• f (X2 |X1 , X3 );

• f (X3 |X1 , X2 ).

4. Give the algorithm by Gibbs sampling.

5. Verify your draws by Gibbs sampling are correct.