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

Статистический анализ данных

Дмитрий В. Наумов (ОИЯИ)

Зачем нужен Монте-Карло

Когда аналитики недостаточно

  • Распределение результата может быть неизвестно в замкнутой форме.
  • Преобразование измеряемых величин может быть сильно нелинейным.
  • Модель может включать пороги, отборы и физические границы.
  • Экспериментальная установка может иметь сложный многомерный отклик.
  • Интеграл может иметь большую размерность.

Метод Монте-Карло заменяет сложное аналитическое вычисление статистикой большого числа искусственно сгенерированных реализаций.

Основная идея

Пусть требуется узнать свойства случайной величины \[ Y=g(X), \qquad X\sim f_X(x). \]

Алгоритм Монте-Карло:

  1. Сгенерировать независимые значения \[ X_1,\ldots,X_N\sim f_X. \]
  2. Для каждого вычислить \[ Y_k=g(X_k). \]
  3. Использовать выборку \(Y_1,\ldots,Y_N\) для оценки распределения, среднего, дисперсии или вероятности события.

Три применения одного метода

Численное интегрирование

\[ \int h(x)f(x)\,dx = \mathbb E[h(X)]. \]

Распространение ошибок

\[ \mathbf X\sim f(\mathbf x) \quad\Longrightarrow\quad Y=g(\mathbf X). \]

Псевдоэксперименты

Повторить весь эксперимент много раз и изучить распределение результата.

Что Монте-Карло не делает

  • Не исправляет неверную физическую или статистическую модель.
  • Не устраняет систематические неопределённости автоматически.
  • Не превращает малую реальную выборку в большую.
  • Не даёт точного ответа при конечном числе симуляций.
  • Не заменяет проверку аналитических предельных случаев.

Монте-Карло точно решает ту задачу, которую мы запрограммировали, а не обязательно ту задачу, которую хотели решить.

Генерация случайных величин

Псевдослучайные числа

Компьютерный генератор создаёт детерминированную последовательность \[ U_1,U_2,\ldots, \qquad U_k\in[0,1), \] которая статистически похожа на независимую равномерную выборку.

  • Начальное состояние генератора задаётся специальным начальным числом.
  • Одинаковое начальное число воспроизводит ту же последовательность.
  • Разные вычислительные потоки должны использовать независимые последовательности.
  • Для научной воспроизводимости начальное число и алгоритм генератора необходимо сохранять.

Равномерное распределение как исходный материал

Большинство генераторов непосредственно создаёт \[ U\sim\mathrm U(0,1). \]

Из равномерных чисел можно получить другие распределения:

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

Метод обратной функции распределения

Пусть \[ U\sim\mathrm U(0,1), \] а \(F_X(x)\) — функция распределения требуемой случайной величины.

Тогда \[ \boxed{ X=F_X^{-1}(U) } \] имеет распределение \(F_X\).

Действительно, \[ P(X\le x) = P(F_X^{-1}(U)\le x) = P(U\le F_X(x)) = F_X(x). \]

Пример: экспоненциальное распределение

  • Для \[ f_X(x)=\lambda e^{-\lambda x}, \qquad x\ge0, \]

  • Функция распределения: \[ F_X(x)=\int\limits_{0}^{x}\, dx\, \lambda e^{-\lambda x} = 1-e^{-\lambda x}. \]

  • Следовательно, \[ \boxed{ X = -\frac{1}{\lambda}\ln(1-U). } \]

  • Поскольку \(1-U\) также равномерно, часто используют \(X=-\ln U/\lambda\).

Метод отбора

Пусть требуется генерировать плотность \(f(x)\), для которой \[ f(x)\le M g(x), \] а из \(g(x)\) генерировать легко.

  1. Сгенерировать кандидата \(X\sim g(x)\).
  2. Сгенерировать \(U\sim\mathrm U(0,1)\).
  3. Принять \(X\), если \[ U\le\frac{f(X)}{Mg(X)}. \]
  4. Иначе повторить попытку.
  • Средняя доля принятых событий равна \(1/M\), если \(f\) и \(g\) нормированы.

Метод отбора. Пример

  • Хотим генерировать двумерное распределение, ограниченное: \[ f(x)=6x(1-x), \qquad 0\le x\le1, \]
  • Генерируем точки равномерно под прямоугольником высоты \(M=1.5\).
  • Используем метод отбора

Интерактивный пример

Выбор метода генерации

Метод Преимущество Ограничение
Обратная функция распределения каждое число используется нужна вычислимая \(F^{-1}\)
Отбор простой и универсальный может быть неэффективным
Специальный алгоритм высокая скорость отдельная реализация для каждого семейства
Марковские цепи сложные многомерные плотности коррелированные выборки и диагностика сходимости

Интегрирование методом Монте-Карло

Интеграл как математическое ожидание

Пусть \[ I=\int_a^b h(x)\,dx. \]

Для \(U\sim\mathrm U(a,b)\): \[ \mathbb E[h(U)] = \frac{1}{b-a} \int_a^b h(x)\,dx. \]

Поэтому \[ \boxed{ \widehat{I}_N = \frac{b-a}{N} \sum_{k=1}^{N}h(U_k). } \]

Оценка ошибки вычисления интеграла методом Монте-Карло

Если \[ \operatorname{Var}[h(U)]=\sigma_h^2, \] то \[ \operatorname{Var}(\widehat I_N) = \frac{(b-a)^2}{N^2}\operatorname{Var}\left[\sum_{k=1}^N h(U_k)\right] = \frac{(b-a)^2\sigma_h^2}{N}. \]

Следовательно, \[ \boxed{ \sigma_{\widehat I_N} \propto \frac{1}{\sqrt N}. } \]

  • Чтобы уменьшить ошибку в 10 раз, требуется в 100 раз больше точек.
  • Сходимость медленная, но почти не зависит от размерности интеграла.

Пример: площадь четверти круга

Равномерно генерируем точки в единичном квадрате: \[ (X,Y)\sim\mathrm U([0,1]^2). \]

Вероятность попасть внутрь четверти круга: \[ p=P(X^2+Y^2\le1)=\frac{\pi}{4}. \]

Если \(N_{\mathrm{in}}\) точек попали внутрь, то \[ \boxed{ \widehat\pi = 4\frac{N_{\mathrm{in}}}{N}. } \]

Интерактив: оценка числа \(\pi\)

Почему Монте-Карло полезен в большой размерности

Для регулярной сетки с \(m\) точками по каждой координате число вычислений равно \[ N_{\text{сетка}}=m^d. \]

В Монте-Карло статистическая ошибка обычно ведёт себя как \[ \sigma_{\widehat I}\propto N^{-1/2}, \] почти независимо от размерности \(d\).

Монте-Карло особенно полезен для многомерных интегралов, где регулярная сетка страдает от «проклятия размерности».

Не все точки одинаково полезны

Если основной вклад в интеграл приходит из малой области, равномерная генерация тратит большинство точек впустую.

Вместо \[ I=\int h(x)f(x)\,dx \] можно генерировать \(X\sim q(x)\): \[ \boxed{ I = \mathbb E_q\!\left[ h(X)\frac{f(X)}{q(X)} \right]. } \]

Это называется выборкой по важности.

Веса событий

При выборке по важности каждое событие получает вес \[ w_k=\frac{f(X_k)}{q(X_k)}. \]

Оценка интеграла: \[ \widehat I = \frac1N\sum_{k=1}^Nw_kh(X_k). \]

  • Хорошая \(q(x)\) чаще генерирует важные области.
  • Сильно различающиеся веса увеличивают дисперсию.
  • Вес не превращает одно событие в большое число независимых событий.

Почему сетка быстро проигрывает

В большой размерности даже грубая сетка быстро становится дорогой. Монте-Карло платит не за полную сетку, а за число реально брошенных точек.

Редкая область: равномерная генерация тратит точки

Выборка по важности: бросаем точки туда, где вклад

Веса: одно большое событие не заменяет статистику

При выборке по важности мы меняем способ генерации точек.

Событие получает вес, который компенсирует эту замену.

Хорошая выборка по важности делает веса умеренными. Плохая — создаёт редкие огромные веса.

хорошие веса плохие веса много сопоставимых вкладов дисперсия мала оценку держат редкие точки дисперсия велика

VEGAS: адаптивная выборка по важности

Идея VEGAS: не угадывать хорошее распределение заранее.

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

На следующих итерациях выборка становится ближе к важностной.

до адаптации после адаптации равные ячейки мелкие ячейки около пика

Тестовый интеграл: узкий пик в большой размерности

Возьмём нормированный гауссов пик в единичном гиперкубе:

\[ f(\mathbf x) = \prod_{i=1}^d \frac{1}{\sqrt{2\pi}\sigma} \exp\left[-\frac{(x_i-\mu_i)^2}{2\sigma^2}\right], \qquad 0\le x_i\le1 . \]

Интеграл почти равен единице, но вклад сосредоточен в малой области.

Этот пример контролируем: точное значение получается как произведение одномерных интегралов через erf.

Простой Монте-Карло

Равномерная генерация честная, но не умная: она почти не знает, где находится пик.

Для этого примера точное значение: \[ I_{\text{точн.}}=0.9978196359 . \]

Типичный запуск:

\(N\) простой Монте-Карло
\(10^4\) \(0.946\pm0.48\)
\(10^5\) \(0.864\pm0.13\)
\(10^6\) \(1.001\pm0.061\) 
def plain_mc(N, seed=12345):
    rng = np.random.default_rng(seed)
    x = rng.random((N, d))
    y = f_batch(x)
    mean = np.mean(y)
    err = np.std(y, ddof=1) / np.sqrt(N)
    return mean, err

for N in [10_000, 100_000, 1_000_000]:
    val, err = plain_mc(N)
    print(f"простой МК, N={N:>8}: {val:.6g} ± {err:.2g}")

Тензорная квадратура Гаусса–Лежандра

В малой размерности квадратуры прекрасны. В большой размерности проблема не в идее квадратуры, а в росте числа узлов как \(m^d\).

узлов на координату всего узлов результат
3 6 561 0.589
5 390 625 0.722
7 5 764 801 0.963

Даже миллионы регулярных узлов ещё чувствительны к тому, как сетка попала относительно узкого пика.

from numpy.polynomial.legendre import leggauss
import itertools

def tensor_gauss(m):
    nodes, weights = leggauss(m)
    nodes = 0.5 * (nodes + 1.0)
    weights = 0.5 * weights

    total = 0.0
    for idx in itertools.product(range(m), repeat=d):
        x = np.array([[nodes[i] for i in idx]])
        w = np.prod([weights[i] for i in idx])
        total += w * f_batch(x)[0]

    return total, m**d

for m in [3, 5, 7]:
    val, neval = tensor_gauss(m)
    print(f"m={m:>2}, вызовов={neval:>9}: {val:.6g}")

VEGAS в Python

VEGAS не отменяет \(1/\sqrt N\). Он уменьшает коэффициент перед этой ошибкой, если смог найти важные области.

В этом запуске после адаптации: \[ \widehat I_{\mathrm{VEGAS}}=0.99787\pm0.00053, \qquad Q\simeq0.43 . \]

Это меньше \(10^{-3}\) относительной статистической ошибки при меньшем числе вычислений, чем у простого Монте-Карло с \(10^6\) точек.

import runpy

demo = runpy.run_path("../../shared/notebooks/mc_vegas_demo.py")
exact = demo["exact_integral"]()
mc_val, mc_err = demo["plain_mc"](1_000_000)
grid_val, grid_calls = demo["tensor_gauss"](5)
vegas_val = demo["vegas_result"]()

print(f"точное       {exact:.10f}")
print(f"простой МК   {mc_val:.6f} +- {mc_err:.3f}")
print(f"Гаусс m=5    {grid_val:.6f}  ({grid_calls} вызовов)")
print(f"VEGAS      {vegas_val}  Q={vegas_val.Q:.2f}")
точное       0.9978196359
простой МК   1.000935 +- 0.061
Гаусс m=5    0.722156  (390625 вызовов)
VEGAS      0.99787(53)  Q=0.43

Полный воспроизводимый пример: ноутбук mc_vegas_notebook.

Что сравниваем

метод что делает сильная сторона слабое место
сеточная квадратура регулярно покрывает гиперкуб точна в малой размерности число узлов растёт как \(m^d\)
простой Монте-Карло бросает равномерные точки ошибка \(\sim N^{-1/2}\) много точек мимо редкого вклада
выборка по важности бросает точки чаще в важные области меньше дисперсия плохие веса портят оценку
VEGAS строит важностную выборку адаптивно сам ищет пики труднее для диагональных и сложных структур

Главный вывод: выигрыш VEGAS — не магическое изменение закона \(1/\sqrt N\), а уменьшение дисперсии интегранда после адаптации.

Когда VEGAS надо проверять особенно внимательно

  • Очень узкий пик можно просто не увидеть на начальных итерациях.
  • Несколько разнесённых мод могут привести к адаптации вокруг одной из них.
  • Диагональные гребни и сильные корреляции плохо описываются произведением одномерных сеток.
  • Нужно смотреть не только центральное значение, но и итерации, \(Q\)-значение и независимые запуски с разными начальными числами генератора.

Практическая стратегия: сначала простой Монте-Карло или грубый контрольный расчёт, затем VEGAS, затем проверка устойчивости результата.

Псевдоэксперименты

Что такое псевдоэксперимент

Псевдоэксперимент — одна искусственная реализация полного измерения:

  1. Сгенерировать истинные физические события.
  2. Смоделировать отклик детектора.
  3. Применить реконструкцию и отбор.
  4. Выполнить тот же анализ, что и для данных.
  5. Сохранить итоговую статистику или оценку параметра.

Повторение этой процедуры строит распределение возможных результатов эксперимента.

Зачем нужны псевдоэксперименты

  • Проверить смещение и дисперсию оценки параметра.
  • Проверить покрытие доверительного интервала.
  • Найти распределение статистики теста.
  • Оценить ожидаемую чувствительность эксперимента.
  • Проверить устойчивость подгонки.
  • Исследовать влияние систематических неопределённостей.

Пример: счётный эксперимент

Пусть ожидается \[ N\sim\operatorname{Pois}(s+b), \] где \(s\) — сигнал, а \(b\) — известный фон.

Простейшая оценка сигнала: \[ \widehat s=N-b. \]

Она несмещённая: \[ \mathbb E[\widehat s]=s, \] но отдельный эксперимент может дать даже \(\widehat s<0\).

Интерактив: ансамбль счётных экспериментов

Ансамбль и реальные данные

  • Реальный эксперимент даёт одну точку из распределения возможных результатов.
  • Псевдоэксперименты строят это распределение при заданной модели.
  • Они позволяют ответить: «что произошло бы, если повторить эксперимент?»
  • Но ответ всегда условен относительно выбранной модели и её параметров.

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

Моделирование эксперимента

Цепочка моделирования

В физическом анализе Монте-Карло часто проходит несколько уровней:

  1. Генерация процесса: частицы и кинематика.
  2. Взаимодействие с веществом: каскады, потери энергии, вторичные частицы.
  3. Отклик детектора: сигналы сенсоров и электроники.
  4. Реконструкция: восстановленные объекты и параметры.
  5. Отбор: триггеры, критерии качества и анализ.

На каждом уровне модель может вносить систематическую неопределённость.

Моделирование Daya Bay и JUNO в Opticks. Автор: Simon Blyth.

Эффективность отбора

Если из \(N_{\mathrm{gen}}\) сгенерированных событий прошло отбор \(N_{\mathrm{sel}}\), то \[ \widehat\varepsilon = \frac{N_{\mathrm{sel}}}{N_{\mathrm{gen}}}. \]

При независимых событиях \[ N_{\mathrm{sel}}\sim\operatorname{Bin}(N_{\mathrm{gen}},\varepsilon), \] поэтому \[ \boxed{ \sigma_{\widehat\varepsilon} \approx \sqrt{ \frac{\widehat\varepsilon(1-\widehat\varepsilon)} {N_{\mathrm{gen}}} }. } \]

Когда биномиальная ошибка ломается

Если все события прошли отбор, \(N_{\mathrm{sel}}=N_{\mathrm{gen}}\), то наивная формула даёт \[ \widehat\varepsilon=1, \qquad \sigma_{\widehat\varepsilon}=0. \]

Это не означает, что эффективность известна точно. Это означает только, что в данной выборке Монте-Карло не было ни одного отказа.

  • Точечная оценка \(\widehat\varepsilon=N_{\mathrm{sel}}/N_{\mathrm{gen}}\) остаётся разумной как оценка методом максимального правдоподобия.
  • На границах особенно ясно, что вместо дисперсии нужно использовать интервал. Симметричная запись \(\widehat\varepsilon\pm\sigma\) может вывести за границу определения.
  • Частотное решение: интервал Уилсона или интервал Клоппера–Пирсона; численно их можно посчитать в аплетах.
  • Байесовское решение: бета-биномиальное апостериорное распределение и байесовский интервал. Это будет отдельная тема.

Матрица отклика

Истинное значение может быть восстановлено в другом бине: \[ R_{ji} = P(\text{восстановленный бин }j\mid\text{истинный бин }i). \]

Если истинный спектр равен \(\mathbf t\), то ожидаемый восстановленный спектр: \[ \boxed{ \mathbf r=R\mathbf t. } \]

  • Диагональные элементы описывают правильную реконструкцию.
  • Внедиагональные элементы описывают миграции.
  • Потери эффективности уменьшают сумму элементов столбца.
  • Матрица отклика оценивается по Монте-Карло.

Ограниченная статистика моделирования

Монте-Карло имеет собственные статистические флуктуации.

Для взвешенной выборки эффективный размер: \[ \boxed{ N_{\mathrm{eff}} = \frac{\left(\sum_k w_k\right)^2} {\sum_k w_k^2}. } \]

  • Если все веса равны, \(N_{\mathrm{eff}}=N\).
  • Несколько больших весов могут сделать \(N_{\mathrm{eff}}\ll N\).
  • Большое число сгенерированных событий не гарантирует малую ошибку после отбора и перевзвешивания.

Как работать с Монте-Карло

Проверка генератора и реализации

  • Сравнить простые распределения с аналитическими результатами.
  • Проверить средние, дисперсии, корреляции и нормировку.
  • Проверить предельные случаи параметров.
  • Сравнить несколько независимых реализаций.
  • Зафиксировать начальное состояние генератора для воспроизведения ошибки.
  • Отдельно проверить каждый этап цепочки моделирования.

Проверка сходимости

Результат необходимо изучить как функцию числа симуляций \(N\).

  • Статистическая ошибка должна убывать примерно как \(1/\sqrt N\).
  • Независимые запуски должны быть совместимы в пределах ожидаемой ошибки.
  • Редкие области распределения требуют существенно большей статистики.
  • Для взвешенных событий следует контролировать \(N_{\mathrm{eff}}\).

Смещение и дисперсия оценки

Для оценки \(\widehat\theta\) псевдоэксперименты позволяют вычислить \[ \operatorname{bias}(\widehat\theta) = \mathbb E[\widehat\theta]-\theta, \] \[ \operatorname{Var}(\widehat\theta) = \mathbb E\!\left[ (\widehat\theta-\mathbb E[\widehat\theta])^2 \right]. \]

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

Разделять два источника ошибки

  • Физическая статистическая неопределённость: разброс возможных результатов реального эксперимента.
  • Численная ошибка Монте-Карло: неточность оценки этого разброса из-за конечного числа симуляций.

Например, увеличение числа псевдоэкспериментов точнее определяет распределение \(\widehat\theta\), но не уменьшает ожидаемую ошибку одного реального измерения.

Воспроизводимый рабочий процесс

  1. Явно записать вероятностную модель.
  2. Сохранить конфигурацию и версии программ.
  3. Сохранить начальное состояние генератора или правило формирования независимых начальных состояний.
  4. Проверить реализацию на задачах с известным ответом.
  5. Оценить численную ошибку Монте-Карло.
  6. Проверить устойчивость к увеличению статистики.
  7. Сравнить моделирование с контрольными данными.
  8. Документировать веса, отборы и систематические вариации.

Когда нужен другой инструмент

  • Для простого аналитического результата Монте-Карло может скрыть понимание.
  • Для крайне редких событий наивная генерация неэффективна.
  • Для сложного апостериорного распределения могут потребоваться Монте-Карло по марковским цепям или вложенная выборка.
  • Для обратных задач нужны методы регуляризации и развёртки.
  • Для дорогостоящего моделирования могут понадобиться суррогатные модели.

Что важно запомнить

Резюме

  • Монте-Карло строит численный ансамбль возможных реализаций модели.
  • Равномерные псевдослучайные числа преобразуются в требуемые распределения.
  • Интеграл можно представить как математическое ожидание.
  • Статистическая ошибка обычного Монте-Карло убывает как \(1/\sqrt N\).
  • Псевдоэксперименты проверяют всю статистическую процедуру.
  • Численное распространение ошибок сохраняет нелинейность, корреляции и асимметрию.
  • Моделирование имеет собственную конечную статистику и систематические неопределённости.
  • Результату Монте-Карло можно доверять только после проверок сходимости и воспроизводимости.