Description
IE531: Algorithms for Data Analytics
A Multivariate Gaussian RV generator has to generate i.i.d. (d × 1) vectors
x, whose (d × 1) mean is µ and whose (d × d) covariance-matrix is Σ. That is,
x ∼
1
p
(2π)
ddet(Σ)
e
− 1
2
(x−µ)
T Σ−1
(x−µ)
| {z }
=f(x)
(1)
We write the above expression in short-hand as x ∼ N(x, µ, Σ).
In this programming assignment we will write two C++ programs that can
produce a set of vectors {xi}∞
i=1, where xi ∼ N(x, µ, Σ). The approach that uses
the Cholesky decomposition of Σ (see section 1.2.1 of Probability and Statistics
Review, for example) would not work if d is very large. We will cover two (of
many) approaches to overcome this issue that uses concepts from the theory
of Markov-Chain Monte Carlo (MCMC) methods. The first approach uses the
Metropolis-Hasting Algorithm; the second approach uses Gibbs Sampling.
1 Part 1: Metropolis-Hasting Algorithm
The Proposal Distribution is the d-dimensional Multivariate Gaussian with
zero-mean (i.e. µ = 0) and Unit-Covariance (i.e. Σ = I). Assume that for right
now we have generated {xi}
j
i=1, where xi ∼ N(x, µ, Σ). We can generate RVs
xb ∼ N(x, xj , I) as xb = xj + y, where y ∼ N(y, 0, I) (which can be generated
by d-many calls to get gaussian()).
Since the Proposal Distribution is symmetric, following equation 2 of section
4 of my notes on the material in Chapter 4 of the text, we accept xb with
probability
p(xb, xj ) = min
1,
f(xb)
f(xj )
= min (
1,
e
− 1
2
(bx−µ)
T Σ−1
(bx−µ)
e
− 1
2
(x−µ)T Σ−1(x−µ)
)
. (2)
If xb is rejected, the process is repeated with no change. If xb is accepted, then
{xi}
j+1
i=1 ← {xi}
j
i=1 ∪ {xb}, and the process repeats with j ← (j + 1). That is,
we present xb as the (j + 1)-th RV xj+1 ∼ N(x, µ, Σ).
2 Part 2: Gibbs Sampling
For this part, you are going to use Gibbs Sampling to generate Multivariate
Gaussian RVs. The C++ STL has a Univariate Gaussian Generator. That is,
1
it can generate i.i.d. scalars x ∼ N(µ, σ), with mean µ and standard-deviation
σ.
The following web page1 derives the marginal- and conditional distributions
of a Multivariate Gaussian Distribution. Suppose the (d × 1) vector x is partitioned as
x =
x1
x2
where x1 is p × 1 x2 is q × 1 where p + q = d, where
x ∼ N(x, µ, Σ); µ =
µ1
µ2
; and Σ =
Σ11 Σ12
Σ21 Σ22
.
It follows that ΣT = Σ. Then, it is known that the conditional distribution of
xi given the values of xj is also Normally distributed
(xi
| xj ) ∼ N(xi
,(µi
| µj ),(Σi | Σj ))
where
(µi
| µj ) = µi + ΣijΣ
−1
jj (xj − µj )
and
(Σi | Σj ) = Σii − ΣT
jiΣ
−1
jj Σji
In your implementation of a Multivariate Gaussian RV generator that uses
Gibbs Sampling for this programming assignment, I ask that you assume p =
d − 1 and q = 1 in the above interpretation. We know that (x2 | x1) is essentially a Univariate Gaussian (with known mean and variance; that can be
computed using the above formula). By making a call to the C++ STL Univariate Gaussian RV generator, you can find a (random) value for x2. The following
YouTube video describes the MCMC-part of the algorithm. Watch it before you
start this programming assignment.
The Programming Assignment
For the first part of the programming assignment you will write a Multivariate
Gaussian Generator using the MCMC-MH algorithm. I have provided a hint
for the first part on Compass. Although the code is meant to run for highdimensions, I am verifying it for the case when d = 2. You will do the same.
The code runs on command-line, the first variable is the number of trials, the
second is the name of the CSV file that contains the experimental 2D-PDF,
the third contains the theoretical 2D-PDF. You can plot the two figures in
MATLAB, which is used in lieu of a formal verification of the correctness of the
code. Figures 1 and 2 show a sample run on my Mac.
For the second part of the programming assignment you will write a Multivariate Gaussian Generator using Gibbs Sampling. This routine takes as input
1Save a small typo.
2
Figure 1: Sample run of the MCMC-MH based Multivariate Gaussian RV Generator; d = 2 for this illustration; no of trials = 10,000,000.
Figure 2: A comparison of the experimentally observed PDF/histogram plot
(on the left) vs. the theoretical PDF (on the right) for the trial shown in figure
1 (no of trials = 10,000,000).
3
the previous sample (i.e. Previous value); and returns a d-dimensional Multivariate Gaussian RV x, where
xi =
Previous valuei
if i 6= index
∼ N(µ, σ) if i = index; and µ and σ are carefully selected.
Keep in mind, for this to work properly you have to cycle-through all d-dimensions.
That is, index ∈ {1, 2, . . . , d} in a cyclic-fashion.
In my hint.cpp I am verifying the correctness of the code by generating a 2D
Multivariate Gaussian. The code should be self-explanatory. A sample run of
the code is shown in figure 3. The experimental PDF is shown on the left of
figure 4, while the theoretical PDF is shown on the right of figure 4. Your code
should be verifiable for 2D (although it will be written for d-dimensions).
Figure 3: Sample run of the Gibbs Sampling based Multivariate Gaussian RV
Generator; d = 2 for this illustration; no of trials = 10,000,000.
4
Figure 4: A comparison of the experimentally observed PDF/histogram plot
(on the left) vs. the theoretical PDF (on the right) for the trial shown in figure
1 (no of trials = 10,000,000).
5
