Рис. а. Площадь под графиком плотности распределения разделена на 31 часть. Площадь каждой части равна среднему числу вычислений случайной величины с такой плотностью.

Припишем случайный знак ее значению, и это сделает ее действительно нормальной случайной величиной.

Метод прямоугольника-клина-хвоста основан на важных общих технических приемах, которые будут рассмотрены ниже по мере построения алгоритма. Первая ключевая идея - рассматривать F{x) как смесь нескольких других функций, т. е. записать

(х), (13)

где Fi, F2, ..., F„ - подходящие распределения иpi, р2,...,Рп - неотрицательные вероятности, сумма которых равна 1. Если генерировать случайную переменную X, выбирая распределение Fj с вероятностью pj, то легко видеть, что А точно будет иметь F-распределение. С некоторыми распределениями Fj{x), пожалуй, трудно иметь дело, даже труднее, чем с F, но мы обычно устраиваем так, что вероятности Pj в этом случае очень малы. Большинство распределений Fj(x) будут довольно хорошо устроены, поскольку они будут простой модификацией равномерного распределения. Изложение завершается чрезвычайно эффективной программой, поскольку среднее время счета этой программы очень мало.

Рассматриваемый здесь метод легче понять, если работать с производными распределений, а не с самими распределениями. Пусть

f{x)=F{x), fj{x) = Fj{x)

будут плотностями распределений. Тогда равенство (13) можно записать как

f{x)=pifiix)+p2f2ix) + ---+Pnfn{x). (14)

Каждая fj{x) есть > О, и общая площадь под графиком fj{x) равна 1. Поэтому существует подходящий графический метод отображения зависимости (14): площадь под /(х) разделена на п частей и части, соответствующие fj{x), имеют площадь Pj. На рис. 9 иллюстрируется интересующий нас случай при /(х) = F{x) - у/2/пе~ /2. площадь под этой кривой разделена на п = 31 часть. Существует 15 прямоугольников, представляющихpi/i(x), ..., 15/15(3:), 15 клинообразных частей, представляющих Piefiei), Рзо1зо{х), и оставшаяся часть P3i/3i(x) - "хвост", т. е. график /(х) при х > 3.

Рис. 10. Плотность распределения, для которой алгоритм L может использоваться при генерировании случайных чисел.

Прямоугольные части fi{x), fisix) представляют равномерное распределение. Например, /з{х) представляет случайную равномерно распределенную величину, лежащую между и . Высота pj/j(ж) равна f{j/5); следовательно, площадь j-ro прямоугольника равна

257г

для 1 < j < 15.

(15)

Чтобы генерировать распределение, соответствующее таким прямоугольным частям, просто вычислим

X=\U + S, (16)

где и равномерно и 5 принимает значение (j - 1)/5 с вероятностью pj и не зависит от и. Так как pi + + рхь = .9183, можно просто использовать равномерно распределенные случайные величины, подобные этим приблизительно в 92% случаев.

В оставшихся 8% случаев будем обычно генерировать одно из клинообразных распределений Fie, F30. Типичный пример, который показывает, что необходимо делать, представлен на рис. 10. Когда х < \, часть кривой вогнутая, а когда X > 1, она выпуклая, но в каждом случае часть кривой достаточно близка к прямой линии и, как показано, может быть заключена между двумя параллельными линиями.

Чтобы перебрать эти клиновидные распределения, будем использовать другой общий технический прием: метод отбраковки фон Неймана получения сложной плотности из другой плотности, которая ее "заключает". Описанный выше метод полярньпс координат является простым примером такого подхода: шаги Р1-РЗ получают случайную точку внутри единичного круга, генерируя ее в большем круге, отбраковывая ее и начиная снова, если точка была вне круга.

Общий метод отбраковки является даже более сильным, чем этот. Пусть нужно генерировать случайную величину X с плотностью / и пусть д - другая плотность распределения, такая, что

/(<) < cg{t) (17)

для всех t, где с - константа. Тогда генерируем случайную величину X с плотностью д, а также независимую равномерно распределенную случайную величину U. Если и > f{X)/cg{X), отбрасываем X и начинаем снова с другими X и U. Когда

a/fc

Рис. 11. Область "принятия гипотезы" алгоритма L.

условие и < f{X)/cg{X) в конце концов выполняется, на выходе X будет иметь требуемую плотность /. [Доказательство. X < х произойдет с вероятностью Р{) = (5(t)A-/(t)/cp(t))+gp(x), где величина g = Л-(1-/(*)/ср(*))) =

1 - 1/с равна вероятности отбраковки; следовательно, р{х) - f{t) dt.

Техника отбраковки наиболее эффективна, когда с мало, так как должно быть в среднем с итераций, прежде чем значение будет принято (см. упр. 6). В одних случаях f{x)/cg{x) всегда равно О или 1 и нет необходимости генерировать U. В других случаях, если f{x)/cg{x) трудно вычислить, следует постараться "втиснуть" его между двумя более простыми граничными функциями

г{х) < f{x)/cg{x) < s{x)

(18)

и точное значение f{x)/cg{x) не нужно вычислять, если не выполняется неравенство г{х) < и < s{x). Следзющий алгоритм разрешает проблему клина, совершенств}я метод отбраковки.

Алгоритм L (Плотности, близкие к линейным). Этот алгоритм можно использовать для генерирования случайной величины Л" с любым распределением, плотность f(x) которого удовлетворяет следующим условиям (см. рис. 10);

f(x) = О для а; < S и для х > s-\- h;

a-b(x - s)/h < f(x) <b-b(x - s)/h для s < x < s + h.

(19)

LI. [Получить и < v.] Генерировать две независимые случайные величины U и V, разномерно распределенные между О и 1. Если U > V, заменить U -V.

L2. [Простой случай?] Если V < а/Ь, перейти к шагу L4.

L3. [Попытаемся снова?] Если V > U + (l/b)f(s + hU), возвратиться к шагу L1. (Если а/Ъ близко к 1, данный шаг алгоритма нужен не очень часто).

L4. [Вычисление X.] Присвоить X s + hU.

Когда достигнут шаг L4, точка (U, V) является случайной точкой на площади, заштрихованной на рис. И, а именно - О < U < V < U+(l/b)f{s+hU). Условия (19) гарантируют, что

Т <и +lf(s+hU) < 1.