Оганесян Введение в физику тяжелых ионов 2008

Рис. 11.3. Иллюстрация метода обращения функции вероятности

Действительно, P(x < ξ < x + dx) = P{F (x) < γ < F (x + dx)} в силу непрерывности F (x). Так как γ — равномерно распределена в интервале (0, 1), то P(x < ξ < x + dx)= F (x + dx)− F (x)= P(x)dx .

Метод Неймана (алгоритм №3). Теорема. Пусть γ1 и γ2

— не-

зависимые случайные числа и

(b −a)

′

= a + γ1

′

ξ

, η = cγ2 .

Случайная величина ξ, определенная из условия

′

′

′

ξ = ξ , если η

< P(ξ ),

имеет плотность распределения P(x).

В соответствии с формулой полной вероятности:

′

′

′

=

P(ξ < z)= P{ξ < z | η < P(ξ )}

(11.4)

{

′

′

(

′

)}

{

′

(

′

)}

= P ξ < z,η < P

ξ

P η < P ξ .

Знаменатель дроби есть вероятность попадания точки

′

′

Q(ξ ,η ) под

кривую

y = P(x). Плотность

точки

Q

постоянна и

равна

c(b −a) −1 и

361

	b	P(x)	−1	−1
′	′				.
P{η < P(ξ )}= ∫dx ∫			c(b −a)	dy = c(b − a)
	a	0

В то же время
{	′
P ξ < z,
z	P(x)
= ∫dx ∫
a	0

η′ < P(ξ′)} =

c(b −a) −1 dy = c(b − a) −1 ∫z P(x)dx.

a

(11.5)

(11.6)

Окончательно, подставляя выражения (11.5) и (11.6) в (11.4) имеем:

P(ξ < z)= ∫z P(x)dx ,

a

что и требовалось доказать.

Рис. 11.4. Иллюстрация метода Неймана

Резюмируя, можно сказать, что метод Неймана — это метод отбора: из всех точек Q, равномерно распределенных в прямоугольной области (0, с, b, a), отбираются только те, которые лежат под

кривой y = P(x). Абсциссы этих точек и образуют выборку значений случайной величины, распределенной с плотностью P(x).

Моделирование многомерных случайных величин (алго-

ритм №4). Пусть случайные величины ξ1,...ξn в общем случае за-

362

висимы, и надо получить алгоритм моделирования n-мерной точки Q, распределенной с плотностью PQ . Представим PQ в виде:

PQ (x1,..., xn ) = P1 (x1 )P2 (x2 | x1 )P3 (x3 | x1, x2 )...Pn (xn | x1...xn−1 ).

Все входящие в это выражение плотности можно выразить через совместную плотность PQ :

P1 (x1 )= ∫...∫PQ dx1...dxn ,

P2 (x2 | x1) = ∫...∫PQ dx3...dxn P1 (x1 ).

Ниже рассмотрен частный случай моделирования многомерной случайной величины.

Алгоритм выбора случайного направления в пространстве

(алгоритм №5). Рассмотрим сначала алгоритм моделирования точки, равномерно распределенной в шаре, т. е.

P	(x, y, z)=	4	πR3	−1 .
Q
		3

Перейдем к сферическим координатам, соответствующим симметрии задачи:

x= r sin θcos j,

y= r sin θsin j,

z= r cosθ.

Вновых координатах шар превращается в параллелепипед:

0	≤ r < R,
0	≤ θ < π,
0	≤ ϕ < 2π.

Якобиан преобразования ∂(x, y, z)		∂(r,θ,ϕ)= r2 sin θ и плотность
распределения в новых координатах можно записать:
PQ (r,θ,ϕ)=	4 πR3 −1 r2 sin θ.
	3

Легко видеть, что это произведение трех плотностей:

PQ (r,θ,ϕ)= (3r2 R−3 )(2−1 sin θ)(2π)−1 ,

т. е. сферические координаты r, θ, φ точки Q независимы. В рамках метода обращения функции вероятности:

363

r		2	3dr		θ			ϕ
∫3r				= γ1 ,	∫sin θdθ		=1− γ2 ,	∫dϕ		= γ3. .
0	R				0	2		0	2π

Решение этих уравнений приводит к формулам:

r = R 3 γ1 , cosθ = 2γ2 −1, ϕ= 2πγ3.

Собственно алгоритм выбора случайного направления.

Пусть направление задается единичнымGвектором:

ω=G i ω1 + Gjω2 + kω3 ,

ω1 +ω2 +ω3 =1.

Необходимо отыскать случайный вектор ω, такой, чтобы для любого телесного угла Ω

G	(11.7)
P(ω Ω) = Ω 4π .

Если случайная точка Q равномерно распределена в шаре, то радиус-вектор этой точки обладает нужным свойством. Действительно, если Ω1 и Ω2 — два равных телесных угла, то равны и

объемы соответствующих шаровых секторов. В силу этого, вероятности попадания точки Q в эти сектора равны, что доказывает утверждение. Таким образом, формулы для выбора случайного направления:

cosθ = 2γ1 , ϕ= 2πγ2 .

Пример моделирования естественного процесса: задача о поглощении нейтронов. Пусть имеется выпуклая область G, не содержащая делящихся ядер. Считаем известными сечения рассеяния и поглощения

Σ(r, E)= Σs (r, E)+ Σa (r, E),

атакже индикатрису рассеяния ps (r′,Ω′;Ω), которая определяет

вероятность того, что нейтрон в точке r′, летевший в направлении Ω′ , в результате рассеяния полетит в направлении Ω. В точке r0 G расположен источник нейтронов с энергией E0 и равнове-

роятными направлениями Ω0 начальной скорости. Требуется вычислить вероятность pA того, что нейтрон, вылетевший из источ-

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

364

трона в область G, если он ее покинул. Алгоритм моделирования включает следующие шаги.

— Выберем случайное направление скорости нейтрона Ω0 при его вылете из источника (см. алгоритм №5);

— Разыграем для него длину свободного пробега ξ0 в соответ-

ствии с формулой:

ξ = −(1Σ)ln γ .

Она получается, если применить к экспоненциальному распределению, которому подчиняется величина ξ (алгоритм №2).

— Находим точку столкновения нейтрона:

r1 = r0 +ξ0Ω0 .

Если r1 G , то история нейтрона закончилась его вылетом из области G. Такому событию соответствует ηA = 0 .

— При r1 G разыгрываем судьбу нейтрона, т. е. получаем кон-

кретное значение случайной величины, имеющей два значения – s (рассеяние) и а (поглощение), используя алгоритм 1.

—Если нейтрон поглотился, полагаем ηA =1 .

—Если нейтрон рассеялся, то в соответствии с индикатрисой рассеяния ps (r′,Ω′;Ω) разыгрываем новое направление скорости

Ω= Ω1 (алгоритм №4 и далее алгоритм №2 или алгоритм №3) и

переходим на блок программы, где вычисляется длина свободного пробега.

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

N

pˆ A =1 N ∑(ηA )s .

s=1

Или, в другом виде,

pˆ A = NA N ,

где N — полное число разыгранных траекторий, а NA — количест-

во траекторий, завершившихся поглощением.

Дисперсия полученной оценки совпадает с дисперсией величины ηA :

365

D(ηA )= pA (1− pA ).

Рис. 11.5. К задаче о поглощении нейтрона:

вылет из области G маркируется значением η= 0 , поглощение — η=1

11.3. Проверка статистических гипотез

Пусть имеется выборка объема n: x1 , ..., xn , и есть основания предполагать, что она взята из генеральной совокупности, распределенной с плотностью P(x) (будем считать, что параметры рас-

пределения известны). Требуется проверить эту статистическую гипотезу.

Представим выборку в виде полигона накопленных частот (рис. 11.6), т. е. укажем количество элементов выборки, попавших в интервал xi (другими словами, абсолютные частоты ν(xi )). Кро-

ме того, на этом же рисунке отложим значения ожидаемого числа событий в соответствующем интервале, если проверяемая гипотеза верна:

npi = n ∫ P(u)du .

xi

366

Рис. 11.6. Полигон накопленных частот для исходной выборки (сплошные линии) и ожидаемое число событий в интервале x , если выборка взята из генеральной

совокупности, распределенной с плотностью P(x)

В излагаемом подходе дело, по существу, сводится к решению вопроса — похожи ли друг на друга две гистограммы, показанные на рис. 11.6. Чтобы суждение было объективным, необходимо выбрать меру близости между сравниваемыми объектами. Естественной мерой близости представляется величина:

r

S0 = ∑(νi −npi ),

i=1

однако отклонения разных знаков могут скомпенсировать друг друга, так что S0 ~ 0 не будет свидетельством близости сравнивае-

мых гистограмм. Отмеченная неопределенность снимается, если в качестве меры выбрать:

r

S = ∑(νi −npi )2 ,

i=1

где r — число интервалов разбиения, т. е. число столбиков в гистограммах на рис. 11.6.

367

Покажем, что при некоторых условиях немного измененная мера S распределена по χ2 . По определению, распределение χ2 имеет величина:

f
χ2 = ∑ξ2 ,	(11.8)

i=1

где ξ — нормально распределенная со средним «0» и дисперсией «1» случайная величина, f — число независимых слагаемых в сумме, называемое числом степеней свободы.

Определим S1 как:

r
S1 = ∑(νi − npi )2 npi .	(11.9)

i=1

Сформулируем условия, при которых статистика S1 (называемая проверочной статистикой или критерием) распределена по χ2 .

1. Все νi должны быть распределены нормально. Что для этого

требуется? Если рассматривать один отдельный столбик гистограммы, представляющей экспериментальные данные, его высота νi распределена по биному. Действительно, для любого элемента

исходной выборки есть только две возможности: попасть в данный интервал xi либо в любой другой. Такое «испытание» происходит

n раз, число «успехов» равно νi . Как известно, биномиальное распределение может быть аппроксимировано нормальным, если среднее число успехов, в данном случае, νi >>1. На практике это означает, что для любого i должно быть νi ≥10 .

2.	Каждое слагаемое в сумме (11.9) должно иметь нулевое сред-
нее.	Посмотрим,	что вытекает из	такого требования. Если
νi − npi = 0 , νi		= npi , т. е. условие выполняется только для за-
ведомо правильной гипотезы.
3.	Дисперсия каждого слагаемого (без квадрата — см. формулу
(11.8)) должна быть равна единице. Вычислим эту дисперсию:
		D (νi −npi ) npi = D(νi ) npi .
Она равна 1, если D(νi )= npi = νi			. Так будет, если νi распре-

делены по закону Пуассона. Де-факто, как отмечалось выше, νi

368

распределены по биномиальному закону. В тоже время известно, что при выполнении условий νi n → 0 , n → ∞ , биномиальное рас-

пределение можно аппроксимировать распределением Пуассона. В свою очередь, указанные требования выполняются, если объем исходной выборки n велик и велико число интервалов разбиения r.

Таким образом, практическая рекомендация состоит в том, что при выполнении условия νi ≥10 число интервалов разбиения долж-

но быть максимальным.

Резюмируя, можно сказать, что только при выполнении всех трех рассмотренных выше условий статистика S1 (формула (11.9))

распределена по χ2 с r −1 степенью свободы. Условие

r

∑νi = n уменьшило на 1 число независимых слагаемых в выра-

i=1

жении (11.9) и, следовательно, число степеней свободы распределения χ2 .

Расчетное число степеней свободы изменится, если исходная выборка взята из генеральной совокупности с плотностью P(x,αG),

где αG — неизвестный вектор параметров. Он может быть оценен по той же выборке, например, методом максимального правдоподобия. Уравнения правдоподобия (скажем, s штук — по числу неизвестных параметров), уменьшат число независимых слагаемых в формуле (11.9) также на s штук. Окончательно, формула для числа степеней свободы f в общем случае записывается так:

f = r − s −1.

Итак, мы знаем, что проверочная статистика (11.9) распределена по χ2 при выполнении условий 1-3, причем второе из условий вы-

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

сти с плотностью P(x,αG) и она действительно взята из этой гене-

ральной совокупности). Как воспользоваться распределением проверочной статистики (11.9), известным только для заведомо пра-

вильной гипотезы (распределение χ2 ), в реальной задаче, когда

истинное распределение генеральной совокупности неизвестно и мы пытаемся его «угадать»? Подход состоит в следующем.

369

Выделим область значений критерия, вероятность попадания в которую α мала для заведомо правильной гипотезы (типичные значения α = 0.01÷0.05 ) (рис. 11.7). Эта область называется критической (правосторонняя критическая область). Отрезок от начала ко-

ординат до точки χкрит2 называется областью допустимых значений. Пусть найденное по выборке значение критерия равно χexp2 (фор-

мула (11.9)). Если гипотеза верна, а χexp2 > χкрит2 , то произошло маловероятное событие, в силу чего проверяемую гипотезу следует отклонить. И наоборот — если χexp2 < χкрит2 — гипотеза не противоречит экспериментальным данным.

.

Рис. 11.7. Распределение χ2 . Справа от точки χкрит2 лежит правосторонняя критическая область. При этом P(χ2 > χкрит2 )= α .

Из рис. 11.7 видно, что, в принципе, область маловероятных значений критерия можно выделить и в окрестности 0. Как интер-

претировать попадание χexp2 в эту левостороннюю критическую

область? Получается, что эксперимент «слишком хорошо» соответствует ожиданиям (теоретической модели). Во избежание такого искусственного согласия должно быть обеспечено ns >>1. Для

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

370