LinuxQuestions.org

LinuxQuestions.org (/questions/)
-   Programming (http://www.linuxquestions.org/questions/programming-9/)
-   -   Generating Gaussian Random Numbers (http://www.linuxquestions.org/questions/programming-9/generating-gaussian-random-numbers-215942/)

R00ts 08-10-2004 04:37 PM

Generating Gaussian Random Numbers
 
I'm a little stuck trying to write function to generate a Gaussian Random number with range cut-offs. I've google'd my ass off, believe me, and I'm close to my answer but not quite there. Here's what I got so far:

Code:

// Returns a random number between [-1, 1)
float UnitRV() {
  return 2 * ((float) rand()/RAND_MAX) - 1;


// Uses the Polar method to generate a normal distribution with mean 0 and standard deviation 1
// http://ct.radiology.uiowa.edu/~jiang...ml/node41.html
float GaussianValue(int mean, int range) {
  static int iset = 0;
  static float gset;
  float fac, rsq, v1, v2;
  float result;
  int true_range;
 
  if  (iset == 0)  {
    do {
      v1  = 2.0 * UnitRV() - 1.0;
      v2  = 2.0 * UnitRV() - 1.0;
      rsq = v1 * v1 + v2 * v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    fac = sqrt (-2.0 * log (rsq) / rsq);
    gset = v1 * fac;
    iset = 1;
    result = v2 * fac;
  }
  else {
    iset = 0;
    result = gset;
  }
 
  return result;
}

The mean and range arguments in the GaussianValue function are what need to be used. Obviously to change the mean of the normal distribution I just have to add the mean argument to result, that's easy.

The hard thing is that I want to change the standard distribution so that sigma ~= range / 3. (So range will be the boundary of the 3rd standard deviation). I couldn't find any details on how this function computes the standard deviation. Does anyone know how it work? The GaussianValue function is currently almost the exact same as the source where I got it at. Thanks for any help :study:

kev82 08-10-2004 05:40 PM

let X~N(mu,sigma) then z=((x-mu)/sigma)~N(0,1) this is basic normal dist stuff now lets do it backwards ie let z~N(0,1) then X=(sigma*z + mu)~N(mu,sigma) so there ya go. with a N(0,1) generator to make it N(mu, sigma) transform it by sigma*X + mu.

you can prove this without too much trouble by writing E[X] as integral(x*p(x)) (for continuous, obviously sum for discrete) and showing E[aX+b] = aE[X] + b and then considering variance as E[X^2] - E^2[X]

also your code is wrong, it works by writing the normal distribution integral in 2d space and picking a point(v1,v2) in the unit circle which is why you reject if v1^2 + v2^2 > 1. so v1,v2 should lie between -1 and 1 which is why the code calls UnitRV() which should return a value between 0 and 1 however your UnitRV returns a value between -1 and 1 putting v1 and v2 in the range -3 to 1 which shouldnt effect your numbers (i think) but will require more numbers to be generated per each (x1,x2)

<edit>
ive just looked at the maths for generating a N(mu,sigma) distribution closely instead of guessing(as i did originally), you have to do a translated elliptical substitution (x-a=q*cos t, y-b=r*sin t rather than the much simpler polar subsition x=r*cos t,y=r*sin t) into the double integral and its far too much work when you can just multiply by variance and add mean to the resulting N(0,1) numbers.

**i deleted the paragraph that was here cos it was a load of nonsense
</edit>

my masters project is based on random number generation so i do know quite a bit about it and have probably assumed too much knowledge on your part, if theres anything you dont get then please ask.

R00ts 08-10-2004 11:51 PM

Thanks kev, that made perfect sense. :D I really need to review my mathematics, you made it seem so easy. As for your concerns on the function for finding N(1,0), I did change it a little after reading some more (made it simpler) and I've done several test runs on it and it seemed good enough for me. If I introduce eliptical geometry all of those trigonometric calls to the math library are going to consume CPU cycles, and this function doesn't have to be entirely accurate for my purposes. Its accurate enough as it is now, but thanks for your concern. :)


All times are GMT -5. The time now is 05:05 AM.