Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Bear with me as I am not working on all cylinders today.

I am attempting to write a Gaussian random vector generator using the probability distribution function:

$N(x|\mu,\sum) = {1 \over (2\pi)^{D/2}} {1 \over |{\sum}|^{1/2}}exp\{-{1 \over 2}(x - \mu)^{T}\sum^{-1}(x - \mu)\}$

(I should note that I am writing this in python.)

What I do is I create a random vector with $n$ dimensions, then send it into the equation and get a density matrix of $n$ by $n$ (I believe this is what I should get).

This is the point where I am not sure what to do. For the "1 dimensional" it seems simple enough, I get a 1 by 1 matrix that contains the probability that the random vector (value) I created exists in the distribution, and I compare that to another randomly generated probability and if it exceeds that value, then I keep it.

I am not sure however, what to do with the cases of higher dimensions (i.e. when I get a density matrix vs. a single value). Is it just a matter of taking the trace of the matrix? Or is there something else that I am missing?

share|cite|improve this question
up vote 1 down vote accepted

In "I create a random vector with $n$ dimensions", you don't specify what distribution you're sampling the random vector from. It's also not clear to me what you mean when you say that you get an $n\times n$ density matrix. Substituting a vector $x$ into the probability distribution function yields a density, which is a real number; there's no difference in that respect between the one-dimensional and the multi-dimensional cases. If you really wanted to do rejection sampling (you shouldn't; see below), you'd have to first divide that number by the density of the distribution you're sampling from, then divide by the maximal value of that ratio over all vectors, and then compare the result to a random number drawn from the uniform distribution over $[0,1]$.

But it's inefficient to sample from a Gaussian distribution using rejection sampling, anyway, since it can be done directly in closed form; see Wikipedia.

share|cite|improve this answer
This sort of answers the question. My issue was that my math was incorrect, which is basically what you said. Especially the portion in the exp{}. (The end result is a single value as you and I were expecting.) The formula is correct, I was simply using incorrect math library functions that weren't resulting in the correct calculations, so the values were not what I expected. Thanks. – NominSim Oct 27 '12 at 0:32

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.