(315) |

Suppose that we wish to construct a random variable which is uniformly
distributed in the range to . In other words, the probability
density of is

(316) |

x = x1 + (x2 - x1) * double (random ()) / double (RANDMAX);

There are two basic methods of constructing non-uniformly distributed random variables:
*i.e.*, the *transformation method* and the *rejection method*. We shall examine
each of these methods in turn.

Let us first consider the transformation method. Let , where is a known function, and
is a random variable. Suppose that
the probability density of is . What is the probability density, , of ?
Our basic rule is the conservation of probability:

(317) |

(318) |

For example, consider the *Poisson distribution*:

(319) |

(320) |

(321) |

x = double (random ()) / double (RANDMAX); y = - log (x);then

The transformation method requires a differentiable probability distribution function. This is not always practical. In such cases, we can use the rejection method instead.

Suppose that we desire a random variable distributed with density in the range to . Let be the maximum value of in this range (see Fig. 95). The rejection method is as follows. The variable is sampled randomly in the range to . For each value of we first evaluate . We next generate a random number which is uniformly distributed in the range 0 to . Finally, if then we reject the value; otherwise, we keep it. If this prescription is followed then will be distributed according to .

As an example, consider the Gaussian distribution:

(322) |

(323) | |||

(324) |

since there is a negligible chance that lies more than 4 standard deviations from its mean value. It follows that

(325) |

// gaussian.cpp // Function to return random variable distributed // according to Gaussian distribution with mean mean // and standard deviation sigma. #define RANDMAX 2147483646 int random (int = 0); double gaussian (double mean, double sigma) { double ymin = mean - 4. * sigma; double ymax = mean + 4. * sigma; double Pymax = 1. / sqrt (2. * M_PI) / sigma; // Calculate random value uniformly distributed // in range ymin to ymax double y = ymin + (ymax - ymin) * double (random ()) / double (RANDMAX); // Calculate Py double Py = exp (- (y - mean) * (y - mean) / 2. / sigma / sigma) / sqrt (2. * M_PI) / sigma; // Calculate random value uniformly distributed in range 0 to Pymax double x = Pymax * double (random ()) / double (RANDMAX); // If x > Py reject value and recalculate if (x > Py) return gaussian (mean, sigma); else return y; }Figure 96 illustrates the performance of the above function. It can be seen that the function successfully returns a random value distributed according to the Gaussian distribution.