Что плохого в моей попытке сгенерировать случайное число Бокса-Мюллера? - PullRequest
Новая
       41

Что плохого в моей попытке сгенерировать случайное число Бокса-Мюллера?

0 голосов
/ 04 апреля 2019

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

Можете ли вы сказать мне, что я делаю не так?

.stat.GetOneGaussianByBoxMuller:{ 
sq:{
    a:2.0*rand[abs[system"S"]]%abs[system"S"]-1;
    b:2.0*rand[abs[system"S"]]%abs[system"S"]-1;
    sq:(a*a)+(b*b);
    x:sq
   }/[{x>=1};1]; 

  :(2.0*rand[abs[system"S"]]%abs[system"S"]-1)*sqrt[(neg[2]*log[sq])%sq]
 };

1 Ответ

3 голосов
/ 04 апреля 2019

В вашем ответе есть три вопроса: приоритет операций в q, как генерировать случайную переменную в [-1; 1] и условие для переменной sq.

Приоритет операций:

Напомним, что ниже a1 и a2 различны: 

u:rand[1f];
a1:-1+2.0*u;
a2:2.0*u-1;

u:0.5
a1:0f
a2:-1f

u:0
a1:-1f
a2:-2f

Генерация случайной величины в [-1; 1]:

Используйте следующее: 

a:-1+2.0*rand[1.0];
b:-1+2.0*rand[1.0];

Состояние на sq:

Вы действительно приняли во внимание условие sq >= 1, но у вас также есть проблема, если sq = 0, поскольку вы должны делить на sq на последнем шаге.

Кроме того, в вашей реализации вы вычисляете a дважды, что 1) не оптимально и 2) не соответствует методологии, поскольку тот же самый a должен использоваться при вычислении sq в последнем шаг, последний уступает очень большим числам. Я получил некоторое вдохновение со страницы википедии , где они предлагают восстановить sq, если вышеуказанное условие не выполнено. Отсюда рекурсивный вызов функции в приведенной ниже реализации: 

.stat.GetOneGaussianByBoxMuller:{ 
    a:-1+2.0*rand[1.0];
    b:-1+2.0*rand[1.0];
    sq:(a*a)+(b*b);
    if[(sq>=1)|(sq=0);
        :.stat.GetOneGaussianByBoxMuller[];
    ];
    :a*sqrt[(neg[2]*log[sq])%sq];
 };

Теперь вы можете взглянуть на данные, сгенерированные путем построения гистограммы следующего набора данных: 

([]val:{:.stat.GetOneGaussianByBoxMuller[]} each til 100000)

EDIT:

На самом деле вы могли бы иметь более эффективную реализацию, генерируя 2 случайных числа для вызова функции следующим образом: 

.stat.GetOneGaussianByBoxMuller:{ 
    a:-1+2.0*rand[1.0];
    b:-1+2.0*rand[1.0];
    sq:(a*a)+(b*b);
    if[(sq>=1)|(sq=0);
        :.stat.GetOneGaussianByBoxMuller[];
    ];
    :(a*sqrt[(neg[2]*log[sq])%sq];b*sqrt[(neg[2]*log[sq])%sq]);
 };

([]val:raze {.stat.GetOneGaussianByBoxMuller[]} each til 50000)

Эта реализация требует 150 мс для генерации 100000 нормально распределенных случайных чисел, в то время как приведенная выше занимает 245 мс.

...