отношением

Ге-.=1. ,.3)

F1. [Получить и.] Генерировать равномерно распределенное число U-{.ЬоЬх... 6)2, где 6о, 61, bt означают двоичные разряды в бинарной записи. Присвоить ф <г-bo, j ir-I-а а <г-0.

F2. [Найти первый нуль bj] Если bj = 1, присвоить о (- a + dj, j (- j + l и повторить этот шаг. (Если j = t + I, трактовать bj как нуль.)

F3. [Генерировать X.] (Сейчас о = aj-i, и текущее значение j появится с вероятностью и 2~. Будем генерировать X в интервале [aj-i. .Uj), используя метод отбраковки, описанный выше, с h{x) = - = у/2 + ау, где у = X - а. В упр. 12 доказывается, что h{x) < 1, как требуется в (21).) Присвоить Y <- dj, умноженное на (.6+1 ...6)2 и V f- {Y + a)Y. (Так как среднее значение j равно 2, обычно достаточно старших двоичных разрядов в {.bj+i.. .bt)2 для обеспечения приличной точности. Вычисления без труда выполняются арифметикой с фиксированной точкой.)

F4. [Отбраковать?] Генерируем равномерно распределенную случайную величину и. Если V < и, перейти к шагу F5. Иначе - присвоить V новую равномерно распределенную случайную величину. Если сейчас U < V (т. е. если К четное, как обсуждалось выше), возвратиться к шагу F3, иначе - повторить шаг F4.

F5. [Выход из программы по X] Присвоить X i- а + Y. Если ф = I, присвоить X -X. I

Значения dj для 1 < j < 47 появились в статье Аренса и Дитера (Ahrens and Dieter, Math. Сотр. 27 (1973), 927-937). В ней обсуждается усовершенствованный алгоритм, который повышает скорость его выполнения за счет использования больших таблиц. Алгоритм F привлекателен тем, что работает почти с такой же скоростью, как алгоритм М, и его легче выполнить. Среднее число равномерно распределенных случайных величин на каждую нормальную случайную величину равно 2.53947. Р. П. Брент (R. Р. Brent, САСМ 17 (1974), 704-705) показал, как можно сократить это число до 1.37446 с помощью двух вычитаний и одного деления на каждую хранимую равномерно распределенную случайную величину.

4) Отношение равномерно распределенных случайных величин. Существует еще один хороший метод генерирования нормальных случайных величин, открытый в 1976 году А. Дж. Киндерманом (А. J. Kinderman) и Дж. Ф. Монаханом (J. F. Мо-nahan). Его суть - генерирование случайной точки {U, V) в области, определенной как

0<и<1, -2иу/1п{1/и) <v< 2иу/\п{1/и), (24)

и получение отношения X (- V/U. Заштрихованная площадь на рис. 13 - это магическая область неравенства (24), осуществляющая эту работу. Прежде чем изучать теорию, приведем алгоритм, эффективность и простота которого очевидны.

1 = 3

(0,y27i

(0,-y27i

1,27)

1 = 1/3

1 = 0

x = -l/3

(i,-y27i)

Рис. 13. Область "принятия" в методе отношения равномерных случайных величин для нормально распределенных случайных величин. Длины линий с отношением координат X распределены нормально.

1 = -3

1= -1

Алгоритм R [Метод отношений для нормальных случайных величин). Этот алгоритм генерирует нормальные случайные величины X.

R1. [Получить и, V.] Генерировать две независимые равномерно распределенные случайные величины U и V, где U ф О, и присвоить X +- y/S/e {V - ) /U. (Сейчас X равно отношению координат {U, \/SJe {У - )) случайной точки в прямоугольнике, содержащем заштрихованную область на рис. 13. Принимаем А, если соответствующая точка в самом деле находится в заштрихованной области, иначе - начинаем сначала.)

R2. [Необязательная проверка верхней грани.] Если Х < 5 - АеЮ, на выходе - АГ и завершение алгоритма. (Этот шаг можно опустить, если пожелаете; он так или иначе проверяет, находятся ли избранные точки во внешней области рис. 13, делая излишним вычисление логарифма.)

R3. [Необязательная проверка нижней грани.] Если Х" > 4e~-/U + 1.4, возвратиться к шагу R1. (Этот шаг также может быть опущен; он так или иначе проверяет, находятся ли выбранные точки за внешней областью рис. 13, делая ненужным вычисление логарифма).

R4. [Окончательная проверка.] Если АГ < -4 In С/, выход X и завершение алгорти-ма; иначе - возврат к шагу R1.

В упр. 20 и 21 представлен временной анализ; анализируются четыре различных алгоритма, так как шаги R2 и R3 при желании могут быть включены или пропущены. В следующей таблице показано, сколько в среднем времени уходит на

выполнение каждого шага в зависимости от применения необязательной проверки.

Шаг Ни одного Только R2 Только R3 Оба

R1 1.369 1.369 1.369 1.369

R2 О 1.369 О 1.369 (25)

R3 О О 1.369 0.467

R4 1.369 0.467 1.134 0.232

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

Но зачем делать эту работу? По одной единственной причине - чтобы вычислить вероятность того, что X < х, которая оказывается точным значением (10). Но такое вычисление вряд ли будет очень простым, если не удастся найти хороший прием. Во всяком случае, в первую очередь, нужно понять, как может быть открыт такой алгоритм. Киндер.ман (Kinderman) и Монахам (Monahan) открыли его, разрабатывая следующую теорию, которая может использоваться с любой хорошо себя ведущей плотностью f{x) [см. АСМ Trans. Math. Software 3 (1977), 257-260].

Предположим, что точка {U, V) равномерно генерируется в области {и, г;)-плос-кости, определенной неравенствами

и> О, < g{v/u) (26)

для некоторой неотрицательной интегрируемой функции д. Если присвоить X ir-V/U, то вероятность того, что X < х, можно вычислить путем интегрирования по dudv по всей области, определенной двумя соотношениями в (26) и добавочным условием v/u < х с последующим делением на этот же интеграл, но без дополнительного условия. Допустим, V = tu, так что dv = udt. Тогда интеграл имеет вид

ГуДЩ 1

dt udu = - g{t) dt.

оо Jo J-oo

Следовательно, вероятность, что X <x, равна

Г g{t)dt I Г" g{t)dt. (27)

Нормальное распределение возникает, когда д{\) - е~* 1, а условие v? < givju) упрощается в этом случае до {v/uf- < -4 In и. Легко видеть, что множество всех (u,v), удовлетворяющих этому соотношению, полностью содержится в прямоугольнике на рис. 13.

Грани шагов R2 и R3 определяют внутреннюю и внешнюю области простыми граничными уравнениями. Хорошо известное неравенство

еУХх,

справедливое для всех действительных чисел ж, можно использовать, чтобы показать, что неравенства

1-hlnc-cu < -Inu < l/(cu) - 1-hlnc (28)