Метод Монте-Карло. Статистические методы в ядерном эксперименте

Метод Монте-Карло Лекция №4 Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. 1 2 3 4 Введение Моделирование дискретных случайных величин Моделирование непрерывных случайных величин Интегральный метод Метод браковки Неймана Метод существенной выборки Метод суперпозиции Метод частичного аналитического интегрирования Вычисление интегралов методом Монте-Карло Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 1/11 Метод Монте-Карло (I) Метод Монте-Карло – численный метод решения сложных математических задач путем генерирования случайных величин. Позволяет моделировать сложные процессы, на которые влияет множество случайных факторов: рождение и распад частиц, прохождение частиц через вещество и.т.п. Метод Монте-Карло может быть применен и для решения некоторых задач, вовсе не имеющих случайных составляющих, например, вычисление определенных (многомерных) интегралов от сложных функций, решение систем дифференциальных уравнений в частных производных со сложными граничными условиями. Так, например, при равномерном выбрасывании точек N0 раз в объёме гиперкуба (Vhypercube ), вероятность того что точка попадёт (N раз) в область, ограниченную заданными гиперповерхностями, пропорциональна объёму (Vdomain ) области: s „ « N Vhypercube N Vhypercube , σ(Vdomain ) = N 1 − , N0 N0 N0 v„ v„ « « u u V u u hypercube u 1 − NN u Vdomain − 1 t t σ(Vdomain ) 1 = = ∼ p Vdomain N N0 N0 Vdomain = Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 2/11 Метод Монте-Карло (II) Являясь статистическим по своей природе, метод Монте-Карло дает ответ, которому также присуща некоторая статистическая неопределенность - точность результата. При этом точность улучшается с увеличением числа статистических испытаний (событий). Ключевым моментом в применимости метода Монте-Карло является возможность генерации (псевдо)случайных чисел, равномерно распределённых на интервале (0, 1). По способу получения последовательности (псевдо)случайных чисел генераторы делятся на следующие типы: Физические Табличные Алгоритмические Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 3/11 Моделирование дискретных случайных величин Отложим на единичном интервале (0, 1) отрезки P1 , P2 , P3 , ..., Pn . Точка, соответствующая выпавшему равномерно распределённому на этом интервале случайному числу, укажет на один из отрезков Pi , i = 1 ÷ n. Это означает появление соответствующего значения Xi дискретной случайной величины. X1 P1 X2 P2 P1 P3 P1+P2 Xn X3 Pn P1+P2 + P3 1 r r= γ(0,1) Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 4/11 Моделирование непрерывных случайных величин (I) Интегральный метод: Случайная величина x = x(γ), являющаяся решением интегрального уравнения: +∞ Zx Z γ= f (t)dt, f (t)dt = 1, −∞ −∞ обладает ф.п.в. равной f (x). Замечание: Есть случайная величина y с ф.п.в. g(y ) = γ(0, 1), сделаем замену y= Rx −∞ f (t)dt и найдём ф.п.в. для x : g(y (x)) · dy dx = 1 · f (x) = f (x) Пример 1: Случайная величина x = γ(a, b) = (b − a)γ + a Пример 2: Случайная величина x = −λ ln γ распределена по экспоненциальному закону λ1 e−x/λ p Пример 3: Случайная величина x = −2 ln γ1 cos(2πγ2 ) (γ1 и γ2 – 2 независимые случайные величины) распределена по N(0, 1) = √1 e−x /2 (т.н. 2π преобразование Бокса-Мюллера) Пример 4: Если ф.п.в. x равна N(0, 1), то y = σx + µ распределена по N(µ, σ) Пример 5: Случайная величина x = ctg γπ распределена по π1 x 21+1 Пример 6: Сгенерировать интегральным методом случайные величины с ф.п.в.: f (x) = β/x 1−β , x ∈ (0, 1) и g(y ) = β/y 1+β , y ∈ (1, +∞), 0 < β < 1 Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 5/11 Моделирование непрерывных случайных величин (II) Далеко не всегда можно воспользоваться интегральным методом (решить аналитически интегральное уравнение), если ф.п.в. f (x) определена на конечном интервале [a, b] и ограничена сверху f (x) < C, ∀x ∈ [a, b], то используют метод браковки по Нейману. Метод браковки по Нейману: Разыгрываем пару случайных величин (равномерно в прямоугольнике (b − a) × C ): ξ = (b − a)γ1 + a, η = Cγ2 Если η ≤ f (ξ), то сгенерированное x = ξ принимается, в противном случае возвращаемся к предыдущему пункту и повторяем всё сначала C η’ f(ξ ) η f(x) a ξ b Замечание 1: Функцию f (x) не нужно нормировать в методе браковки. Замечание 2: Эффективность метода зависит от того, какую долю от площади прямоугольника (b − a) × C занимает площадь под кривой. Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 6/11 Моделирование непрерывных случайных величин (III) Если невозможно воспользоваться интегральным методом, а также область определения ф.п.в. бесконечна или ф.п.в. неограничена (но интегрируема), используют метод существенной выборки, который является комбинацией интегрального метода и метода браковки. Метод существенной выборки: Заключается в поиске/подборе т.н. мажорирующей функции (мажоранты), M(x), для заданной ф.п.в. f (x), такой что f (x) ≤ M(x) ∀x и для нормированной M(x) применим интегральный метод (если нормированную M(x) рассматривать как R некоторую простую ф.п.в.). Пусть M 0 (x) = const · M(x), такая что M 0 (x)dx = 1. Тогда процедура: Генерируем с.в. ξ согласно ф.п.в. M 0 (x) интегральным методом Генерируем η = M(ξ)γ Если η < f (ξ), то сгенерированное x = ξ принимается, в противном случае возвращаемся к предыдущему пункту и повторяем всё сначала M(ξ ) η’ f(ξ) η M(x) f(x) a ξ b Замечание: Метод существенной выборки применяется также для повышения эффективности генерирования с.в. Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 7/11 Моделирование непрерывных случайных величин (IV) Пример 1: Сгенерировать с.в., распределённую по нормальному закону (N(0, 1)) используя подходящую мажоранту. √ Пример 2: Сгенерировать с.в., распределённую по закону ∼ 1/ sin x , подобрав подходящую мажоранту. Пример 3: Физическая величина, x , распределена согласно ф.п.в. f (x). В реальном эксперименте детектор имеет конечное разрешение, поэтому ф.п.в. для измеренного значения физической величины, x 0 , будет: Z (x 0 −x)2 − 1 2 g(x 0 ) = f (x) √ e 2σx dx 2πσx Если уже есть выборка величины x объёмом n (x1 , ..., xn ), как из неё получить выборку измеренной величины x 0 объёмом n ? Пример 4: Вероятность излучения фотона электроном или позитроном пучка β пропорциональна D(x) = x 1−β h(x), где x – доля энергии пучка, уносимая фотоном, h(x) = 1 + ln(1−x) x)2 ) x 3 4β − β2 12 „ « L 3 + π2 − +x −6 , β = α π (L 47 8 « − 1 1−β (2 2x − 1), L = ln s m2 − x) + β 1−β 8x „ 4(x − 2) ln x − (1 + 3(1 − , s = 4E , E – энергия пучка. Таким 2 образом видимое сечение реакции с учётом т.н. радпоправок равно: Z1 Z1 D(x1 )D(x2 )σ(s(1 − x1 )(1 − x2 ))dx1 dx2 σvis = Сгенерировать с.в. x согласно ф.п.в. D(x)/ R1 D(x)dx . Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 8/11 Моделирование непрерывных случайных величин (V) Метод суперпозиции позволяет генерировать с.в. чьи ф.п.в. составлены из нескольких слагаемых с разным поведением: Z f (x) = c1 g1 (x) + c2 g2 (x), g1,2 (x)dx = 1, c1 + c2 = 1 Введём вспомогательную дискретную с.в. r : P(r = 1) = c1 , P(r = 2) = c2 . Для генерирования исходного распределения сначала разыгрываем r , а затем gr (x) (т.е. либо g1 (x) либо g2 (x)). Достоинством метода является то, что применимость интегрального метода для g1 (x) и g2 (x) в отдельности может быть существенно проще чем для исходной ф.п.в. f (x). 5 Пример: Сгенерировать с.в. с ф.п.в. f (x) = 12 (1 + x 4 ), −1 < x < 1. Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 9/11 Моделирование непрерывных случайных величин (VI) Метод частичного аналитического интегрирования применяется в ситуации, когда нужно сгенерировать многомерную с.в. с некоторой ф.п.в. f (~x |~ α), ~x = (x1 , ..., xn ), α ~ = (α1 , ..., αk ) – набор параметров распределения. Если возможно аналитическое интегрирование f (~x |~ α) по одной или нескольким переменным (x1 , ..., xm , m < n), то генерирование ~x можно проводить в два этапа: генерируем набор R с.в. (xm+1 , ..., xn ) согласно ф.п.в. g(xm+1 , ..., xn ) = f (~x |~ α)dx1 ...dxm генерируем набор оставшихся с.в. (x1 , ..., xm ) согласно f (~x |~ α) = f (x1 , ..., xm |xm+1 , ..., xn , α ~ ), где уже сгенерированные с.в. (xm+1 , ..., xn ) рассматриваются как дополнительные параметры распределения. Замечание: Если m = n − 1, то генерирование ~x распадается на серию из n розыгрышей одномерных с.в. Пример: f (x, y |α) = (α + 1)(α + 2)(1 − x − y )α H(1 − x − y ), 0 < x, y < 1 Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. Лекция №4 10/11 Вычисление интегралов методом Монте-Карло I. Метод Неймана (или геометрический метод) Интеграл от функции (задающей профиль ф.п.в.) рассчитывается как q I = V ε, σI = V ε(1−ε) , ε = NN1 ,, где N0 – полное число испытаний при N0 генерировании выборки из N1 случайных величин по методу Неймана, V – объём гиперкуба (в одномерном случае (b − a) × C ), или интеграл от мажоранты. II. Метод выборки (или обычный метод) I= Rb f (t)dt ≈ a b−a n n P f (xi ), где (x1 , ..., xn ) – выборка с.в. с ф.п.в. i=1 γ(a, b). В общем случае, если есть выборка (x1 , ..., xn ) с.в. с ф.п.в. g(x), то: n R f (t(y )) R R f (t) P f (xi ) f (x) I = f (t)dt = g(t) g(t)dt = g(t(y dy ≈ n1 = h g(x) i , )) g(xi ) i=1 y= Rt g(x) g(z)dz , y = γ(0, 1), x = t(y )–распределена по g(x). −∞ Причём ошибка s интеграла это просто ошибка среднего, т.е. q n P 2 f (xi ) 1 σI = sn = n(n−1) ( g(x − I)2 ∼ √1n i) i=1 Пример: Вычислить интеграл методом Монте-Карло I = Rπ Статистические методы в ядерном эксперименте, ИЯФ, 2019 г. √dx sin x Лекция №4 11/11