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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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 normal distribution.

Answer 4

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).