Реализация генератора случайных чисел Бокса-Мюллера в 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
• fac становится ~ 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, но гауссова случайная переменная (которая генерирует алгоритм Бокса-Мюллера) может находиться где угодно на реальной линии.Подробнее см. wiki / 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;
    }
}