Реализация генератора случайных чисел Бокса-Мюллера в C#

Question

Из этого вопроса: Генератор случайных чисел, который притягивает числа к любому заданному числу в диапазоне? Я провел некоторое исследование, так как я сталкивался с таким генератором случайных чисел прежде. Все, что я помню, было имя "Мюллер", так что, думаю, я нашел его здесь:

• Преобразование Бокса-Мюллера

Я могу найти множество реализаций этого на других языках, но я не могу реализовать это правильно в C#.

На этой странице, например, метод Бокса-Мюллера для генерации гауссовских случайных чисел, говорится, что код должен выглядеть следующим образом (это не C#):

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

double gaussian(void)
{
   static double v, fac;
   static int phase = 0;
   double S, Z, U1, U2, u;

   if (phase)
      Z = v * fac;
   else
   {
      do
      {
         U1 = (double)rand() / RAND_MAX;
         U2 = (double)rand() / RAND_MAX;

         u = 2. * U1 - 1.;
         v = 2. * U2 - 1.;
         S = u * u + v * v;
      } while (S >= 1);

      fac = sqrt (-2. * log(S) / S);
      Z = u * fac;
   }

   phase = 1 - phase;

   return Z;
}

Теперь вот моя реализация выше в C#. Обратите внимание, что преобразование производит 2 числа, отсюда трюк с "фазой" выше. Я просто отбрасываю второе значение и возвращаю первое.

public static double NextGaussianDouble(this Random r)
{
    double u, v, S;

    do
    {
        u = 2.0 * r.NextDouble() - 1.0;
        v = 2.0 * r.NextDouble() - 1.0;
        S = u * u + v * v;
    }
    while (S >= 1.0);

    double fac = Math.Sqrt(-2.0 * Math.Log(S) / S);
    return u * fac;
}

Мой вопрос связан со следующим конкретным сценарием, где мой код не возвращает значение в диапазоне 0-1, и я не могу понять, как оригинальный код тоже может.

• U = 0,5, V = 0,1
• S становится 0.5*0.5 + 0.1*0.1 знак равно 0.26
• лицо становится ~3.22
• Таким образом, возвращаемое значение ~0.5 * 3.22 или ~1.6

Это не в 0 .. 1,

Что я делаю не так / не понимаю?

Если я изменю свой код так, чтобы вместо умножения fac с uЯ умножаю на SЯ получаю значение в диапазоне от 0 до 1, но оно имеет неправильное распределение (кажется, имеет максимальное распределение около 0,7-0,8, а затем сужается в обоих направлениях.)

Answer 1

Ваш код в порядке. Ваша ошибка - думать, что он должен возвращать значения исключительно в пределах [0, 1], (Стандартное) нормальное распределение - это распределение с ненулевым весом по всей вещественной линии. То есть значения за пределами [0, 1] возможны На самом деле, ценности внутри [-1, 0] так же вероятно, как значения в [0, 1]и, кроме того, дополнение [0, 1] имеет около 66% веса нормального распределения. Таким образом, в 66% случаев мы ожидаем значение вне [0, 1],

Кроме того, я думаю, что это не преобразование Бокса-Мюллера, а полярный метод Марсальи.

Answer 2

Я не математик и не статистик, но если бы я подумал об этом, я бы не ожидал, что распределение Гаусса будет возвращать числа в точном диапазоне. Учитывая вашу реализацию, среднее значение равно 0, а стандартное отклонение равно 1, поэтому я ожидаю, что значения распределены по кривой колокола с 0 в центре, а затем уменьшаются по мере отклонения чисел от 0 с обеих сторон. Таким образом, последовательность определенно охватывает оба +/- числа.

Тогда, поскольку он статистический, почему его трудно ограничить -1..1 только потому, что std.dev равен 1? Статистически может быть какая-то игра с любой стороны и все же выполнить статистическое требование.

Answer 3

Равномерная случайная переменная действительно находится в пределах 0..1, но гауссова случайная переменная (которую генерирует алгоритм Бокса-Мюллера) может находиться где угодно на реальной линии. Смотрите вики /NormalDistribution для получения более подробной информации.

Answer 4

Я думаю, что функция возвращает полярные координаты. Таким образом, вам нужны оба значения, чтобы получить правильные результаты.

Кроме того, распределение Гаусса не между 0 .. 1, Это может легко закончиться как 1000, но вероятность такого появления чрезвычайно низка.

Answer 5

Это метод Монте-Карло, поэтому вы не можете зафиксировать результат, но вы можете игнорировать сэмплы.

// return random value in the range [0,1].
double gaussian_random()
{
    double sigma = 1.0/8.0; // or whatever works.
    while ( 1 ) {
        double z = gaussian() * sigma + 0.5;
        if (z >= 0.0 && z <= 1.0)
            return z;
    }
}