Прямоугольник?

(МЗ. Клин или хвост?) >

Хвост

j;;yib. Простой случай?)->

ДаМб. Попытаемся снова?У

М7. Получить мажорирую- М8. Отбраковать?Л

щую случайную величину Ч, \" /

--1 уДа

М9. Присоединить знак

Т"

Рис. 12. Алгоритм прямоугольника-клина-хвоста для генерирования нормальных случайных величин.

Сейчас вероятность того, что X <s + hx для О < х < 1, - это площадь, лежащая слева от вертикальной линии U = х (см. рис. 11) и деленная на общую площадь, т. е.

f{s + hu)du 11 f{s + hu)du =

f{v)dv.

Следовательно, X имеет нужное распределение.

С подходящими константами uj , bj и Sj алгоритм L можно применять к клинообразным плотностям /,+15 (см. рис. 9) для 1 < j < 15. Последнее распределение F31 нуждается в обработке приблизительно один раз из 370; оно используется тогда, когда получаем результат АГ > 3. В упр. И показано, что стандартная схема отбраковки может использоваться для этого "хвоста". Рассмотрим процедуру в целом.

Алгоритм М (Метод прямоугольника-клина-хвоста для нормально распределенных случайных величин). Для этого алгоритма (рис. 12) используем вспомогательные таблицы {Ро,...,Рз1), (Qi,...,Qi5), {Yo,...,Y3i), {Zo,...,Z3i), (5i,... ,Si6), (J9i6,... ,1>зо), (Eie,... ,Езо), построенные так, как в упр. 10; примеры выбираются из табл. 1. Предполагается, что используется бинарный компьютер (подобные процедуры могут быть построены для десятичных машин).

Ml. [Получить и.] Генерируем равномерно распределенное случайное число U = (.606162.. .6()2. (Здесь bj - двоичные разряды в бинарном представлении U. Для приемлемой точности t должно быть по крайней мере равно 24.) Присвоить ф -(г- bo. (Позже ф будет использовано для определения знака результата.)

М2. [Прямоугольник?] Присвоить j (6162636465)2 (двоичное число определяется старшими двоичными разрядами U) и присвоить / i- (.бебу... bt)2 (дробная часть определяется оставшимися двоичными разрядами). Если f > Pj, присвоить X (г- Yj + fZj и перейти к шагу М9. Иначе, если j < 15 (т. е. 61 =0), присвоить X <г- Sj + fQj и перейти к шагу М9. (Это использование метода псевдонимов Уолкера (Walker) (3).)

МЗ. [Клин или хвост?] (Сейчас 15 < j < 31 и каждое определенное значение j появляется с вероятностью pj.) Если j = 31, перейти к шагу М7.

М4. [Получить и < V.] Генерировать две новые независимые равномерно распределенные случайные величины U иУ; если U > V, заменить U V. (Сейчас выполняется алгоритм L.) Присвоить X i- Sj-15 + \U.

М5. [Простой случай?] Если V < Dj, перейти к шагу М9.

Мб. [Попытаемся снова?] Если V > U -\- Ej{ef-*~-> - 1), вернуться к шагу М4; иначе - перейти к шагу М9. (Вероятность, что этот шаг нужно будет выполнять, мала.)

М7. [Получить мажорирующую стучайную величину.] Генерировать две новые независимые равномерно распределенные случайные величины (7 и F и присвоить X <г- v/9-21nl/.

М8. [Отбраковать?] Если UX >Ъ, возвратиться к шагу М7. (Это будет случаться приблизительно в одном из двенадцати случаев после достижения шага MS.)

М9. [Присоединить знак.] Если ф = 1, присвоить Л (--X

Этот алгоритм является очень миленьким примером тесного переплетения математической теории с изобретательностью программирования - прекрасная иллюстрация искусства программирования! На выполнение шагов Ml, М2 и М9 уходит почти все время; остальные шаги осуществляются намного быстрее. Впервые метод прямоугольника-клина-хвоста опубликовали Дж. Марсалья, Дж. Марсалья, М. Д. Мак-Ларен и Т. А. Брей (см. работы G. Marsaglia, Annals Math. Stat. 32 (1961), 894-899; G. Marsaglia, M. D. MacLaren, and T. A. Bray, CACM 7 (1964), 4-10). В дальнейшем алгоритм М усовершенствовали Дж. Марсалья, К. Ананханараянан и Н. Дж. Поль (см. работу G. Marsaglia, К. Ananthanarayanan, and N. J. Paul, Inf. Proc. Letters 5 (1976), 27-30).

3) Метод "чет-нечет" принадлежит Дж. Э. Форсайту (G. Е. Forsythe). Поразительно простая техника генерирования случайных величин с обычной экспоненциальной плотностью вида

f{x) = Ce-[a<x<b], (20)

О < Н{х) < 1 для а < ж < 6, (21)

была открыта Джоном фон Нейманом и Дж. Е. Форсайтом приблизительно в 1950 году. Идея основана на методе отбраковки, описанном ранее. Предпсшага-ется, что д{х) - это равномерное распределение на интервале [о..6). Присвоим X а -\- (Ь - a)U, где и - равномерно распределенная случайная величина, и примем X с вероятностью е"*-). Последняя операция может быть осуществлена путем сравнения e~ с V либо h{X) с -InV, когда V - другая равномерно

Таблица 1

ПРИМЕРЫ ТАБЛИЦ, ИСПОЛЬЗУЕМЫХ С АЛГОРИТМОМ М*

*Практически эти данные могут быть приведены с намного большей точностью; в таблице приведено только достаточное количество цифр, чтобы интересующиеся читатели могли более точно проверить собственные алгоритмы вычисления значений.

распределенная случайная величина, но работу можно проделать без применения каких-либо трансцендентных функций следующим интересным способом. Присвоим Vq i- h{x) и будем генерировать обычные равномерно распределенные случайные величины Vi, V2, пока не найдутся К > I с Vk-i < V. Для фиксированных x и к вероятность того, что h{x) > Vi > • • • > V, равна 1/к\, умноженному на вероятность того, что ma.x{Vi,... ,Vk) < h{x), т. е. h{x)/k\. Следовательно, вероятность того, что К = к, равна h{x)~/{k - 1)! - h{x)/k\, а вероятность того, что К - нечетное число, равна

нечетное, fc>l

Следовательно, мы отбраковываем x и начинаем снова, если К - четное число. Принимаем x как случайную величину с плотностью распределения (20), если К - нечетное число. Обычно мы не генерируем много значений Vj, чтобы определить К, поскольку среднее значение К (при заданном x) равно k>o Рг(- > к) = E,>oWA! = eW<e.

Несколькими годами позже Форсайт показал, что этот подход приводит к получению эффективного метода вычисления нормальных случайных величин без использования вспомогательных программ для вычисления квадратных корней или логарифмов, как в алгоритмах Р и М. Его процедуру с улучшенным выбором интервалов [а..Ь), осуществленную И. Г. Аренсом (J. Н. Ahrens) и У. Дитером (и. Dieter), можно представить в виде алгоэитма.

Алгоритм F (Метод "чет-нечет" для нормальных случайных величин). Этот алгоритм генерирует нормальные случайные величины на бинарном компьютере; предполагается, что имеется приблизительно t-\-l двоичных разрядов точности. Он нуждается в таблице значений dj - Oj - Oj-i для I < j <t+l, где Oj определяется