1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 25
Текст из файла (страница 25)
подразделы 6.4 и 13.1); для дискретных случайных величин это метод Уолкера (алгоритм 10.6);– специальные методы, основанные на свойствах конкретных распределений; для непрерывных случайных величин (как примеры) это методы моделирования гамма- бета- и гауссовского распределений (см. разделы 12, 13); для дискретных случайных величин этоалгоритмы 10.5 и 10.8.Напомним, что непрерывные случайные величины принимают континуум значений из множеств типа интервал, прямоугольный параллелепипед, шар и т. п.В этом разделе мы рассмотрим более «простые» (с точки зрениятеории вероятностей, см., например, [14]) дискретные случайные величины, имеющие конечные или счетные множества значений. Необ152ходимость численного моделирования выборочных значений таких случайных величин возникает, в частности, при моделировании случайныхсобытий (см., например, алгоритм 2.7) и при реализации упомянутого здесь метода дискретной суперпозиции (см., в частности, алгоритмы 4.2 и 11.1).Итак, рассмотрим дискретную случайную величину ξ, принимающую значения x1 , .
. . , xM и имеющую распределениеP{ξ = xi } = pi ; pi > 0,MXpi = 1; M ≤ ∞.(10.1)i=1В книгах [5, 9, 23] отмечено, что вместо «табличного» представления(10.1) распределения дискретной случайной величины ξ можно рассматривать дельта-плотность видаfξ (u) =MXpi δ(u − xi ),(10.2)i=1Rгде δ(u) – дельта-функция Дирака: y(u)δ(u − xi ) du = y(xi ) для любойнепрерывной в точке xi функции y(u).Более того, в этих же книгах отмечено, что представленный в следующем подпункте алгоритм моделирования дискретной случайной величины можно трактовать как специальный случай метода обратнойфункции распределения (см. алгоритм 2.4), если вместо функции Fξ−1 (x)рассматривать функцию G(x) = inf y:x<Fξ (y) y (ведь функция распределения случайной величины ξ с распределением (10.1) – кусочнопостоянна).Однако такая трактовка годится только для частного случаяx1 < x2 < . .
. < xM ,(10.3)а это условие не всегда обеспечивает экономичность алгоритма моделирования (см. далее утверждение 10.1).Таким образом, для изучения вопросов моделирования дискретнойслучайной величины ξ удобнее использовать «табличное» представление распределения (10.1). Дельта-плотность (10.2) пригодится нам вразделе 11 для представления метода дискретной суперпозиции.15310.2. Стандартный алгоритм моделирования дискретнойслучайной величины для случая M < ∞ и его трудоемкость.Рассмотрим сначала случай, когда число значений M случайной величины ξ с распределением (10.1) является конечным и «умеренно малым»; это означает, что числа {x1 , .
. . , xM } и {p1 , . . . , pM } заданы явно.Случай M 1 (вплоть до M = ∞) рассмотрен далее в подразделе 10.4.В силу соотношений (10.1), численное моделирование того или иноговыборочного значения ξ0 = xm случайной величины ξ означает розыгрыш события с вероятностью pm . Кроме того, из соотношений (10.1)следует, что интервал (0, 1) можно разбить на полуинтервалы ∆m длины pm :∆m = [Rm−1 , Rm );Rm =mXpi ; m = 1, 2, . .
. , M ;i=1для m = 1 полагаем Rm−1 = R0 = 0.Используя соответствующий генератор, численно моделируем выборочное значение α0 ∈ U (0, 1) стандартного случайного числа. В силусоотношения (9.6) (или (2.14)) для α ∈ U (0, 1) имеем P{α ∈ ∆m } = pm .Таким образом, если α0 ∈ ∆m , то для данного испытания следуетположить, что выборочное значение ξ0 случайной величины ξ равноξ0 = xm .Технически определение того номера m полуинтервала ∆m , в который попало выборочное значение α0 , осуществляется последовательнымвычитанием из α0 сумм Rm для m = 1, 2, . .
. до тех пор, пока разностьα0 − Rm не станет отрицательной (рис. 10.1).Рис. 10.1. Стандартный алгоритм моделирования дискретной случайнойвеличиныОписанную операцию можно осуществлять без непосредственноговычисления сумм Rm , используя операцию переприсваивания.154АЛГОРИТМ 10.1 (см., например, [9, 13]). Моделируем выборочноезначение Q := α0 ∈ U (0, 1) (т. е. Q := RAN DOM ) и полагаем m := 1.Производим переприсваиваниеQ := Q − pm(10.4)(т. е. заносим новое значение (α0 − p1 ) в ячейку Q). Если новое Qне положительно, то в качестве m выбираем текущее его значениеи полагаем ξ0 = xm , в противном случае делаем переприсваиванияm := m + 1 и (10.4), производим проверку Q на положительность ит.
д.Важным примером дискретных случайных величин являются целочисленные случайные величины η, для которыхxi = i, P{η = i} = pi ; i = 1, 2, . . . , M.(10.5)Алгоритм 10.1 в этом случае имеет следующий вид.АЛГОРИТМ 10.2 (см., например, [9, 13]). Моделируем значениеQ := α0 ∈ U (0, 1) и полагаем m := 1. Производим переприсваивание(10.4). Если новое значение Q не положительно, то полагаем η0 = m,в противном случае делаем переприсваивания m := m + 1 и (10.4), производим проверку Q на положительность и т.
д.В качестве целочисленных можно также рассматривать случайныевеличины, принимающие всевозможные (не только натуральные,а еще 0 и отрицательные) значения из множества целых чисел Z; см.,в частности, примеры и рассуждения подраздела 10.3.
Однако чащевсего значения целочисленной случайной величины имеют смысл натурального числа (например, числа ветвлений траекторий – см. подразделы 3.4 и 10.3) или номера (как, например, в методе дискретной суперпозиции), т. е. используется распределение (10.5) и соответствующий алгоритм 10.2.Оценим трудоемкость (затраты) стандартного алгоритма 10.1. Прежде всего сформулируем следующее замечание.ЗАМЕЧАНИЕ 10.1. В подразделе 1.9 указано, что в формулахs = t × n и S = t × Dζ (см. также соотношения (1.16), (1.18)) для затрат итрудоемкости метода Монте-Карло (1.1) величина t является средним временем вычисления одного выборочного значения ζj .
Слово среднее здесьпринципиально, так как затраты на моделирование конкретных выборочных значений ζj чаще всего являются случайными.155Это проявляется уже на простейшем примере случайной величины ξ сдискретным распределением (10.1) (если ставится задача вычисления маPMтематического ожидания Eξ = i=1 xi pi методом Монте-Карло, например, с целью тестирования генераторов стандартных случайных чисел).Легко видеть, что в случае, когда ξ0 = xm , приходится осуществлятьm проверок Q на положительность.
Трудоемкость δ алгоритма 10.1 включает затраты на моделирование одного стандартного случайного числа α0 ∈ U (0, 1) (эти затраты обозначим a) и на реализациюсравнений Q с нулем (затраты на каждое сравнение обозначим b).Тогда средняя трудоемкость t алгоритма 10.1 равна!MXt = Eδ = a +m pm × b.(10.6)m=1Заметим, что для целочисленной случайной величины η с распределением (10.5) величина t из (10.6) представима в видеt1 = a + b Eη.(10.7)При реализации алгоритмов 10.1 и 10.2 целесообразно (если это возможно) располагать вероятности pi в порядке их убывания.УТВЕРЖДЕНИЕ 10.1 [9, 23]. Оптимальным порядком нумерациизначений {x1 , .
. . , xM } и вероятностей {p1 , . . . , pM } из (10.1), при котором средние затраты t из (10.6) минимальны, является такой, прикоторомp1 ≥ p2 ≥ . . . ≥ pM .(10.8)ДОКАЗАТЕЛЬСТВО. Рассуждаем от противного. Пусть для неко00торого распределения вероятностей P{p1 , . . . , pM } соотношение (10.8) неM00выполнено, а величина t = a +m=1 mpm × b минимальна.
Тогданайдутся такие номера s и k, для которых одновременно выполненоs < k и p0s < p0k .(10.9)Рассмотрим новое распределение вероятностей {p001 , . . . , p00M }, котороеполучено из распределения {p01 , . . . , p0M } перестановкой вероятностей p0sи p0k , т. е. p00j = p0j при j 6= s, j 6= k и p00s = p0k , p00k = p0s .Рассмотрим разностьt0 −t00 = (sp0s +kp0k −sp00s −kp00k )b = (sp0s +kp0k −sp0k −kp0s )b = (k−s)(p0k −p0s )b.156Из соотношения (10.9) следует, что t0 − t00 > 0 и t00 < t0 .
Получилипротиворечие с тем, что величина t0 минимальна. Утверждение 10.1доказано.Таким образом, использование аналога метода обратной функциираспределения для случая (10.3) (см. подраздел 10.1) может дать неоптимальный вариант алгоритма 10.1, так как в этом случае могут бытьне выполнены неравенства (10.8).10.3. Случай малого числа значений. Теорема оминимальной дисперсии целочисленной случайной величины.Отметим, что в алгоритмах 10.1, 10.2 возможна следующая модификация.Если α0 −RM −1 > 0, то последнее вычитание можно не производить,так как в этом случае α0 ∈ ∆M . Соответствующий аналог формулы(10.6) имеет вид!M−1Xt=a+m pm + (M − 1) pM × b.m=1Однако выигрыш от этой модификации сказывается только для малых M , так как при ее применении в алгоритмах 10.1, 10.2 требуетсяпроводить сравнение текущего m c (M − 1).Рассмотрим, в частности, случай M = 2.
Здесь модифицированныйалгоритм 10.1 приобретает следующий простой вид.АЛГОРИТМ 10.3 (см., например, [9, 13]). Моделируем значениеα0 ∈ U (0, 1). Если α0 < p1 , то ξ0 = x1 , иначе ξ0 = x2 .Справедлива также формула С. А. Роженко [23]:ξ0 = x[2−p1 +α0 ] ,(10.10)где [A] обозначает целую часть числа A. Действительно, при α0 < p1число 2 − p1 + α0 принадлежит полуинтервалу [1, 2) и [2 − p1 + α0 ] = 1;иначе 2 − p1 + α0 ∈ [2, 3) и [2 − p1 + α0 ] = 2.В случае, когда x1 = 1 и x2 = 0 величина ξ является бернуллиевской случайной величиной с вероятностью «успеха» p1 (см., например, [14]).Здесь уместно напомнить, что бернуллиевская случайная величинавводится для математической формализации изучения случайных событий.
Если придать x1 и x2 «словесные» значения типаx1 = {событие A случилось} и x2 = {событие A не случилось}, то157алгоритм 10.3 определяет численное моделирование события A.Частным случаем такого метода является алгоритм 2.7, где(j−1) A = {произошел обрыв траектории} и p1 = p(a) ξ0.Наконец, важный случай появления дискретных случайных величин с двумя целыми значениями (типа бернуллиевских) возникает вследующей задаче оптимизации алгоритмов, связанных с расщеплением (ветвлением траекторий) – см.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
