Д. Кнут - Искусство программирования том 2 (3-е издание) - 2001 (Часть 1) (1119452), страница 39
Текст из файла (страница 39)
[Вычисление Я,) Присвоить Я ! У»» + У»», РЗ. !Проверить Я > 17) Если 5 > 1, возврат к шагу Р1. (Шаги Р1-РЗ выполняются в среднем 1.27 раз со среднеквадратичным отклонением, равным 0.587; см. упр. 6.) Р4, !Вычисление Хм Х».) Присвоить Х» и Л» следующие значения: -21пЯ -21п5 Х1+ У»Ч Х»+ У»Ч я 1/ я (11) Это требуемые нормапьно распределенные случайные величины. $ Х, = 4:2!пЯг!п8. Х» = ~/-2!плсов8, Используя также полярные координаты Х» = В'соз8' и Х» = В'г!п8', получим 8' = 8 и В' = ~/-2 !и л. Ясно, что В' и 8' независимы, поскольку В и 8 независимы в единичном круге.
К тому же 8' равномерно распределено между 0 и 2я, и вероятность того, что В' < г, равна вероятности, что -2 1п Я < г», т. е. вероятности, что Я > е '7». Эта вероятность равна 1- е '~», так как 5 = В» равномерно распределено мвклу 0 и 1. Вероятность того, что В' лежит между г и г+Й', поэтому с»» '',»» равна дифференциалу от 1 — е ' 7», т.
е. ге ' 7» Й, Аналогично вероятность того, что 8' лежит между д и й+ ар, равна !1(2я) Ю Совместная вероятность того, что Х» < х» и Х» < х», равняется — е '~ гЙФ -'» !!кг)! сьвгйьь иьгьг») 2я е !'+" )7 Ихпу 2к /!(,.„)1,,(„и г(„) — е *~»йх — е г )~Ну Это и доказывает, что Х» и Л» независимы и нормально распределены. 2) Мепьк) прямоугольника-клина-хаос»иа предложен Дж.
Марсалья. Здесь используется функция Р(х) =ег!!х/й2) = ~ — )) е '7»й, х>0, о !12) которая является функцией распределения абсолюя»ного гиачгиил нормальной случайной величины. Затем Х вычисляется в соответствии с распределением (12). Для доказательства законности данного метода используем элементарную аналитическую геометрию и вычисления: если на шаге РЗ Я < 1, точка плоскости с декартовыми координатами Я, У») является случайной точкой, равномерно распределенной внутри единичного круга. Перейдя к полярным координатам Ъ~ = В сог 8, У» -- Вз!п 8, получим о.э 0.8 О.е о.з Ол ол од ол о.о о Рнс.
9. Площадь пад графиком плотности распределения разделенана31 часть. Площадь каждой части равна среднему числу вычислений случайной величины с такой плотностью. Припишем случайный знак ее значению, и это сделает ее действительно нормальной случайной величиной. Метод прямоугольника-клина-хвоста основан на важных общих технических приемах, которые будут рассмотрены ниже по мере построения алгоритма. Первая ключевая идея — рассматривать Е(х) кэк смесь нескольких других функций, т. е.
записать Г(х) = р1 Г~(х) + рзРз(х) + ° ° ° + р„г'„(х), (13) где г м гз ° °, гя — подходящие распределения и рз, рз, ..., р„— неотрицательные вероятности, сумма которых равна 1. Если генерировать случайную переменную Л, выбирая распределение Ру с вероятностью рЗ, то легко видеть, что Х точно будет иметь Г-распределение. С некоторыми распределениями Рз(х), пожалуй, трудно иметь дело, даже труднее, чем с Р, но мы обычно устраиваем так, что вероятности рз в этом случае очень малы.
Большинство распределений Р (х) будут довольно хороню устроены, поскольку они будут простой модификацией равномерного распределения. Изложение завершается чрезвычайно эффективной программой, поскольку среднее время счета этой программы очень малб. Рассматриваемый здесь метод легче понять, если работать с производнмми распределений, а не с самими распределениями. Пусть Дх) = Р'(х), Д(х) = г)'(х) будут плопзносщлми распределений. Тогда равенство (13) можно записать как 1(х) = рзл(х) + ртах) + - ° +р„У„(х). (14) Каждая Д(х) есть > О, и общая площадь под графиком Ях) равна 1.
Поэтому существует подходящий графический метод отображения зависимости (14): плошадь под Дх) разделена на и частей и части, соответствующие Д(х), имеют площадь р,. На рис. 9 иллюстрируется интересуюший нас случай при з(х) = Г'(х) ~/2/ге * /з:, площадь под этой кривой разделена на и = 31 часть. Существует Уб пРЯмоУгольников, пРедставлающих Рз уз (х),..., Р| 8~~ 8 (х), 13 клинообРазных частей, представляющих р18Лз(х)., рзоХза(х), и оставшаяся часть рз1 Ьз(х) — "хвост", т. е. график Дх) при х > 3.
о Рис. 10. Плотность распределения, лля которой алгоритм 1. может использоваться прн генерированны случайнык чисел. Прямоугольные части /~(з), ..., /гв(я) представляют равномерное распределение. Например, /э(к) представляет случайную равномерно распределенную величину, лежащую между ь н з.
Высота р /1(л) равна /(у/5); следовательно, площадь у-го прямоугольника равна рз — — -/(у/5) = ~( —. е 11ЬЕ для 1 (/ < 15. -Р ьо 5 з' 25я (15) Чтобы генерировать распределение, соответствующее таким прямоугольным частям, просто вычислим Х= ~и+Ю, (16) где У равномерно и Я принимает значение (/ — 1)/5 с вероятностью рт и не зависит от У. Так как р1 + "° + р~ь = .9183, можно просто использовать равномерно распределенные случайные величины, подобные этим приблизительно в 92% случаев. В оставшихся 8% случаев будем обычно генерировать одно из клинообразных распределений Гш, ..., Гзе.
Типичный пример, который показывает. что необходимо делать, представлен на рнс. 10. Когда я к 1, часть кривой вогнутая, а когда я > 1, она выпуклая, но в каждом случае часть кривой достаточно близка к прямой линии н, Как показано, может быть заключена между двумя параллельнымн линиями. Чтобы перебрать эти клиновидные распределения, будем использовать другой общий технический прием: мешод ошбракоекп фон Неймана получения сложной плотности из другой плотности, которая ее "заключает". Описанный выше метод полярных координат является простым примером такого подхода: шаги Р1.-РЗ получают случайную точку внутри единичного круга, генерируя ее в большем круге, отбраковывая ее и начиная снова, если точка была вне круга. Общий метод отбраковки является даже более сильным, чем этот.
Пусть нужно генерировать случайную величину Х с плотностью / и пусть д — другая плотность распределения, хакан, что /(1) < сд(Г) (17) для всех 1, где с — константа. Тогда генерируем случайную величину Х с плотностью д, а также независимую равномерно распределенную случайную величину Г Если У > /(Х)/сд(Х), отбрасываем Х и начинаем снова с другими Х н Г Когда Рис.
11. Область "принятия гипотезы" алгоритма 1.. условие (? < /(Х)/су(Х) в конце концов выполняется, на выходе Х будет иметь требуемую плотность /. [Докаоательсшоо. Х < х произойдет с вероятностью р(х) = ) (у(Ь) аь /(Ф)/су(Ф))+ур(х), где величина д = ) (д(1) аь (1-/(г)/сд(1))) = 1 — 1/с равна вероятности отбраковки; следовательно, р(х) = ) /($) аь,) Техника отбраковки наиболее зффективна, когда с малб, так как должно быть в среднем с итераций., прежде чем значение будет принято (см. упр. 6). В одних случаях /(х)/су(х) всегда равно О или 1 и нет необходимости генерировать (/. В других случаях, если /(х)/сд(х) трудно вычислить, следует постараться "втиснуть" его между двумя более простыми граничными функциями г(х) < /(х)/су(х) < о(х) и точное значение /(х)/сд(х) не нужно вычислять, если не выполняется неравенство г(х) < П < о(х).
Следующий алгоритм разрешает проблему клина, совершенствуя метод отбраковки. Алгоритм Е (Плошиости, близкие к линейным). Этот алгоритм можно использовать для генерирования случайной величины Х с любым распределением, плотность /(х) которого удовлетворяет следующим условиям (см, рис. 10): /(х) = 0 для х < о и для х > о+ Ь; (19) а — Ь(х — а)/й < /(х) < Ь - Ь(х — о)/Ь для э < х < г + 6. Е1. (Получить 1? < К) Генерировать две независимые случайные величины 1? и !', равномерно распределенные между 0 и 1, Если (? > г', заменить (/ т+ 1'. 1 2. (Простой случайу) Если 1' < а/Ь, перейти к шагу 1.4. 13. [Попьпаемся снова?] Если 1' > (/+ (1/Ь)/(о+ ЬП), возвратиться к шагу 1.1. (Если а/Ь близко к 1, данный шаг алгоритма нужен не очень часто).
Е4. (Вычисление Х.) Присвоить Х ~- з+ ЬГ. $ еогда достигнут шаг 14, точка (1/, 1') является случайной точкой на плошади, заштрихованной на рис. 11, а именно — 0 < П < г' < (?+(1/Ь)/(о+ЬП). Условия (19) гарантируют, что а 1 — < (/+ -/(о+ Ь(?) < 1. Ь Ь Рнс.