1626435387-5d0ec7dd22f55dfcc65f3209f497fb0b (844202), страница 35
Текст из файла (страница 35)
Этивеличины относительно невелики (близки к единице).Например, для µ ≈ ν ≈ 2, 5 имеемs(2) ≈ s(3) ≈=Γ(5)Γ(2, 5)Γ(2)×=Γ(4, 5)Γ(2, 5)Γ(2, 5)15364!=≈ 1, 55;(7/2)(5/2)(3/2)2 (1/2)2 π315π√здесь использованы соотношения (12.6), (12.10), (12.11) и Γ(1/2) = π.(1)Однако следует учитывать, что при реализации величины ξ0 впунктах 1 алгоритмов 12.6, 12.7 требуется [ν] (или [µ]) обращений кдатчику случайных чисел (см. замечание 9.2).Отметим также, что для нецелых значений µ и ν вычисление зна(1) чений g ξ0в алгоритмах 12.5–12.7 может оказаться трудоемким. Вэтом случае можно использовать двусторонний метод исключения (алгоритм 11.13) с кусочно-постоянными мажорантой и минорантой.21013. Моделирование гауссовcких случайныхвеличин, векторов, процессов и полей13.1. Использование изотропного вектора случайнойдлины.
Рассмотрим вопрос о моделировании выборочных значений2(m,σ 2 )ξiслучайной величины ξ (m,σ ) , имеющей плотность гауссовского (нормального) распределения1(m,σ 2 )−(u−m)22σ 2√fξe(u) =, −∞ < u < +∞(13.1)σ 2π22с параметрами (m, σ 2 ) : m = Eξ (m,σ ) , σ 2 = Dξ (m,σ ) .(m,σ 2 )Заметим, что для получения значений ξi, распределенных со(0,1)гласно плотности (13.1), достаточно моделировать значения ξiстандартной гауссовской (нормальной) случайной величины ξ (0,1) ,имеющей плотность распределения(0,1)fξ21e−u /2 , −∞ < u < +∞,(u) = √2π(m,σ 2 )(0,1)(13.2)а затем использовать преобразование ξi= m + σξi .Заметим также, что, в отличие от бета- и гамма-распределений (12.1),(12.2), ни при каких значениях параметров m и σ 2 (в том числе дляm = 0 и σ 2 = 1) плотность (13.1) не является элементарной (см. рассуждения подраздела 2.6), и приходится конструировать специальные(отличные от стандартного метода обратной функции распределения– см.
алгоритм 2.4) алгоритмы моделирования выборочных значений(0,1)ξi .Для построения таких алгоритмов используется следующий факт.УТВЕРЖДЕНИЕ 13.1 (см., например, [9]). Если ξ – d-мерный изотропный случайный вектор (см. определение 6.1), квадрат длины которого имеет χ2 -распределение с d степенями свободы (см. формулу(12.8)), то его компоненты ξ (1) , . . . , ξ (d) независимы и нормальны с параметрами (0, 1) (т. е. распределены согласно плотности (16.2)).ДОКАЗАТЕЛЬСТВО.Рассмотримd-мерную «сферическую» систему координат r, w(1) , ..., w(d−1) , в которой компонента r описываетрасстояние от описываемой точки до начала декартовых координат,а компоненты {w(1) , ..., w(d−1) } описывают точку d-мерной единичнойсферы.
В качестве примера такой системы координат можно рассмотреть, в частности, гиперсферические координаты (6.19).211В рассматриваемой системе координат случайный вектор ξ имееткомпоненты χ̂(d) , φ̃(1) , . . . , φ̃(d−1) , где φ̃(1) , . . . , φ̃(d−1) – случайная точка, равномернораспределенная на d-мерной единичной сфере,pа χ̂(d) = χ2(d) . Здесь χ2(d) – случайная величина, имеющая χ2 -распределение с d степенями свободы (см.
формулу (12.8)).Функция распределения случайной величины χ̂(d) имеет вид:Fχ̂(d) (r) = P χ̂(d) < r = P χ2(d) < r2 =Zr20=1d/22 Γ(d/2)Ztd/2−1 e−t/2dt =2d/2 Γ(d/2)r2td/2−1 e−t/2 dt,r ≥ 0.0Дифференцируя полученную функцию по r, получаем плотность2112 d/2−1 −r 2 /2rd−1 e−r /2 .(2r) = d/2−1re2d/2 Γ(d/2)2Γ(d/2)В силу изотропности вектора χ̂(d) , φ̃(1) , . .
. , φ̃(d−1) и очевидной независимости компоненты χ̂(d) и вектора φ̃(1) , . . . , φ̃(d−1) , плотность распределения вектора ξ в рассматриваемой d-мерной «сферической» системе координат имеет видfχ̂(d) (r) =f(1)ξr, w(1) , .., w(d−1) =12d/2−1 Γ(d/2)rd−1 e−r2/2×1Ŝd=21rd−1 e−r /2 ,d/2(2 π)d/2где Ŝd = 2πΓ(d/2) – «площадь» поверхности d-мерной единичнойсферы (на самом деле это объем размерности (d −1)).Якобиан перехода от сферических координат r, w(1) , ..., w(d−1) кдекартовым координатам x = x(1) , .
. . , x(d) равен ∂ r(x), w(1) (x), .., w(d−1) (x) 11=h, = d−1(1)(d)2 i(d−1)/2r(x)2∂ x , .., xx(1) + .. + x(d)так как из соображений подобия бесконечно малый «декартовый» объемdV (x) выражается в «сферических» координатах следующим образом:dV (x) = rd−1 (x) dr(x)dw(1) (x) . . . dw(d−1) (x).212Используя утверждение 6.2, получаем плотность распределения вектора ξ в декартовых координатах: ∂ r(x), w(1) (x), .., w(d−1) (x) (1)fξ (x) = fr(x), w(1) (x), . .
. , w(d−1) (x)) =(1)(d)ξ∂ x , .., xh22 i h i(d−1)/2 − x(1) +..+ x(d)21(1) 2(d) 2=x+ .. + xe×d/2(2 π) d Y1 − x(i) 2 21√ e=.×h22 i(d−1)/22πi=1x(1) + .. + x(d)Из последнего соотношения следует, что ξ (1) , . . . , ξ (d) независимы иодинаково распределены согласно плотности (13.2). Утверждение 13.1доказано.(0,1)Таким образом, для получения d независимых значений ξiможноиспользовать следующий алгоритм.АЛГОРИТМ 13.1 (см., например, [9, 23]). 1. Моделируем выборочное значение ρ̃0 случайной величины ρ̃, имеющей гамма-распределениес параметрами λ = 1/2 и ν = d/2, а также выборочное значение(1)(d) ω0 , .
. . , ω0единичного d-мерного изотропного случайного вектора(см. определение 6.1).√(i)(0,1)2. Полагаем ξi= ρ̃0 × ω0 ; i = 1, . . . , d.Алгоритм 13.1 имеет ограниченное применение по следующим причинам.Во-первых,алгоритмнеэффективендлянечетныхd = 2k − 1, k = 1, 2, ... по той причине, что затруднено моделированиегамма-распределения с нецелым параметром (см. подразделы 12.3, 12.4).Во-вторых, для больших четных d = 2k (начиная с k = 2) имеютсятрудности с моделированием единичного изотропного вектора.Как показано в подразделах 6.3, 6.4, здесь, согласно утверждению 6.1, построение соответствующего алгоритма сводится к решениюпроблемы 6.1, связанной с моделированием точки, равномерно распределенной в d-мерном шаре B (d,0,R) фиксированного радиуса R (см. соотношение (6.7)).
В этих же подразделах 6.3 и 6.4 показано, что проблема 6.1 эффективно разрешима только для размерностей d = 1, 2, 3 (с помощью перехода к полярным координатам – для d = 2 и к сферическимкоординатам – для d = 3).213Поэтому алгоритм 13.1 по сути применяется только для размерности d = 2.Здесь распределение случайной величины ρ̃ является экспоненциальным (это гамма-распределение с параметром ν = 1 – см.
формулы(2.18), (12.2)) для λ = 1/2 и моделируется по формуле (12.5) (или (2.20)):ρ̃0 = −2 ln α0 ; α0 ∈ U (0, 1).В свою очередь, моделирующие формулы для двумерного единичного изотропного случайного вектора имеют вид:(1)(2)ω0 = cos 2πα0 ; α0 ∈ U (0, 1)(13.3)√(см. формулы (6.30) для случая ρ0 = R α1 = 1).Таким образом, в двумерном случае алгоритм 13.1 сводится к знаменитым (и наиболее употребимым) формулам Бокса – Мюллера(см., например, [9, 13]):pp(0,1)(0,1)ξ1= −2 ln α1 sin 2πα2 , ξ2= −2 ln α1 cos 2πα2 ; α1 , α2 ∈ U (0, 1),(13.4)дающим сразу пару выборочных значений стандартной гауссовской (нормальной) случайной величины ξ (0,1) .13.2.
Решение проблемы 6.1. Альтернативные алгоритмымоделирования d-мерного изотропного случайного вектора истандартной гауссовской (нормальной) случайной величины.Полученные нами формулы (13.4) помогают решить проблему 6.1 о построении алгоритма моделирования выборочного значения ξ 0 случайной точки ξ, равномерно распределенной в d-мерном шаре B (d,0,R) (см.соотношение (6.7)) для d ≥ 4 (в подразделе 6.4 для этого случая построить такой алгоритм с помощью перехода к гиперсферическим координатам не удалось).Из рассуждений подраздела 6.4 и доказательства утверждения 13.1следует, что моделируемый вектор ξ можно представить в видеω0 = sin 2πα0 ,ξ = ρ × ω,где ω – единичный изотропный случайный вектор, а случайная величина ρ (длина вектора), распределенная согласно степенной плотностиfρ (r) = drd−1 /Rd ; при этом величины ω и ρ независимы.Используя известное из теории вероятностей представление222(d)χ= ξ 1(0,1) + .
. . + ξ d(0,1) , несложно доказать обратное (по отношению к утверждению 13.1) утверждение.214УТВЕРЖДЕНИЕ 13.2 (см., например, [9, 24]). Если случайные величины η (1) , . . . , η (d) независимы и распределены нормально с параметра- ми (0, 1) (т. е. согласно плотности (13.2)), то вектор η = η (1) , . . . , η (d)является изотропным.Используя это утверждение, а также формулы Бокса – Мюллера(13.4) и аналог формулы (2.22), получаем искомый алгоритм, решающий проблему 6.1 для случая d ≥ 4.АЛГОРИТМ 13.2 (см., например, [9, 24]). Пусть d = 2k; k = 2, 3, ....Моделируем k пар значений по формулам Бокса – Мюллера (13.4)qq(i)(i)(1,i)(i)(i)(2,i)η0= −2 ln α1 sin α2 , η0= −2 ln α1 cos α2 ; i = 1, .