std.random.normal(), std.random.fastNormal()

[dlangBugzillaToGithub]

bearophile_hugs reported this on 2011-04-27T12:15:44Z

Transfered from https://issues.dlang.org/show_bug.cgi?id=5901

CC List

• greeenify
• lt.infiltrator

Description

I suggest to add a very commonly useful random generator with normal distribution to std.random. This is a possible API and some implementation ideas (not tested much but there is debug code that performs a visual test. This doesn't replace more rigorous tests but it's useful as sanity test for them):

// for std.random

import std.math: sqrt, log, cos, PI, isnan;
import std.traits: isFloatingPoint, CommonType;
import std.random: rndGen, Xorshift;


/**
Generates a number with normal distribution,
with specified mean, standard deviation and random generator
(mean and standard deviation must be floating point values).
*/
CommonType!(T1,T2) normal(T1, T2, UniformRandomNumberGenerator)
(T1 mean, T2 stdDev, ref UniformRandomNumberGenerator urng)
    if (isFloatingPoint!T1 & isFloatingPoint!T2) {
        alias typeof(return) T;
        enum T fmax = cast(T)(typeof(urng.front()).max);
        T r1 = void;
        do {
            r1 = urng.front() / fmax;
            urng.popFront();
        } while (r1 <= 0.0);
        T r2 = urng.front() / fmax;
        urng.popFront();
        return mean + stdDev * sqrt(-2 * log(r1)) * cos(2 * PI * r2);
    }

/// As above, but uses the default generator rndGen.
CommonType!(T1,T2) normal(T1, T2)(T1 mean, T2 stdDev)
    if (isFloatingPoint!T1 & isFloatingPoint!T2) {
        return normal(mean, stdDev, rndGen);
    }

/// As above, but uses a less precise algorithm, and a fast random generator.
CommonType!(T1,T2) fastNormal(T1, T2)(T1 mean, T2 stdDev)
    if (isFloatingPoint!T1 & isFloatingPoint!T2) {
        // this is meant to be as fast as possible ...
        // uses Xorshift too
    }


// There is a faster algorithm to compute normal values,
// for fastNormal() too, but it requires a static variable.
// For a normalDistribution() it doesn't need a static variable

/*
// 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)
// Code adapted from Python random module.
alias typeof(return) T;
enum T fmax = cast(T)(typeof(urng.front()).max);
static double gauss_next; // nan
auto z = gauss_next;
gauss_next = double.init; // nan
if (isnan(z)) {
    T x2pi = (urng.front() / fmax) * PI * 2;
    urng.popFront();
    T g2rad = sqrt(-2.0 * log(1.0 - (urng.front() / fmax)));
    urng.popFront();
    z = cos(x2pi) * g2rad;
    gauss_next = sin(x2pi) * g2rad;
}
return mu + z * sigma;
*/


debug { // Debug of normal()
    import std.stdio;

    void main() {
        {
            writeln("Debug of normal(,):");
            auto bins = new int[30];

            foreach (i; 0 .. 500) {
                auto r = cast(int)normal(bins.length / 2.0f, 5.0f);
                if (r < 0)
                    r = 0;
                if (r > (bins.length - 1))
                    r = bins.length - 1;
                bins[r]++;
            }

            foreach (count; bins) {
                write(">");
                foreach (i; 0 .. count)
                    write("*");
                writeln();
            }
            writeln("

");
        }

        {
            writeln("Debug of normal(,,Xorshift):");
            auto gen = Xorshift(1);
            auto bins = new int[30];

            foreach (i; 0 .. 500) {
                auto r = cast(int)normal(bins.length / 2.0L, 5.0L, gen);
                if (r < 0)
                    r = 0;
                if (r > (bins.length - 1))
                    r = bins.length - 1;
                bins[r]++;
            }

            foreach (count; bins) {
                write(">");
                foreach (i; 0 .. count)
                    write("*");
                writeln();
            }
            writeln();
        }
    }
}

[dlangBugzillaToGithub]

lt.infiltrator commented on 2014-03-19T17:48:25Z

Could you generalise this to a gaussian function which takes a, b, c, and d and make a PR?