Random Gaussian Variables

C#.Net

C# Problem Overview


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

C# Solutions


Solution 1 - C#

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)

Solution 2 - C#

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):

enter image description here

Solution 3 - C#

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

Solution 4 - C#

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;
    }
}

Solution 5 - C#

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.

Solution 6 - C#

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

Solution 7 - C#

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));
    }
}

Solution 8 - C#

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();
    }
}

Solution 9 - C#

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.

Solution 10 - C#

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;
}

Solution 11 - C#

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.

Solution 12 - C#

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;
   }
} 

Attributions

All content for this solution is sourced from the original question on Stackoverflow.

The content on this page is licensed under the Attribution-ShareAlike 4.0 International (CC BY-SA 4.0) license.

Content TypeOriginal AuthorOriginal Content on Stackoverflow
QuestionSebastian M&#252;llerView Question on Stackoverflow
Solution 1 - C#yoyoyoyosefView Answer on Stackoverflow
Solution 2 - C#SuperbestView Answer on Stackoverflow
Solution 3 - C#Gordon SlyszView Answer on Stackoverflow
Solution 4 - C#Drew NoakesView Answer on Stackoverflow
Solution 5 - C#DoomjunkyView Answer on Stackoverflow
Solution 6 - C#Jason DeFontesView Answer on Stackoverflow
Solution 7 - C#user8262209View Answer on Stackoverflow
Solution 8 - C#Hameer AbbasiView Answer on Stackoverflow
Solution 9 - C#NeilView Answer on Stackoverflow
Solution 10 - C#Andrew AllbrightView Answer on Stackoverflow
Solution 11 - C#Aaron StainbackView Answer on Stackoverflow
Solution 12 - C#Daniel HowardView Answer on Stackoverflow