Generating a random gaussian number

[tony873004]

I want write a short block of computer code that allows me to generate a random gassuian number. For example, if the inputs are mean = 100 and sd = 10, if I generage 1000 random numbers, I'd like 67% of them to be between 90 and 110, 98% of them to be be between 80 and 120. But I only know how to generate random numbers between 2 different values. For example I can generate random numbers between 90 and 110, but 90 would have just as as much probability of being generated as 100. How can I generate numbers that would make a nice bell curve if I generated enough of them?

Thanks.

 

[CompuChip]

Something like this? I don't know what language you are looking for, but Googling "generate gaussian random number" gives me tons of results (about 1.6 million, actually)

 

[D H]

The standard approach to generating normal deviates is the Box-Müller algorithm. The polar form of the algorithm is very robust and has withstood many tests of normality.

Here is a non-reentrant version of the algorithm. It assumes a function random_uniform that returns a random deviate drawn from U(0,1).

extern double random_uniform();
double random_gaussian() {
   static int have_deviate = 0;
   static double u1, u2;
   double  x1, x2, w;

   if (have_deviate) {
      return u2;
      have_deviate = 0;
   } else {
      do {
         x1 = 2.0 * random_uniform() - 1.0;
         x2 = 2.0 * random_uniform() - 1.0;
         w  = x1 * x1 + x2 * x2;
      } while ( w >= 1.0 );
      w  = sqrt ((-2.0 * ln (w))/w);
      u1 = x1 * w;
      u2 = x2 * w;
      return u1;
      have_deviate = 1;
   }
}

 

[tony873004]

The standard approach to generating normal deviates is the Box-Müller algorithm. The polar form of the algorithm is very robust and has withstood many tests of normality.

Here is a non-reentrant version of the algorithm. It assumes a function random_uniform that returns a random deviate drawn from U(0,1).

extern double random_uniform();
double random_gaussian() {
   static int have_deviate = 0;
   static double u1, u2;
   double  x1, x2, w;

   if (have_deviate) {
      return u2;
      have_deviate = 0;
   } else {
      do {
         x1 = 2.0 * random_uniform() - 1.0;
         x2 = 2.0 * random_uniform() - 1.0;
         w  = x1 * x1 + x2 * x2;
      } while ( w >= 1.0 );
      w  = sqrt ((-2.0 * ln (w))/w);
      u1 = x1 * w;
      u2 = x2 * w;
      return u1;
      have_deviate = 1;
   }
}

Thanks, CompuChip and DH.

In the above code, what is the purpose of computing u2 if only u1 is returned and u2 isn't used again? I took CompuChip's advice and Googled this method to read a little more. It seems it assumes a mean of 0 and a standard deviation of 1. How can I modify this so I can specify the mean and standard deviation? And I'm only expecting one output. Is that in u1, u2, or do I have to perform an additional step once this function has been executed?

Also, if I'm getting data from a table and it reads mean=100, uncertainty (1sd)= 10, does 10 represent the distance from the mean to 1sd, or the distance from -sd to sd. It's been years since I've taken a stats class. Thanks!

 

[tony873004]

I figured it out. Multiply u1 by standard deviation and adding mean to it gives me my random output. I'm still not sure what good it did to compute u2, but if I ignore it everything works fine. Now I can plot beautiful gaussian curves. Thanks for your help!

 

[D H]

The routine I supplied computes two guassian deviates, u1 and u2 and returns one of them (u1) on the first time it is called. On the second call, it simply returns u2 without computing anything. On the third call it computes two new gaussian deviates and returns one of them.

 

[EnumaElish]

The general formula for generating a random variable distributed F is y = F^-1(x) where x is the uniform random variable 0 < x < 1; i.e. f_x(x) = 1. This follows from a theorem that states that for any random variable y distributed F, the random variable x = F(y) is the uniform random variable 0 < x < 1.

 

[Pere Callahan]

Maybe, I can add another question here.

What if I want to generate a multivariate Gaussian with singular covariance matrix...in this case there does not even exist a probability density function, hence no integral representation for the cumulative distribution function.

Does somebody know what kind of algorithms is used to tackle this problem?

Thanks

Pere

 

[EnumaElish]

This reduces to a vector of N independent random Gaussian variates (N being the "multi" in "multivariate"); so you can use the univariate algorithm N times.

 

[Pere Callahan]

This reduces to a vector of N independent random Gaussian variates (N being the "multi" in "multivariate"); so you can use the univariate algorithm N times.

Does "singular covariance matrix" imply "independent"...?? I don't think so, but I may be mistaken ...

 

[D H]

What if I want to generate a multivariate Gaussian with singular covariance matrix...in this case there does not even exist a probability density function, hence no integral representation for the cumulative distribution function.

A singular covariance matrix means some of your random variables aren't random. Instead, some of your parameters are either deterministic or can be expressed as a linear combination of other parameters. I will assume your NxN covariance matrix has rank M > 1. (If the rank is one, you only have one random variable. If the rank is zero, you have no random variables.) This means that you can find some subset of size M of your N parameters such that the covariance matrix of these M parameters has rank M. Call this covariance matrix S_M. This matrix by definition is positive definite so it has a (non-unique) square root, \sqrt{S_M}.

Generate an M-vector x_M of standard gaussian deviates by drawing each element from N(0,1). The M-vector x_M = \sqrt{S_M}\,u_M will have covariance S_M. Simply add the mean vector to this M-vector if the selected variables have non-zero mean.

What about the remaining N-M parameters? They are completely determined by the already computed M random variables.

 

[EnumaElish]

Does "singular covariance matrix" imply "independent"...?? I don't think so, but I may be mistaken ...

You are right; I interpreted (read) "singular covariance matrix" as "diagonal covariance matrix" (Cov[X,Y] = 0 for X \ne Y).