Normal distributed rand

bearophile bearophileHUGS at lycos.com
Wed Apr 8 17:53:17 PDT 2009


nobody:
> Does anyone know of a normal/gaussian distributed random number generator 
> for D1 + phobos? 

I have created this for you, not much (unit)tested yet, adapted from Python random module, I'll put it in my dlibs. Maybe there are ways to do this in a faster way (especially if you just need an approximate normal distribution).
Also take a look at Tango source code, it too surely has a similar function/method.

import std.math: sin, cos, log, sqrt, PI, isnan;
import std.random: rand;

double random() {
    return rand() / cast(double)uint.max;
}

/*************************************
Normal distribution.
mu is the mean, and sigma is the standard deviation. Not thread safe.
*/
double normal(double mu=0.0, double sigma=1.0) {
    // When x and y are two variables from [0, 1), uniformly
    // distributed, then
    //
    //    cos(2*pi*x)*sqrt(-2*log(1-y))
    //    sin(2*pi*x)*sqrt(-2*log(1-y))
    //
    // are two *independent* variables with normal distribution
    // (mu = 0, sigma = 1).
    // (Lambert Meertens)
    static double gauss_next; // nan
    auto z = gauss_next;
    gauss_next = double.init; // nan
    if (isnan(z)) {
        double x2pi = random() * PI * 2;
        double g2rad = sqrt(-2.0 * log(1.0 - random()));
        z = cos(x2pi) * g2rad;
        gauss_next = sin(x2pi) * g2rad;
    }
    return mu + z * sigma;
}


//--------------------------------

import std.string: repeat;
import std.stdio: put = writef, putr = writefln;

void main() {
    auto bins = new int[30];

    for (int i; i < 10000; i++) {
        auto r = cast(int)normal(bins.length / 2, 5);
        if (r < 0)
            r = 0;
        if (r > (bins.length - 1))
            r = bins.length - 1;
        bins[r]++;
    }

    foreach (count; bins)
        putr(">", "*".repeat(count / 20));
}

Bye,
bearophile


More information about the Digitalmars-d-learn mailing list