Generating Random Numbers

Random numbers are important resources for scientific applications, education, game development and visualization.

The standard RTL function random generates random numbers that fulfill a uniform distribution. Uniformly distributed random numbers are not useful for every application. In order to create random numbers of other distributions special algorithms are necessary.

Normal (Gaussian) Distribution

One of the more common algorithms to produce normally distributed random numbers from uniformly distributed random numbers is the Box-Müller approach. The following function calculates Gaussian-distributed random numbers:

 function rnorm (mean, sd: real): real;
 {Calculates Gaussian random numbers according to the Box-Müller approach}
  var
   u1, u2: real;
 begin
   u1 := random;
   u2 := random;
   rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
  end;

The same algorithm is used by the randg randg function from the RTL math unit:

function randg(mean,stddev: float): float;

Exponential Distribution

An exponential distribution occurs frequently in real-world problems. A classical example is the distribution of waiting times between independent Poisson-random events, e.g. the radioactive decay of nuclei [Press et al. 1989].

The following function delivers a single real random number out of an exponential distribution. Rate is the inverse of the mean and the constant RESOLUTION determines the granularity of generated random numbers.

  function randomExp(rate: real): real;
  const
    RESOLUTION = 1000;
  var
    unif: real;
  begin
    if rate = 0 then
      randomExp := NaN
    else
    begin
      repeat
        unif := random(RESOLUTION) / RESOLUTION;
      until unif <> 0;
      randomExp := -ln(unif / rate) / rate;
    end;
  end;

Gamma Distribution

The gamma distribution is a two-parameter family of continuous random distributions. It is a generalization of both the exponential distribution and the Erlang distribution. Possible applications of the gamma distribution include modelling and simulation of waiting lines, or queues, and actuarial science.

The following function delivers a single real random number out of a gamma distribution. The shape of the distribution is defined by the parameters a, b and c. The function makes use of the function randomExp as defined above.

  function randomGamma(a, b, c: real): real;
  const
    RESOLUTION = 1000;
    T = 4.5;
    D = 1 + ln(T);
  var
    unif: real;
    A2, B2, C2, Q, p, y: real;
    p1, p2, v, w, z: real;
    found: boolean;
  begin
    A2 := 1 / sqrt(2 * c - 1);
    B2 := c - ln(4);
    Q := c + 1 / A2;
    C2 := 1 + c / exp(1);
    found := false;
    if c < 1 then
    begin
      repeat
        repeat
          unif := random(RESOLUTION) / RESOLUTION;
        until unif > 0;
        p := C2 * unif;
        if p > 1 then
        begin
          y := -ln((C2 - p) / c);
          if unif <= power(y, c - 1) then
          begin
            randomGamma := a + b * y;
            found := true;
          end;
        end
        else
        begin
          y := power(p, 1 / c);
          if unif <= exp(-y) then begin
            randomGamma := a + b * y;
            found := true;
          end;
        end;
      until found;
    end
    else if c = 1 then
    { Gamma distribution becomes exponential distribution, if c = 1 }
    begin
      randomGamma := randomExp(a, 1/b);
    end
    else
    begin
      repeat
        repeat
          p1 := random(RESOLUTION) / RESOLUTION;
        until p1 > 0;
        repeat
          p2 := random(RESOLUTION) / RESOLUTION;
        until p2 > 0;
        v := A2 * ln(p1 / (1 - p1));
        y := c * exp(v);
        z := p1 * p1 * p2;
        w := B2 + Q * v - y;
        if (w + D - T * z >= 0) or (w >= ln(z)) then
        begin
          randomGamma := a + b * y;
          found := true;
        end;
      until found;
    end;
  end;

See also

• Dev random
• Functions for descriptive statistics

References

1. G. E. P. Box and Mervin E. Muller, A Note on the Generation of Random Normal Deviates, The Annals of Mathematical Statistics (1958), Vol. 29, No. 2 pp. 610–611
2. Dietrich, J. W. (2002). Der Hypophysen-Schilddrüsen-Regelkreis. Berlin, Germany: Logos-Verlag Berlin. ISBN 978-3-89722-850-4. OCLC 50451543.
3. Press, W. H., B. P. Flannery, S. A. Teukolsky, W. T. Vetterling (1989). Numerical Recipes in Pascal. The Art of Scientific Computing, Cambridge University Press, ISBN 0-521-37516-9.