Д. Кнут - Искусство программирования том 2 (3-е издание) - 2001 (Часть 1) (1119452), страница 40
Текст из файла (страница 40)
12. Алгоритм прямоугольника-клина-хвоста для генерирования нормальных случай- ных величин. Сейчас вероятность того, что Х < з+ Ьх для О < х < 1, — это площадь, лежащая слева от вертикальной линии У = х (см. рис. 11) и деленная на общую площадь, т. е, / ь l г'1 Гььь* -1(з+ йи) Нн/ / -у(з+ Ьи) Ип = / У(и) до. о е Следовательно, Х имеет нужное распределение. С подходящими константами а, Ьу и з. алгоритм 1 можно применять к клинооб- разным плотностям Узчлз (см. рнс. 9) для 1 < у < 15. Последнее распределение Гзг нуждается в обработке приблизительно один раз из 370; оно используется тогда, когда получаем результат Х > 3, В упр.
11 показано, что стандартная схема отбраковки может использоваться для этого "хвоста". Рассмотрим процедуру в целом. Алгоритм М (Метлод прямоугольника-клина-хеостла длл нормально распределенных случайных величин). Для этого алгоритма (рис. 12) используем вспомогательныетаблнцы(Ро,. Рз~) (9г Яы), (уо*,Уз~), (Хо,-.,Хзг) (%,,~и) Фзе, >1Узе) (Еы, . Езе), построенные так, как в упр.
10; примеры выбираются из табл. 1. Предполагается, что используется бинарный компьютер (подобные процедуры могут быть построены для десятичных машин). М1. [Получить Г) Генерируем равномерно распределенное случайное число П = ( ЬоЬ|Ьз . Ьс)з. (Здесь Ьу — двоичные разряды в бинарном представлении Г Для приемлемой точности 1 должно быть по крайней мере равно 24.) Присвоить ~4 +- Ье. (Позже 4 будет использовано для определения знака результата.) М2.
(Прямоугольник?) Присвоить у' +- (Ь«Ь«Ь«Ь«Ь«)з (двоичное число определяется старшими двоичными разрядами Ц и присвоить 1 +- (.Ь«Ьт... Ь«)з (дробная часть определяется оставшимися двоичными разрядами). Если у > Р,, присвоить Х «- У'+ ~Я. и перейти к шагу М9. Иначе, если у < 15 (т. е. Ь« = О), присвоить Х «- Я~ + Я н перейти к шагу М9. (Это использование метода псевдонимов Ъолкера (%аИ«ег) (3).) МЗ.
(Клин или хвост?) (Сейчас 15 < у < 31 и каждое определенмое значение у появляется с вероятностью р ..) Если у = 31, перейтн к шагу М7. М4. (Получить У < Г) Генерировать две новые независимые равномерно распределенные случайные величины П и Ъ', если У > Ь', заменить (7 «-> )г.
(Сейчас выполняется алюрнтм 1..) ПРисВОить Х «- 91-««+ «(? М5. (Простой случай?) Если Ь' < 17~, перейти к шагу М9. Мб. (Попытаемся снова?) Если Ъ' > У+ Еу(е«г~-и ~ !?~ — 1), вернуться к шагу М4; иначе — перейти к шагу М9. (Вероягмость, что этот шаг нужно будет выполнять, мала.) М7. (Получить мажорируюшую случайную величину.) Генерировать две новые независимые равномерно распределенные случайные величины У и Ъ'и присвоить Х 9 — 11~'.
М8. [Отбрвковать?) Если УХ > 3, возвратиться к шагу М7. (Это будет случаться приблизительно в одном из двенадцати случаев после достижения шага М8.) М9. (Присоединить знак.) Если ф = 1, присвоить Х «- -Х. 1 Этот алгоритм является очень миленьким примером тесного переплетения мате««атмческой теории с изобретательностью программирования — прекрасная иллюстрация искусства программирования! На выполнение шагов М1, М2 и М9 уходит почти все время; остальные «наги осуществляются намного быстрее. Впервые метод прямоугольника-клнна-хвоста опубликовали Дж, Марсэлья, Дж. Марсалья, М.
Д. Мак-Ларси и Т. А. Брей (см. работы 6. Мвгэай!!а, Аппа)э Ыа«Ь. 5«ак 32 (1961), 894 — 899; 6. 54ап«ай!мн М, П. Мас1 агеп, ап«( Т, А. Вгау, САСМ 7 (1964), 4-10). В дальнейшем алгоритм М усовершенствовали Дж. Марсэлья, К. Ананханараянан и Н. Дж. Поль (см. работу 6. Махваййа, К. Анап«Ьапагауапап, ап«1 Х. 3. Рап(, 1пй Ргос. Ьеггегэ 5 (1976), 27 — 30). 3) Мепюд "чет-мечеть принадлежит Дж. Э. Форсайту (6.
Е. рог«у«пе). Поразительно простая техника генерирования случайных величин с обычной экспоненциальной плотностью вида 1(х) = Се ьы!)а<х<Ь), (20) где 0<Ь(х)<1 дляа<х<Ь, (21) была открыта Джоном фон Неймамом и Дж. Е. Форсайтом приблизительно в 19о0 году. Идея основана на к«етоде отбраковки, описанном ранее. Предполагается, что д(х) — это равномерное распределение на интервале (а..Ь). Присвоим Х «- а + (Ь вЂ” а)У, где П -- равмомерно распределенная случайная величина, и примем Х с вероятностью е "«Я!.
Последняя операция может быть осуществлена путем сравнения е «~! с Ь' либо Ь(Х) с -!пЪ; когда $' — другая равномерно Таблица 1 ПРИМЕРЫ ТАБЛИЦ, ИСПОЛЬЗУЕМЫХ С АЛГОРИТМОМ М* 27+и 51-.! Пг+ ЕЗ+ые еПрактически вти данные могут быть цривелеиы с намного большей точностью; в таблице лри- ведено только достаточное количество цифр, чтобы ннтересуюшиеги читатели могли более точно лровернть собственные алгоритмы вычислении значений. распределенная случайная величина, но работу можно проделать без применения каких-либо трансцендентных функций следующим интересным способом.
Присвоим УО 4 — й(Х) и будем генерировать обычные равномерно распределенные случайные величины 1"1„172, ..., пока не найдутея К > 1 е 1'к ! < Ук. Для фиксированных Х и !г вероятность того, что й(Х) > 1'! » 1е„равна 1/!41, умноженному на веРоатность того, что шах(1'1,..., Уа) < гг(Х), т. е. )!(Х)е//г!. Следовательно, вероятность того, что К = к, равна А(Х)" !/(74 — Ц! — !г(Х)е/!г(, а вероятность того, что К вЂ” нечетное число, равна гмх!'-' г(хГ) (, (й- ц! й! енечегиее,а>! (22) Следовательно, мы отбраковываем Х и начинаем снова, если К вЂ” четное число. Принимаем Х как случайную величину с плотностью распределения (20), если К— нечетное число. Обычно мы не генерируем много значений 1ху, чтобы определить К, поскольку среднее значение К (при заданном Х) равно 2 е>ОРг(К > !г) = ',Се>0 Й(Х)е/74! = еа(х! < е.
Несколькими годами позже Форсайт показал, что этот подход приводит к получению эффективного метода вычисления нормальных случайных величин без использования вспомогательных программ для вычисления квадратных корней нли логарифмов, как в алгоритмах Р и М. Его процедуру е улучшенным выбором интервалов (а..Ь), осуществленную И. Г.
Аренеом (Я. Н. А)!геле) и У. Дитером (П. 1!!еьег). можно представить в виде алго)!итма. Алгоритм Р (Метод "чет-нечет" для нормальных случайных величин), Этот алгоритм генерирует нормальные случайные величины на бинарном компьютере: предполагается, что имеетея приблизительно 1+ 1 двоичных разрядов точности.
Он нуждаетея в таблице значений ду = а — а ! для 1 < у < 1+ 1, где а определяется О .000 .067 ! .849 .161 2 .970 .236 З .855 .285 4 .994 .308 5 .995 .304 6 .933 .280 7 .923 .241 8 .727 .197 9 1.000 .152 10 .691 .112 1! .454,079 12 .287 .052 13 . ! 74 .033 14 .101 020 15 .057 .086 0.00 .236 — 0.92 .206 — 5.86 .234 — 0.58 .201 -33.13 .201 — 39.56 .214 — 2.57 ,217 — !.6! .275 0.67 .200 0.00 .289 0.35 .440 — О. 17 .698 0,92 1.150 0.36 1.974 — 0.02 3.526 О.!9 0.59 0.20 0,21 0.0 0,96 1.32 0.24 0.2 -0.06 6.66 0.26 0.4 0.12 !.38 0.28 0.6 1.31 34,93 0.29 0.8 0,3! 41,35 0.29 1.0 1.12 2.97 0.28 1,2 0.54 2,61 0,26 1.4 0.75 0.73 0.25 1.6 056 0,00 0.24 1.8 0.17 0.65 0.23 2,0 0.38 0.37 0.22 2.2 -0.0! 0.28 0.21 2.4 0.39 0.24 0,21 2.6 0.20 0.22 0.20 2.8 0.78 0.21 0.22 3.0 .505 2е.ОО .773 ! 2.50 .876 8.33 .939 6.25 .986 5.00 Л95 4.06 .987 3.37 .979 2.86 ,972 2.47 .966 2.16 .960 1,92 .954 1.71 .948 1Л4 .942 !.40 .936 1.27 отношением — е '~ Ых=23.
(23) Р1, [Получить Ц Генерировать равномерно распределенное число У=(.ЬоЬ1... Ь )ш где Ьо, Ь|, ..., Ь| означают двоичные разряды в бинарной записи. Присвоить Ф -Ьв,/ -1иа -О. Р2. [Найти первый нуль Ь .) Если Ь, = 1, присвоить в < — а+ду,,у +- 7'+ 1 и повторить этот шаг. (Если у = г+ 1, трактовать Ьз как нуль.) РЗ. [Генерировать Х.) (Сейчас а = а. ы и текущее значение у появится с вероятностью ~ 2 ~. Будем генерировать Х в интервале [ау 1..ау), используя метод отбраковки, описанный выше, с Ь(х) = хт/2 — аз/2 = рт/2+ ау, где р = х — а. В упр. 12 доказывается, что 5(х) < 1, как требуется в (21).) Присвоить У ~- д, умноженное на (.Ьчю ...Ьс)з и Р е- (11'+ а)У. (Так как среднее значение уб равно 2, обычно достаточно старших двоичных разрядов в (.Ьую ...Ь~)т для обеспечения приличной точности.
Вычисления без труда выполняются арифметикой с фиксированной точкой,) Р4. [Отбраковатьу[ Генерируем равномерно распределенную случайную величину С. Если Р < (/, перейти к шагу Р5. Иначе — присвоить Р новую равномерно распределенную случайную величину. Если сейчас У < Р (т. е. если К четное, как обсуждалась выше), возвратиться к шагу Р3, иначе — повторить шаг Р4. Рб. [Выход из программы по Х.[ Присвоить Х е- а + 1'. Если в' = 1, присвоить Х+- — Х.
! 4) Отношение равномерно распределенных случайных величав. Существует еще один хороший метод генерирования нормальных случайных величии„открытый в 1976 году А. Дж. Киндерманом (А. Л. К1пйегшап) и Дж. Ф. Монаханом (3. Р. Мопайап). Его суть — генерирование случайной точки (У, Р) в области, определенной как О < а < 1, -2а~/)п(1/а) < е < 2а~/Гп(1/а), (24) и получение отношения Х +- Р/(7. Заштрихованная площадь на рис. 13 — это магическая область неравенства (24), осуществляющая эту работу. Прежде чем изучать теорию, приведем алгоритм, эффективность и простота которого очевидны.
Значения ду для 1 < / < 47 появилясь в статье Аренса и Дитера (АЬгепв апб Вйегег, Маей. Сотр. 27 (1973), 927-937), В ней обсуждается усовершенствованный алгоритм, который повышает скорость его выполнения за счет использования больших таблиц. Алгоритм Р привлекателен тем, что работает почти с такой же скорость1о, как алгоритм М, и его легче выполнить. Среднее число равномерно распределенных случайных величин на каждую нормальную случайную величину равно 2.53947. Р. П. Брент (В.
Р. Вгепс, САСМ 17 (1974), 704-705) показал, как можно сократить это число до 1.37446 с помощью двух вычитаний и одного деления на каждую хранимую равномерно распределенную случайную величину. 10„/3/е) т/е я=1/3 я=о Рне. 13. Область "принятия" в методе отношения равномерных случайных величин для нормально распределенных случайных величин.