Random Gaussian Variables

Question

Is there a class in the standard library of .NET that gives me the functionality to create random variables that follow Gaussian distribution?

Answer 1

Jarrett's suggestion of using a Box-Muller transform is good for a quick-and-dirty solution. A simple implementation:

Random rand = new Random(); //reuse this if you are generating many
double u1 = 1.0-rand.NextDouble(); //uniform(0,1] random doubles
double u2 = 1.0-rand.NextDouble();
double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
             Math.Sin(2.0 * Math.PI * u2); //random normal(0,1)
double randNormal =
             mean + stdDev * randStdNormal; //random normal(mean,stdDev^2)

Answer 2

This question appears to have moved on top of Google for .NET Gaussian generation, so I figured I'd post an answer.

I've made some extension methods for the .NET Random class, including an implementation of the Box-Muller transform. Since they're extensions, so long as the project is included (or you reference the compiled DLL), you can still do

var r = new Random();
var x = r.NextGaussian();

Hope nobody minds the shameless plug.

Sample histogram of results (a demo app for drawing this is included):

Answer 3

Math.NET provides this functionality. Here's how:

double mean = 100;
double stdDev = 10;

MathNet.Numerics.Distributions.Normal normalDist = new Normal(mean, stdDev);
double randomGaussianValue=   normalDist.Sample();

You can find documentation here: http://numerics.mathdotnet.com/api/MathNet.Numerics.Distributions/Normal.htm

Answer 4

I created a request for such a feature on Microsoft Connect. If this is something you're looking for, please vote for it and increase its visibility.

https://connect.microsoft.com/VisualStudio/feedback/details/634346/guassian-normal-distribution-random-numbers

This feature is included in the Java SDK. Its implementation is available as part of the documentation and is easily ported to C# or other .NET languages.

If you're looking for pure speed, then the Zigorat Algorithm is generally recognised as the fastest approach.

I'm not an expert on this topic though -- I came across the need for this while implementing a particle filter for my RoboCup 3D simulated robotic soccer library and was surprised when this wasn't included in the framework.

---

In the meanwhile, here's a wrapper for Random that provides an efficient implementation of the Box Muller polar method:

public sealed class GaussianRandom
{
    private bool _hasDeviate;
    private double _storedDeviate;
    private readonly Random _random;

    public GaussianRandom(Random random = null)
    {
        _random = random ?? new Random();
    }

    /// <summary>
    /// Obtains normally (Gaussian) distributed random numbers, using the Box-Muller
    /// transformation.  This transformation takes two uniformly distributed deviates
    /// within the unit circle, and transforms them into two independently
    /// distributed normal deviates.
    /// </summary>
    /// <param name="mu">The mean of the distribution.  Default is zero.</param>
    /// <param name="sigma">The standard deviation of the distribution.  Default is one.</param>
    /// <returns></returns>
    public double NextGaussian(double mu = 0, double sigma = 1)
    {
        if (sigma <= 0)
            throw new ArgumentOutOfRangeException("sigma", "Must be greater than zero.");

        if (_hasDeviate)
        {
            _hasDeviate = false;
            return _storedDeviate*sigma + mu;
        }

        double v1, v2, rSquared;
        do
        {
            // two random values between -1.0 and 1.0
            v1 = 2*_random.NextDouble() - 1;
            v2 = 2*_random.NextDouble() - 1;
            rSquared = v1*v1 + v2*v2;
            // ensure within the unit circle
        } while (rSquared >= 1 || rSquared == 0);

        // calculate polar tranformation for each deviate
        var polar = Math.Sqrt(-2*Math.Log(rSquared)/rSquared);
        // store first deviate
        _storedDeviate = v2*polar;
        _hasDeviate = true;
        // return second deviate
        return v1*polar*sigma + mu;
    }
}

Answer 5

Here is another quick and dirty solution for generating random variables that are normal distributed. It draws some random point (x,y) and checks if this point lies under the curve of your probability density function, otherwise repeat.

Bonus: You can generate random variables for any other distribution (e.g. the exponential distribution or poisson distribution) just by replacing the density function.

    static Random _rand = new Random();

    public static double Draw()
    {
        while (true)
        {
            // Get random values from interval [0,1]
            var x = _rand.NextDouble(); 
            var y = _rand.NextDouble(); 

            // Is the point (x,y) below the graph of the density function?
            if (y < f(x)) 
                return x;
        }
    }

    // Probability density function of the normal "Gaussian" distribution
    public static double f(double x, double μ = 0.5, double σ = 0.5)
    {
        return 1d / Math.Sqrt(2 * σ * σ * Math.PI) * Math.Exp(-((x - μ) * (x - μ)) / (2 * σ * σ));
    }

Important: Select the interval of y and the parameters σ and μ so that the curve of the function is not cutoff at it's maximum/minimum points (e.g. at x=mean). Think of the intervals of x and y as a bounding box, in which the curve must fit in.

Answer 6

Math.NET Iridium also claims to implement "non-uniform random generators (normal, poisson, binomial, ...)".

Answer 7

Expanding off of @Noakes and @Hameer's answers, I have also implemented a 'Gaussian' class, but to simplify memory space, I made it a child of the Random class so that you can also call the basic Next(), NextDouble(), etc from the Gaussian class as well without having to create an additional Random object to handle it. I also eliminated the _available, and _nextgauss global class properties, as I didn't see them as necessary since this class is instance based, it should be thread-safe, if you give each thread its own Gaussian object. I also moved all of the run-time allocated variables out of the function and made them class properties, this will reduce the number of calls to the memory manager since the 4 doubles should theoretically never be de-allocated until the object is destroyed.

public class Gaussian : Random
{

    private double u1;
    private double u2;
    private double temp1;
    private double temp2;

    public Gaussian(int seed):base(seed)
    {
    }

    public Gaussian() : base()
    {
    }

    /// <summary>
    /// Obtains normally (Gaussian) distrubuted random numbers, using the Box-Muller
    /// transformation.  This transformation takes two uniformly distributed deviates
    /// within the unit circle, and transforms them into two independently distributed normal deviates.
    /// </summary>
    /// <param name="mu">The mean of the distribution.  Default is zero</param>
    /// <param name="sigma">The standard deviation of the distribution.  Default is one.</param>
    /// <returns></returns>

    public double RandomGauss(double mu = 0, double sigma = 1)
    {
        if (sigma <= 0)
            throw new ArgumentOutOfRangeException("sigma", "Must be greater than zero.");

        u1 = base.NextDouble();
        u2 = base.NextDouble();
        temp1 = Math.Sqrt(-2 * Math.Log(u1));
        temp2 = 2 * Math.PI * u2;

        return mu + sigma*(temp1 * Math.Cos(temp2));
    }
}

Answer 8

The other implementations are a little too heady for me. My implementation relies on the fact that throwing three six sided dice (D6) and recording the sum of the dice will creates a normal distribution. Just replace the three D6 with three calls to your your language's "rand()" function and you're just about there. In the implementation below, I subtract by two divided by the number of "dice," so 1.5 in this case.

using System;

public class CustomMath
{
    private static readonly Random _random = new Random();
    public static double GaussianRandom() =>
        _random.NextDouble() + _random.NextDouble() + _random.NextDouble() - 1.5;
}

Answer 9

I'd like to expand upon @yoyoyoyosef's answer by making it even faster, and writing a wrapper class. The overhead incurred may not mean twice as fast, but I think it should be almost twice as fast. It isn't thread-safe, though.

public class Gaussian
{
     private bool _available;
     private double _nextGauss;
     private Random _rng;

     public Gaussian()
     {
         _rng = new Random();
     }

     public double RandomGauss()
     {
        if (_available)
        {
            _available = false;
            return _nextGauss;
        }

        double u1 = _rng.NextDouble();
        double u2 = _rng.NextDouble();
        double temp1 = Math.Sqrt(-2.0*Math.Log(u1));
        double temp2 = 2.0*Math.PI*u2;

        _nextGauss = temp1 * Math.Sin(temp2);
        _available = true;
        return temp1*Math.Cos(temp2);
     }

    public double RandomGauss(double mu, double sigma)
    {
        return mu + sigma*RandomGauss();
    }

    public double RandomGauss(double sigma)
    {
        return sigma*RandomGauss();
    }
}

Answer 10

Expanding on Drew Noakes's answer, if you need better performance than Box-Muller (around 50-75% faster), Colin Green has shared an implementation of the Ziggurat algorithm in C#, which you can find here:

http://heliosphan.org/zigguratalgorithm/zigguratalgorithm.html

Ziggurat uses a lookup table to handle values that fall sufficiently far from the curve, which it will quickly accept or reject. Around 2.5% of the time, it has to do further calculations to determine which side of the curve a number is on.

Answer 11

This is my simple Box Muller inspired implementation. You can increase the resolution to fit your needs. Although this works great for me, this is a limited range approximation, so keep in mind the tails are closed and finite, but certainly you can expand them as needed.

    //
    // by Dan
    // islandTraderFX
    // copyright 2015
    // Siesta Key, FL
    //    
// 0.0  3231 ********************************
// 0.1  1981 *******************
// 0.2  1411 **************
// 0.3  1048 **********
// 0.4  810 ********
// 0.5  573 *****
// 0.6  464 ****
// 0.7  262 **
// 0.8  161 *
// 0.9  59 
//Total: 10000

double g()
{
   double res = 1000000;
   return random.Next(0, (int)(res * random.NextDouble()) + 1) / res;
}

public static class RandomProvider
{
   public static int seed = Environment.TickCount;

   private static ThreadLocal<Random> randomWrapper = new ThreadLocal<Random>(() =>
       new Random(Interlocked.Increment(ref seed))
   );

   public static Random GetThreadRandom()
   {
       return randomWrapper.Value;
   }
} 

Answer 12

You could try Infer.NET. It's not commercial licensed yet though. Here is there link

It is a probabilistic framework for .NET developed my Microsoft research. They have .NET types for distributions of Bernoulli, Beta, Gamma, Gaussian, Poisson, and probably some more I left out.

It may accomplish what you want. Thanks.