Generating points from bi variate normal mixture densities

Question

I am asked to generate 200 and 1000 points from a bi-variate normal mixture densities. I am trying to understand the algorithm, not just the matlab code (I have to write it, not use an existing function). I found a code on mathworks: http://www.mathworks.com/help/matlab/ref/randn.html#bufqioz-2 on the second example, but I do not understand the math behind it. If someone can help me out with that, or explain it regardless to the code, or if you have another algorithm which can be explained, I will be grateful. Thanks a lot!

Answer 1

MATLAB very probably uses the Box-Mueller transformation which maps a pair of independent deviates uniformly distributed on (0,1) to a pair of independent unit normal deviates. The principle is this:

Consider taking uniform $x_1, x_2$ and forming $$ y_1 = \sqrt{-2 \ln x_1} cos (2\pi x_2) \\ y_1 = \sqrt{-2 \ln x_1} sin (2\pi x_2) $$

Then each of $y_1, y_2$ will be unit normal distributed, and they will be independent. The proof:

It is easy to express $x_1, x_2$ in terms of $x_1, x_2$. $$ x_1 = e^{-\frac{1}{2}(y_1^2 + y_2^2)} \\ x_2 = \tan^{-1} \frac{y_2}{y_1} $$

From that, it is easy to find the Jacobean of he transformartion (if you are unfamilear with Jacobeans, you won't be able to follow any other proof either) which is $$ -\left[ \frac{1}{\sqrt{2\pi}} e^{-y_1^2/2}\right] \left[ \frac{1}{\sqrt{2\pi}} e^{-y_1^2/2}\right] $$ which obviously factors into a term depending only on $y_1$ and a term depending only on $y_2$, so they are independent, and each has the normal distribution function as shown.

This is not the only possibility. Codes involving possible but rare rejection of the input deviates can be faster; see sections 7.3.7 and 7.3.9 of Numerical Recipes 3rd edition.

Warning For serious numerical work, the uniform random generator used by MATLAB in rand is fairly weak, and there is danger in long projects of subtle correlations causing incorrect overall results. Google the "Diehard" randoms test (created by Marsaglia) for information and more references. This is a subtle field of study.

Answer 2

that is a great answer, but i would point out that the Box-Muller transform is:

$$ R = \sqrt{-2 \ln x_1} \\ \theta =2\pi x_2 $$

and the Jacobean actually makes sense.

also, note a little confusion- of course the second line of the transformation equaion was meant to be:

$$y_2 = \sqrt{-2 \ln x_1} sin (2\pi x_2) $$

the same change in the scond part of the Jacobean.