Случайные гауссовские переменные

Question

Есть ли в стандартной библиотеке .NET класс, который дает мне функциональность для создания случайных величин, следующих за распределением по Гауссу?

Answer 1

Предложение Джарретта об использовании преобразования Бокса-Мюллера хорошо для быстрого и грязного решения. Простая реализация:

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

Этот вопрос, похоже, переместился на вершину Google для поколения .NET Gaussian, поэтому я решил опубликовать ответ.

Я создал несколько методов расширения для класса .NET Random , включая реализацию преобразования Бокса-Мюллера. Поскольку они являются расширениями, при условии, что проект включен (или вы ссылаетесь на скомпилированную DLL), вы все равно можете сделать

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

Надеюсь, никто не возражает против бесстыдной вилки.

Пример гистограммы результатов (демонстрационное приложение для рисования включено):

Answer 3

Math.NET обеспечивает эту функцию. Вот как:

double mean = 100;
double stdDev = 10;

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

Вы можете найти документацию здесь: http://numerics.mathdotnet.com/api/MathNet.Numerics.Distributions/Normal.htm

Answer 4

Я создал запрос на такую ​​функцию в Microsoft Connect. Если это то, что вы ищете, пожалуйста, проголосуйте за это и увеличьте его наглядность.

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

Эта функция включена в Java SDK. Его реализация доступна как часть документации и легко переносится на C # или другие языки .NET.

Если вы ищете чистую скорость, то алгоритм 1010 * Zigorat обычно считается самым быстрым.

Хотя я не эксперт по этой теме - я столкнулся с необходимостью в этом при реализации фильтра частиц для моей RoboCup 3D-моделируемой библиотеки футбольного робота и был удивлен когда это не было включено в рамки.

---

А пока вот оболочка для Random, которая обеспечивает эффективную реализацию полярного метода Бокса Мюллера:

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

Math.NET Iridium также заявляет о реализации "неоднородных случайных генераторов (нормальных, пуассоновых, биномиальных, ...)".

Answer 6

Вот еще одно быстрое и грязное решение для генерации случайных переменных, которые нормально распределены . Он рисует некоторую случайную точку (x, y) и проверяет, лежит ли эта точка под кривой вашей функции плотности вероятности, в противном случае повторите.

Бонус: вы можете генерировать случайные переменные для любого другого распределения (например, экспоненциальное распределение или распределение Пуассона ), просто заменив функцию плотности.

    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) under the curve of the density function?
            if (y < f(x)) 
                return x;
        }
    }

    // Normal (or gauss) distribution function
    public static double f(double x, double μ = 0.5, double σ = 0.5)
    {
        return 1d / Math.Sqrt(2 * σ * σ * Math.PI) * Math.Exp(-((x - μ) * (x - μ)) / (2 * σ * σ));
    }

Важно: выберите интервал y и параметры σ и μ , чтобы кривая функции не была обрезана в максимальных / минимальных точках (например, при х = среднее). Представьте себе интервалы x и y как ограничивающую рамку, в которую должна вписываться кривая.

Answer 7

Я бы хотел расширить ответ @ yoyoyoyosef, сделав его еще быстрее и написав класс-оболочку. Издержки могут означать не в два раза быстрее, но я думаю, что они должны быть в почти раза выше. Это не потокобезопасный, хотя.

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 8

Расширяя ответы @Noakes и @ Hameer, я также реализовал класс 'Gaussian', но чтобы упростить пространство памяти, я сделал его дочерним по отношению к классу Random, чтобы вы также могли вызывать базовый Next (), NextDouble () и т.д. из гауссовского класса, а также без необходимости создания дополнительного объекта Random для его обработки. Я также исключил свойства глобального класса _available и _nextgauss, поскольку не считал их необходимыми, поскольку этот класс основан на экземплярах, поэтому он должен быть потокобезопасным, если вы предоставляете каждому потоку свой гауссов объект. Я также переместил все переменные, выделенные во время выполнения, из функции и сделал их свойствами класса, это уменьшит количество обращений к диспетчеру памяти, поскольку теоретически 4 двойных числа никогда не должны перераспределяться, пока объект не будет уничтожен.

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 9

В продолжение ответа Дрю Ноакса, если вам нужна лучшая производительность, чем у Бокса-Мюллера (примерно на 50-75% быстрее), Колин Грин поделился реализацией алгоритма Зиккурата в C #, который вы можете найти здесь:

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

Ziggurat использует справочную таблицу для обработки значений, которые находятся достаточно далеко от кривой, которые он быстро примет или отклонит. Примерно в 2,5% случаев приходится выполнять дальнейшие вычисления, чтобы определить, на какой стороне кривой находится число.

Answer 10

Это моя простая реализация, вдохновленная Box Muller. Вы можете увеличить разрешение в соответствии с вашими потребностями. Хотя это прекрасно работает для меня, это приближение в ограниченном диапазоне, поэтому имейте в виду, что хвосты замкнуты и конечны, но, конечно, вы можете расширить их по мере необходимости.

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