# How to sample from a normal distribution with known mean and variance using a conventional programming language?

I've never had a course in statistics, so I hope I'm asking in the right place here.

Suppose I have only two data describing a normal distribution: the mean $\mu$ and variance $\sigma^2$. I want to use a computer to randomly sample from this distribution such that I respect these two statistics.

It's pretty obvious that I can handle the mean by simply normalizing around 0: just add $\mu$ to each sample before outputting the sample. But I don't see how to programmatically generate samples to respect $\sigma^2$.

My program will be in a conventional programming language; I don't have access to any statistical packages.

-
Does your language has random number generator? Is this generator from uniform distribution only or it can generate from normal distribution too? –  ttnphns Oct 1 '11 at 21:17
@ttnphns: Pretty much every computer language comes with a random number generator. They are overwhelmingly uniform generators on some finite domain. –  Fixee Oct 2 '11 at 0:02

If you can sample from a given distribution with mean 0 and variance 1, then you can easily sample from a scale-location transformation of that distribution, which has mean $\mu$ and variance $\sigma^2$. If $x$ is a sample from a mean 0 and variance 1 distribution then $$\sigma x + \mu$$ is a sample with mean $\mu$ and variance $\sigma^2$. So, all you have to do is to scale the variable by the standard deviation $\sigma$ (square root of the variance) before adding the mean $\mu$.

How you actually get a simulation from a normal distribution with mean 0 and variance 1 is a different story. It's fun and interesting to know how to implement such things, but whether you use a statistical package or programming language or not, I will recommend that you obtain and use a suitable function or library for the random number generation. If you want advice on what library to use you might want to add specific information on which programming language(s) you are using.

Edit: In the light of the comments, some other answers and the fact that Fixee accepted this answer, I will give some more details on how one can use transformations of uniform variables to produce normal variables.

• One method, already mentioned in a comment by VitalStatistix, is the Box-Muller method that takes two independent uniform random variables and produces two independent normal random variables. A similar method that avoids the computation of two transcendental functions sin and cos at the expense of a few more simulations was posted as an answer by francogrex.
• A completely general method is the transformation of a uniform random variable by the inverse distribution function. If $U$ is uniformly distributed on $[0,1]$ then $$\Phi^{-1}(U)$$ has a standard normal distribution. Though there is no explicit analytic formula for $\Phi^{-1}$, it can be computed by accurate numerical approximations. The current implementation in R (last I checked) uses this idea. The method is conceptually very simple, but requires an accurate implementation of $\Phi^{-1}$, which is probably not as widespread as the (other) transcendental functions log, sin and cos.
• Several answers mention the possibility of using the central limit theorem to approximate the normal distribution as an average of uniform random variables. This is not generally recommended. Arguments presented, such as matching the mean 0 and variance 1, and considerations of support of the distribution are not convincing. In Exercise 2.3 in "Introducing Monte Carlo Methods with R" by Christian P. Robert and George Casella this generator is called antiquated and the approximation is called very poor.
• There is a bewildering number of other ideas. Chapter 3 and, in particular, Section 3.4, in "The Art of Computer Programming" Vol. 2 by Donald E. Knuth is a classical reference on random number generation. Brian Ripley wrote Computer Generation of Random Variables: A Tutorial, which may be useful. The book mentioned by Robert and Casella, or perhaps Chapter 2 in their other book, "Monte Carlo statistical methods", is also recommended.

At the end of the day, a correctly implemented method is not better than the uniform pseudo random number generator used. Personally, I prefer to rely on special purpose libraries that I believe are trustworthy. I almost always rely on the methods implemented in R either directly in R or via the API in C/C++. Obviously, this is not a solution for everybody, but I am not familiar enough with other libraries to recommend alternatives.

-
Maybe you can clarify your first sentence, especially since it could give the wrong impression to a novice. :) –  cardinal Oct 1 '11 at 21:38
@cardinal, done. You were right, and I knew it. –  NRH Oct 1 '11 at 23:01
(+1) Good answer and advice for the OP. –  cardinal Oct 1 '11 at 23:02
I am not sure if I am making an unnecessary comment here, but, if you have only access to a uniform random number generator, then you can use the Box-Muller Transform to generate independent N(0,1) random numbers. In a nutshell, if U_1 and U_2 are independent draws from the Uniform(0,1) distribution then $$\sqrt{-2 \log(U_1) } \cos(2\pi U_2)$$ and $$\sqrt{-2 \log(U_1) }\sin(2\pi U_2)$$ are distributed as independent N(0,1) random variables. The basic idea –  VitalStatistix Oct 1 '11 at 23:12
@Vital: Not an unnecessary comment; a good one. The Box-Muller transform is probably the very easiest to program with minimal chance of inadvertently doing something bad. It's not the fastest, but it's competitive enough. That said, using an established code library is probably safer still, especially since the place where one is most likely to make a misstep is in how the uniform random variate inputs are generated! –  cardinal Oct 2 '11 at 0:07

This is really a comment on Michael Lew's answer and Fixee's comment, but is posted as an answer because I don't have the reputation on this site to comment.

The sum of twelve independent random variables uniformly distributed on $[0, 1]$ has mean $6$ and variance $1$. In other words, $$E\left [\sum_{i=1}^{12} X_i\right ] = \sum_{i=1}^{12} E[X_i] = 12\times \frac{1}{2} = 6$$ and $$\text{var} \left [\sum_{i=1}^{12} X_i\right ] = \sum_{i=1}^{12} \text{var}[X_i] = 12\times \frac{1}{12} = 1.$$ The CLT can then be used to assert that the distribution of $\sum_{i=1}^{12} X_i - 6$ is approximately a standard normal distribution. Compared to the ten variables considered by Michael Lew and Fixee, two additional calls to the random number generator are required, but we avoid division by $\sqrt{10/12}$ to get the desired unit variance. It is also worth remembering that $\sum_{i=1}^{12} X_i - 6$ can take on values only in the range $[-6, 6]$ and thus extreme (very low-probability) values differing from the mean by more than $6$ standard deviations will never occur. This is often a problem in simulations of computer and communication systems where such very low probability events are of much interest.

-

In addition to the answer by NRH, if you still have no means to generate random samples from a "standard normal distribution" N(0,1), below is a good and simple way (since you mention you don't have a statistical package, the functions below should be available in most standard programming languages).

1. Generate u and v as two uniformly distributed random numbers in the range from -1 to 1 by
u = 2 r1 - 1 and v = 2 r2 - 1

2.calculate w = u^2 + v^2 if w > 1 the go back to 1

3.return u*z and y= v*z with z= sqrt(-2ln(w)/w) A sample code would look like this:

u = 2 * random() - 1;
v = 2 * random() - 1;
w = pow(u, 2) + pow(v, 2);
if (w < 1) {
z = sqrt((-2 * log(w)) / w);
x = u * z;
y = v * z;
}


then use what MHR has suggested above to obtain the random deviates from N(mu, sigma^2).

-
When I posted my answer above I didn't notice that @vitalStatistix gave you the Box-Muller Transform algorithm. The one I give above is also as good I suppose. –  francogrex Oct 2 '11 at 9:29
Could you please explain the reason for generating normal variates from uniform distribution (other than from an algorithmic perspective) and not just using the pdf of a Gaussian/Normal distribution directly? Or is it totally wrong? –  Arun Oct 2 '11 at 10:14
@Arun One reason: The Marsaglia's polar method is useful when you only have a RNG that generates uniform deviates. –  chl Oct 2 '11 at 10:42
@Arun it is the easiest way. You can also generate from the pdf directly using for example the "acceptance rejection" method. I posted for you a simple example on my site (because not enough space in the comment box here). –  francogrex Oct 2 '11 at 13:10
francogrex, chl, thank you! –  Arun Oct 2 '11 at 15:25

The normal distribution emerges when one adds together a lot of random values of similar distribution (similar to each other, I mean). If you add together ten or more uniformly distributed random values then the sum is very nearly normally distributed. (Add more than ten if you want it to be even more normal, but ten is enough for almost all purposes.)

Say that your uniform random values are uniformly distributed between 0 and 1. The sum will then be between 0 and 10. Subtract 5 from the sum and the mean of the resulting distribution will be 0. Now you divide the result by the standard deviation of the (near) normal distribution and multiply the result by the desired standard deviation. Unfortunately I'm not sure what the standard deviation of the sum of ten uniform random deviates is, but if we are lucky someone will tell us in a comment!

I prefer to talk to students about the normal distribution in these terms because the utility of the assumption of a normal distribution in many systems stems entirely from the property that the sums of many random influences leads to a norml distribution.

-
You are using the Central Limit Thm here (that a bunch of iid random variables sum to a normal random variable). I didn't consider this because I thought it would be too slow, but you say 10 is sufficient?! This is better than computing a log and a sin/cos and a sqrt! –  Fixee Oct 3 '11 at 0:01
Also, the mean of the uniform r.v. on [0,1] is 0.5 with variance 1/12. If you sum 10 of these you get a mean of 5 and variance 10/12 = 5/6. –  Fixee Oct 3 '11 at 0:03
From a pedagogical standpoint this method provides for a nice, useful discussion and demonstration. However, I would strongly discourage anyone from using this approach in practice. –  cardinal Oct 3 '11 at 18:41
@Fixee: You need to be sure and balance the computation of $\log$, $\sin$, $\cos$ and the square-root against the generation of additional uniform random variates. For example, Intel CPUs have all four of these functions as built-in operations performed in hardware. The square-root is a fundamental "arithmetic" operation according to the IEEE 754 standards. –  cardinal Oct 3 '11 at 19:01
@cardinal: Good point! But in some processors you don't have these instructions (and it's a lot of work to do the Taylor Series to compute them) whereas a Linear Congruential Generator is very cheap/simple to implement. So I guess it depends on your setting, as you say. :) –  Fixee Oct 3 '11 at 21:22

Using random number generator of your programming language, try this -

1. Obtain the pdf of normal dist. with that mean and variance. Write it as a function of p(x). Example -

x = f(p(x))

1. Choose a random number between 0 and 1. Say this is p. What value of x will p(x)=p.(You can feed the equation). Now you have x. You know that x has p% chances of getting selected.

2. Choose another random number r between 0 and 1. If r < p, then choose x or else ignore it.

3. Keep iterating till you get desired number of numbers.

-
Welcome to the site, @AdityaKumarSingh. I don't quite follow your answer. It seems like you are talking about both using the inverse distribution function to output a variate given a probability, & talking about an accept / reject algorithm, as though they are part of the same process. I wonder if you can clarify your answer. Also, since you are new to the site, be sure to read our FAQ. –  gung Oct 4 '12 at 2:50