должно произойти со следующими 16 вероятностями.

j = О 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 - J J 36 36 36 36 36 36 36 36 36 36 36

Это можно осуществить, используя (3), если А; = 16 и Xj+i = j при О < j < 16 и если таблицы для Р и Y имеют следующий вид.

j = О 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 i.-O О i I 1 I 1 1 1 i i I I О О О Fj = 5974*6***847 10 678

(Когда Pj = 1, Yj не используются.) Например, значение 7 встречается с вероятностью ] • ((1 - Р2) -I- Рт + (1 - Рц) + (1 - Р14)) = ! как и требуется. Это необычный способ бросания игральных костей, но результаты получаются такие же, как и в реальной ситуации.

Вероятности pj, безусловно, могут быть представлены неотрицательными весами wi, W2, •.., lUk, если обозначить сумму весов через W, то pj = Wj/W. В разных применениях отдельные веса весьма изменчивы. Матиас, Виттер и Ни (см. работу Matias, Vitter, and Ni, SODA 4 (1993), 361-370) показали, как изменять веса и генерировать X с постоянным средним временем.

В. Общие методы для непрерывных распределений. В общем случае распределение действительных чисел может быть выражено в терминах "функции распределения" F{x), которая точно определяет вероятность того, что случайная величина X не превысит значение х:

F{x) = Рг(Х < х). (4)

Эта функция всегда монотонно возрастает от О до 1, т. е.

F{xi)<F{x2), еслих1<Х2; F(-oo) = О, F(-i-c») = 1. (5)

Примеры функций распределения приведены в разделе 3.3.1 (см. рис. 3). Если F{x) непрерывна и строго возрастающая (так что F{xi) < F(x2), когда xi < Х2), то она принимает все значения между О и 1 и существует обратная функция F[~il(y), такая, что для О <у < 1

у = F{x) тогда и только тогда, когда х = F~\y). (6)

В большинстве случаев, когда F(x) непрерывна и строго возрастающая, можно вычислить случайную величину X с распределением F{x), полагая

Х = р[-11(С7), (7)

где и - равномерно распределенная случайная величина. Действительно, вероятность того, что X < X, равна вероятности, что Ff~4(f/) < х, а именно - вероятности того, что и < F{x), т. е. F{x).

Теперь проблема сводится к решению задачи численного анализа - к нахождению хороших методов вычисления Ff-il([7) с требуемой точностью. Численный

анализ в этой книге о получисленных алгоритмах не рассматривается, одаако существует ряд важных методов, способных улучшить общий подход (7), и здесь они будут рассмотрены.

Заметим, что если Xi - случайная величина, имеющая функцию распределения Fi(x), и если Х2 - независимая от Ху случайная величина с функцией распределения F2{x), то

max(Xi,X2) имеет распределение Fi{x)F2{x),

min(Xi,J\r2) имеет распределение Fi(x) -f 2(х) - Fi(x)F2(x).

(См. упр. 4.) Например, равномерно распределенная случайная величина U имеет распределение F{x) - х для О < х < 1; если Ui, U2, , Ut - независимые равномерно распределенные случайные величины, то max(l[/i, f/2, ,Ut) имеет функцию распределения F{x) = х* при О < х < 1. Эта формула является основой критерия "максимум-", описанного в разделе 3.3.2. Обратная функция равна F~{y) = Vy-В частном случае при t = 2 получаем, следовательно, что формулы

X = VU и X = max{Ui,U2) (9)

дадут одинаковое распределение случайной величины X, хотя, на первый взгляд, это не очевидно. Нет необходимости извлекать квадратный корень из равномерно распределенной случайной величины.

Количество подобных хитростей бесконечно: любой алгоритм, использующий случайные числа на входе, дает на выходе случайные величины с некоторым распределением. Задача состоит в нахождении общих методов составления алгоритма, обеспечивающего заданную функцию распределения на выходе. Вместо того чтобы рассматривать подобные методы в исключительно абстрактных терминах, изучим, как они могут применяться в важных случаях.

С. Нормсьпьное распределение. Возможно, наиболее значительным неравномерным, непрерывным распределением является нормальное распределение с нулевым средним значением и среднеквадратичным отклонением, равным единице:

F{x) = Г e-l4t. (10)

Значительность данного распределения показана в разделе 1.2.10. В нашем случае обратную функцию F-I не так легко вычислить; но, как мы увидим, существует несколько технических приемов моделирования этого распределения.

1) Метод полярных координат, предложенный Дж. Э. П. Боксом, М. Э. Мюллером и Дж. Марсалья [см. G. Е..Р. Box, М. Е. Muller, and G. Marsaglia, Annals Math. Stat. 29 (1958), 610-611, и Boeing Scientific Res. Lab. report Dl-82-0203 (1962)].

Алгоритм P (Метод полярных координат для нормальных случайных величин). Этот алгоритм вычисляет две независимые нормально распределенные случайные величины: Xi и Х2.

Р1. [Получение равномерно распределенных случайных величин.] Генерируем две независимые случайные величины Ui и С/2, равномерно распределенные между О

и 1. Присвоить Vi 2Ui - 1, V2 -t- 2U2 - 1- (Здесь Vi и V2 равномерно распределены между -ё1 и -1-1. На большинстве компьютеров предпочтительнее представление T/i и V2 в виде чисел с плавающей точкой.)

Р2. [Вычисление S.] Присвоить S -i-V + V.

РЗ. [Проверить 5 > 1?] Если 5 > 1, возврат к шагу Р1. (Шаги Р1-РЗ выполняются в среднем 1.27 раз со среднеквадратичным отклонением, равным 0.587; см. упр. 6.)

Р4. [Вычисление Xi,X2-] Присвоить Xi и Х2 следующие значения:

/-21nS /-2lnS

XiVi-, X2V2-. (11)

Это требуемые нормально распределенные случайные величины.

Для доказательства законности данного метода используем элементарную аналитическую геометрию и вычисления: если на шаге РЗ S < 1, точка плоскости с декартовыми координатами (Vi, Т/г) является случайной точкой, равномерно распределенной внутри единичного крут. Перейдя к полярным координатам Vi = Д cos 0, V2 = i?sin 0, получим

S = R\ Xi =\/-21n5cos0, X2 = \/-21n5sin0.

Используя также полярные координаты Xi = RcosQ и Х2 = i?sin0, получим 0 = 0 и i? = \/-2lnS. Ясно, что R и 0 независимы, поскольку RhQ независимы в единичном круге. К тому же 0 равномерно распределено между О и 27г, и вероятность того, что R < г, равна вероятности, что -2lnS < г, т. е. вероятности, что S > е~/2. Эта вероятность равна 1 - е"/ так как S = R равномерно распределено между О и 1. Вероятность того, что R лежит между г и r+dr, поэтому равна дифференциалу от 1 - е~ Z, т. е. ге"*" dr. Аналогично вероятность того, что 0 лежит между 9 и 9 + d9, равна (1/27г) d9. Совместная вероятность того, что Ai < xi и 2 < Х2, равняется

l-e-/\drd9

{(r,fl)rcosfl<xi, rsine<x2}

27Г /(,

Это и доказывает, что Xi и Х2 независимы и нормально распределены.

2) Метод прямоуголъника-клина-хвоста предложен Дж. Марсалья. Здесь используется функция

F{x) = erf(x/V) = \[\ e-dt, х > О,

(12)

которая является функцией распределения абсолютного значения нормальной случайной величины. Затем X вычисляется в соответствии с распределением (12).