Issue 5901 - std.random.normal(), std.random.fastNormal()
Summary: std.random.normal(), std.random.fastNormal()
Status: NEW
Alias: None
Product: D
Classification: Unclassified
Component: phobos (show other issues)
Version: D2
Hardware: All All
: P4 enhancement
Assignee: Seb
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2011-04-27 12:15 UTC by bearophile_hugs
Modified: 2024-12-01 16:14 UTC (History)
2 users (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this issue.
Description bearophile_hugs 2011-04-27 12:15:44 UTC
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("\n\n");
        }

        {
            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();
        }
    }
}
Comment 1 Infiltrator 2014-03-19 17:48:25 UTC
Could you generalise this to a gaussian function which takes a, b, c, and d and make a PR?
Comment 2 dlangBugzillaToGithub 2024-12-01 16:14:06 UTC
THIS ISSUE HAS BEEN MOVED TO GITHUB

https://github.com/dlang/phobos/issues/9903

DO NOT COMMENT HERE ANYMORE, NOBODY WILL SEE IT, THIS ISSUE HAS BEEN MOVED TO GITHUB