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

Аннотация

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

8.1 Физическая задача: слабый сигнал на фоне

Эксперимент ищет редкий процесс. За время набора данных ожидается в среднем $s=5$ сигнальных и $b=10$ фоновых событий. Полное число зарегистрированных событий имеет распределение Пуассона:

\[ N\sim\operatorname{Pois}(s+b). \tag{8.1}\]

Будем считать фон известным точно и оценим число сигнальных событий простым вычитанием:

\[ \widehat s=N-b. \tag{8.2}\]

До проведения эксперимента можно задать несколько вполне физических вопросов. Какое значение $\widehat s$ мы в среднем получим? Каков разброс оценки? Как часто флуктуация фона приведёт к отрицательному $\widehat s$?

Первые два ответа следуют из свойств распределения Пуассона:

\[ \mathbb E[\widehat s]=s, \qquad \operatorname{Var}(\widehat s)=s+b. \tag{8.3}\]

Для $s=5$ и $b=10$ стандартное отклонение равно $\sqrt{15}\simeq3.87$ события. Отрицательная оценка получается при $N<10$. Её вероятность равна

\[ P(\widehat s<0) =P(N\leq9\mid s+b=15) =\sum_{n=0}^{9}e^{-15}\frac{15^n}{n!} =0.0699. \tag{8.4}\]

Эту задачу можно решить и численно. Сгенерируем много значений $N_k\sim\operatorname{Pois}(15)$, для каждого вычислим $\widehat s_k=N_k-10$ и посмотрим на полученное распределение. Его среднее, дисперсия и доля отрицательных оценок приблизятся к значениям в формулах 8.3 и 8.4.

В этой игрушечной модели аналитический ответ помещается в несколько строк. Настоящий анализ использует энергетический спектр, отклик детектора, отбор событий и мешающие параметры. Распределение результата такой процедуры обычно уже нельзя выписать одной формулой. Тогда мы повторяем эксперимент на компьютере. Это и есть метод Монте-Карло.

8.2 Что вычисляет метод Монте-Карло

Пусть случайная величина $X$ имеет плотность $f_X(x)$, а интересующая нас величина определяется функцией $g$:

\[ X\sim f_X(x), \qquad Y=g(X). \tag{8.5}\]

Численный алгоритм состоит из трёх действий:

сгенерировать независимые значения $X_1,\ldots,X_N$ из распределения $f_X$;
вычислить $Y_k=g(X_k)$ для каждого события;
по выборке $Y_1,\ldots,Y_N$ найти среднее, дисперсию, вероятность нужного события или всё распределение $Y$.

В короткой записи

\[ X_k\sim f_X, \qquad Y_k=g(X_k), \qquad k=1,\ldots,N. \tag{8.6}\]

Один и тот же принцип используется в нескольких задачах. Интеграл можно записать как математическое ожидание,

\[ \int h(x)f_X(x)\,dx=\mathbb E[h(X)]. \tag{8.7}\]

При распространении ошибок генерируется случайный вектор $\mathbf X$, после чего для каждого события вычисляется $Y=g(\mathbf X)$. Псевдоэксперимент воспроизводит весь анализ целиком.

Метод был сформулирован как общий численный подход Метрополисом и Уламом в 1949 году [1]. Его точность определяется числом реализаций и дисперсией вычисляемой величины.

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

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

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

\[ U_1,U_2,\ldots, \qquad U_k\in[0,1), \tag{8.8}\]

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

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

Для параллельной генерации нужны независимые последовательности. Если тысячу процессов запустить с одним начальным состоянием, каждый из них создаст одну и ту же тысячу событий. В файле окажется миллион строк, а независимых событий останется тысяча. Поэтому вместе с конфигурацией анализа сохраняют алгоритм генератора и правило формирования начальных состояний.

Большинство генераторов непосредственно создаёт

\[ U\sim\mathrm U(0,1). \tag{8.9}\]

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

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

Пусть $F_X(x)$ — функция распределения требуемой величины:

\[ F_X(x)=P(X\leq x) =\int_{-\infty}^{x}f_X(t)\,dt. \tag{8.10}\]

Она монотонно возрастает от нуля до единицы. Если обратную функцию можно вычислить, берём $U\sim\mathrm U(0,1)$ и полагаем

\[ \boxed{X=F_X^{-1}(U).} \tag{8.11}\]

8.4.0.1 Почему это работает

Монотонность $F_X$ даёт

\[ \begin{aligned} P(X\leq x) &=P\!\left(F_X^{-1}(U)\leq x\right)\\ &=P\!\left(U\leq F_X(x)\right)\\ &=F_X(x). \end{aligned} \tag{8.12}\]

В последней строке использовано $P(U\leq u)=u$ для равномерной величины на $[0,1]$. Получена требуемая функция распределения. $\blacksquare$

8.4.1 Пример: время между событиями

Пусть события образуют пуассоновский процесс с интенсивностью $\lambda$. Время $X$ до следующего события имеет экспоненциальную плотность

\[ f_X(x)=\lambda e^{-\lambda x}, \qquad x\geq0, \tag{8.13}\]

и функцию распределения

\[ F_X(x)=1-e^{-\lambda x}. \tag{8.14}\]

Решая уравнение $U=F_X(X)$ относительно $X$, получаем

\[ X=-\frac{1}{\lambda}\ln(1-U). \tag{8.15}\]

Величина $1-U$ также равномерна на $[0,1]$, поэтому в программе часто пишут

\[ \boxed{X=-\frac{1}{\lambda}\ln U.} \tag{8.16}\]

8.5 Метод отбора

Обратную функцию удаётся найти не для каждого распределения. Пусть плотность $f(x)$ ограничена сверху другой плотностью $g(x)$:

\[ f(x)\leq M g(x), \tag{8.17}\]

причём из $g(x)$ генерировать удобно. Алгоритм отбора устроен так:

сгенерировать кандидата $X\sim g(x)$;
независимо сгенерировать $U\sim\mathrm U(0,1)$;
принять $X$, если \[ U\leq\frac{f(X)}{M g(X)}; \tag{8.18}\]
при невыполнении условия повторить попытку.

Если обе плотности нормированы, средняя доля принятых событий равна $1/M$. Чем ближе оболочка $Mg(x)$ к $f(x)$, тем меньше вычислений заканчивается отказом.

8.5.1 Пример: плотность $6x(1-x)$

Рассмотрим плотность

\[ f(x)=6x(1-x), \qquad 0\leq x\leq1. \tag{8.19}\]

Её максимум достигается при $x=1/2$ и равен $3/2$. Равномерно генерируем точки в прямоугольнике $0\leq x\leq1$, $0\leq y\leq3/2$ и принимаем точку, если $y\leq f(x)$.

viewof rejectN = Inputs.range([100, 6000], {value: 1200, step: 100, label: "N точек"})
viewof regenerateReject = Inputs.button("Новая выборка", {
  value: 0,
  reduce: value => value + 1
})

seededUniformReject = seed => {
  let state = seed >>> 0;
  return () => {
    state = (state + 0x6D2B79F5) >>> 0;
    let value = state;
    value = Math.imul(value ^ (value >>> 15), value | 1);
    value ^= value + Math.imul(value ^ (value >>> 7), value | 61);
    return ((value ^ (value >>> 14)) >>> 0) / 4294967296;
  };
}

rejectSeed = (20260908 + 1009 * rejectN + 104729 * regenerateReject) >>> 0

rejectPoints = {
  const random = seededUniformReject(rejectSeed);
  const values = [];
  for (let i = 0; i < rejectN; ++i) {
    const x = random();
    const y = 1.5 * random();
    values.push({x, y, accepted: y <= 6 * x * (1 - x)});
  }
  return values;
}

rejectCurve = Array.from({length: 301}, (_, i) => {
  const x = i / 300;
  return {x, y: 6 * x * (1 - x)};
})

acceptedReject = rejectPoints.filter(d => d.accepted)
rejectRate = acceptedReject.length / rejectPoints.length

rejectPlot = Plot.plot({
  width: 820,
  height: 430,
  marginLeft: 65,
  x: {label: "x", domain: [0, 1], grid: true},
  y: {label: "y", domain: [0, 1.55], grid: true},
  marks: [
    Plot.dot(rejectPoints.filter(d => !d.accepted), {x: "x", y: "y", r: 2.2, fill: "#777", fillOpacity: 0.35}),
    Plot.dot(acceptedReject, {x: "x", y: "y", r: 2.5, fill: "#4ea1ff", fillOpacity: 0.65}),
    Plot.line(rejectCurve, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 4})
  ]
})

html`
<div class="book-applet-panel">
  ${rejectPlot}
  <div class="book-applet-summary">
    <span>Принято: ${acceptedReject.length} из ${rejectPoints.length};</span>
    <span>эффективность: ${(100 * rejectRate).toFixed(1)}%;</span>
    <span>теоретическое значение: 66.7%.</span>
  </div>
</div>`

Метод отбора для плотности $6x(1-x)$ — Метод отбора для плотности 6x(1 − x)

8.5.2 Как выбирать способ генерации

Таблица 8.1: Основные способы генерации случайных величин

Метод	Сильная сторона	Ограничение
Обратная функция распределения	используется каждое сгенерированное число	нужна вычислимая $F_X^{-1}$
Метод отбора	применим к широкому классу плотностей	плохая оболочка снижает эффективность
Специальный алгоритм	может быть очень быстрым	нужна отдельная реализация
Марковская цепь	работает со сложными многомерными плотностями	значения коррелированы, сходимость требует проверки

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

Рассмотрим одномерный интеграл

\[ I=\int_a^b h(x)\,dx. \tag{8.20}\]

Если $U\sim\mathrm U(a,b)$, то

\[ \mathbb E[h(U)] =\frac{1}{b-a}\int_a^b h(x)\,dx. \tag{8.21}\]

Математическое ожидание заменим выборочным средним:

\[ \boxed{ \widehat I_N =\frac{b-a}{N}\sum_{k=1}^{N}h(U_k). } \tag{8.22}\]

8.6.1 Статистическая ошибка интеграла

8.6.1.1 Дисперсия оценки $\widehat I_N$

Обозначим

\[ \operatorname{Var}[h(U)]=\sigma_h^2. \tag{8.23}\]

Для независимых точек дисперсии слагаемых складываются:

\[ \begin{aligned} \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}. \end{aligned} \tag{8.24}\]

Следовательно,

\[ \boxed{ \sigma_{\widehat I_N} =\frac{(b-a)\sigma_h}{\sqrt N}. } \tag{8.25}\]

$\blacksquare$

Чтобы уменьшить ошибку в десять раз, число точек придётся увеличить в сто раз. Это медленная сходимость. Её важное достоинство проявляется в большой размерности: показатель $N^{-1/2}$ не превращается в $N^{-1/d}$. Коэффициент перед ним, конечно, зависит от конкретной функции и выбранного распределения точек.

8.7 Проверка на числе $\pi$

Равномерно сгенерируем точки в единичном квадрате:

\[ (X,Y)\sim\mathrm U([0,1]^2). \tag{8.26}\]

Точка лежит внутри четверти единичного круга при $X^2+Y^2\leq1$. Поэтому

\[ P(X^2+Y^2\leq1)=\frac{\pi}{4}, \tag{8.27}\]

и доля попавших внутрь точек даёт оценку

\[ \boxed{ \widehat\pi=4\frac{N_{\mathrm{in}}}{N}. } \tag{8.28}\]

viewof piN = Inputs.range([100, 12000], {value: 1200, step: 100, label: "число точек N"})
viewof regeneratePi = Inputs.button("Новая выборка", {
  value: 0,
  reduce: value => value + 1
})

seededUniformPi = seed => {
  let state = seed >>> 0;
  return () => {
    state = (state + 0x6D2B79F5) >>> 0;
    let value = state;
    value = Math.imul(value ^ (value >>> 15), value | 1);
    value ^= value + Math.imul(value ^ (value >>> 7), value | 61);
    return ((value ^ (value >>> 14)) >>> 0) / 4294967296;
  };
}

piSeed = (20260908 + 1009 * piN + 104729 * regeneratePi) >>> 0

piPoints = {
  const random = seededUniformPi(piSeed);
  const values = [];
  for (let i = 0; i < piN; ++i) {
    const x = random();
    const y = random();
    values.push({x, y, inside: x * x + y * y <= 1});
  }
  return values;
}

piInside = piPoints.filter(d => d.inside)
piEstimate = 4 * piInside.length / piPoints.length
piSe = 4 * Math.sqrt((piInside.length / piPoints.length) * (1 - piInside.length / piPoints.length) / piPoints.length)
quarterCircle = Array.from({length: 301}, (_, i) => {
  const x = i / 300;
  return {x, y: Math.sqrt(1 - x * x)};
})

piPlot = Plot.plot({
  width: 620,
  height: 560,
  marginLeft: 55,
  x: {label: "x", domain: [0, 1], grid: true},
  y: {label: "y", domain: [0, 1], grid: true},
  marks: [
    Plot.dot(piPoints.filter(d => !d.inside), {x: "x", y: "y", r: 2, fill: "#ff6b6b", fillOpacity: 0.45}),
    Plot.dot(piInside, {x: "x", y: "y", r: 2, fill: "#4ea1ff", fillOpacity: 0.50}),
    Plot.line(quarterCircle, {x: "x", y: "y", stroke: "#ffb347", strokeWidth: 4})
  ]
})

html`
<div class="book-applet-panel">
  ${piPlot}
  <div class="book-applet-summary">
    <span>Оценка π: ${piEstimate.toFixed(5)};</span>
    <span>статистическая ошибка: ${piSe.toFixed(5)};</span>
    <span>отклонение от π: ${(piEstimate - Math.PI).toFixed(5)}.</span>
  </div>
</div>`

Оценка числа $\pi$ методом Монте-Карло — Оценка числа π методом Монте-Карло

Если единственная цель состоит в вычислении числа $\pi$, калькулятор справится лучше. Этот пример ценен другим: точный ответ известен, геометрия видна на рисунке, а ошибка подчиняется закону $1/\sqrt N$ из формулы 8.25. Такие задачи нужны для проверки реализации метода.

8.8 Высокая размерность

Регулярная сетка с $m$ узлами по каждой координате в $d$ измерениях содержит

\[ N_{\mathrm{grid}}=m^d \tag{8.29}\]

точек. Сто узлов по каждой из шести координат дают $10^{12}$ вычислений функции. Сетка ещё довольно грубая, а вычислительная задача уже нет.

function mcRange(label, domain, options) {
  const [min, max] = domain;
  const value = options.value;
  const step = options.step ?? 1;
  const control = html`<label class="mc-range-control">
    <span class="mc-range-label">${label}</span>
    <input class="mc-range-number" type="number" min=${min} max=${max} step=${step} value=${value}>
    <input class="mc-range-slider" type="range" min=${min} max=${max} step=${step} value=${value}>
  </label>`;
  const number = control.querySelector(".mc-range-number");
  const slider = control.querySelector(".mc-range-slider");
  control.value = Number(value);
  const sync = event => {
    const next = Number(event.target.value);
    control.value = next;
    number.value = next;
    slider.value = next;
    control.dispatchEvent(new CustomEvent("input", {bubbles: true}));
  };
  number.addEventListener("input", sync);
  slider.addEventListener("input", sync);
  return control;
}

viewof grid_m = mcRange("точек по каждой координате", [2, 12], {value: 5, step: 1})
viewof grid_d = mcRange("размерность", [1, 12], {value: 8, step: 1})

gridPoints = Array.from({length: 12}, (_, index) => {
  const dimension = index + 1;
  return {
    dimension,
    points: Math.pow(grid_m, dimension)
  };
})

gridPlot = Plot.plot({
  width: 820,
  height: 390,
  marginLeft: 82,
  x: {
    label: "размерность d",
    domain: [1, 12],
    ticks: 12,
    grid: true
  },
  y: {
    label: "число узлов mᵈ",
    type: "log",
    grid: true
  },
  marks: [
    Plot.line(gridPoints, {
      x: "dimension",
      y: "points",
      stroke: "#216b73",
      strokeWidth: 3
    }),
    Plot.dot(gridPoints.filter(point => point.dimension === grid_d), {
      x: "dimension",
      y: "points",
      fill: "#ff8c42",
      r: 6
    })
  ]
})

gridResult = {
  const points = Math.pow(grid_m, grid_d);
  return html`<div class="book-applet-panel">
    ${gridPlot}
    <div class="book-applet-summary">
      <span>${grid_m} узлов по каждой координате;</span>
      <span>размерность ${grid_d};</span>
      <span>всего ${points.toLocaleString("ru-RU")} вычислений.</span>
    </div>
  </div>`
}

Рисунок 8.1: При двенадцати узлах по каждой из восьми координат число вычислений возрастает до $12^8=429\,981\,696$

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

8.9 Выборка по важности

Пусть требуется вычислить

\[ I=\int h(x)f(x)\,dx. \tag{8.30}\]

Равномерная генерация не учитывает форму интегранда. Если основной вклад сосредоточен в малой области, большая часть вычислений почти ничего не добавит к результату. Выберем плотность $q(x)$, положительную везде, где $h(x)f(x)\neq0$, и перепишем интеграл:

\[ \boxed{ I=\mathbb E_q\!\left[ h(X)\frac{f(X)}{q(X)} \right]. } \tag{8.31}\]

Точка $X_k\sim q(x)$ получает вес

\[ w_k=\frac{f(X_k)}{q(X_k)}, \tag{8.32}\]

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

\[ \widehat I_N =\frac1N\sum_{k=1}^{N}w_k h(X_k). \tag{8.33}\]

Хорошая плотность $q(x)$ часто приводит нас в область большого вклада и даёт веса одного порядка. Если редкое событие получает огромный вес, оно может определить весь ответ и его дисперсию. Число строк в файле тогда плохо характеризует статистическую точность. Для положительных весов полезно вычислять эффективный размер выборки

\[ N_{\mathrm{eff}} =\frac{\left(\sum_k w_k\right)^2}{\sum_k w_k^2}. \tag{8.34}\]

При одинаковых весах $N_{\mathrm{eff}}=N$. Сильно различающиеся веса дают $N_{\mathrm{eff}}\ll N$.

8.9.1 Редкая область

Следующий интерактив показывает одну и ту же малую область при двух способах генерации. Равномерная выборка покрывает весь квадрат. Выборка по важности сосредоточивает точки вблизи области, которая даёт основной вклад. При численном интегрировании изменение плотности обязательно компенсируется весом из формулы 8.32.

viewof importanceN = Inputs.range([200, 6000], {
  value: 1500,
  step: 100,
  label: "число точек"
})

viewof importanceSigma = Inputs.range([0.02, 0.18], {
  value: 0.06,
  step: 0.005,
  label: "ширина важной области"
})

function rng32(seed) {
  return function() {
    seed |= 0; seed = seed + 0x6D2B79F5 | 0;
    let t = Math.imul(seed ^ seed >>> 15, 1 | seed);
    t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
    return ((t ^ t >>> 14) >>> 0) / 4294967296;
  }
}

rare_canvas = {
  const width = 720, height = 410;
  const canvas = DOM.canvas(width, height);
  const ctx = canvas.getContext("2d");
  const rnd = rng32(20260618);
  const L = 320, x0 = 50, y0 = 48;
  const cx = x0 + 0.68 * L, cy = y0 + 0.42 * L, r = importanceSigma * L;
  ctx.clearRect(0, 0, width, height);
  ctx.strokeStyle = "#666"; ctx.lineWidth = 2; ctx.strokeRect(x0, y0, L, L);
  ctx.beginPath(); ctx.arc(cx, cy, r, 0, 2 * Math.PI); ctx.strokeStyle = "#2163ff"; ctx.stroke();
  let hit = 0;
  for (let i = 0; i < importanceN; ++i) {
    const x = rnd(), y = rnd();
    const px = x0 + x * L, py = y0 + y * L;
    const inside = (px-cx)**2 + (py-cy)**2 < r**2;
    if (inside) hit += 1;
    ctx.beginPath(); ctx.arc(px, py, inside ? 2.5 : 1.4, 0, 2*Math.PI);
    ctx.fillStyle = inside ? "#2163ff" : "rgba(100,100,100,0.28)"; ctx.fill();
  }
  ctx.fillStyle = "#222"; ctx.font = "20px system-ui";
  ctx.fillText(`равномерно: попало ${hit} из ${importanceN}`, 395, 180);
  return canvas;
}

importance_canvas = {
  const width = 720, height = 410;
  const canvas = DOM.canvas(width, height);
  const ctx = canvas.getContext("2d");
  const rnd = rng32(424242);
  function normal() {
    const u1 = Math.max(rnd(), 1e-12), u2 = rnd();
    return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
  }
  const L = 320, x0 = 50, y0 = 48;
  const cx = x0 + 0.68 * L, cy = y0 + 0.42 * L, r = importanceSigma * L;
  ctx.clearRect(0, 0, width, height);
  ctx.strokeStyle = "#666"; ctx.lineWidth = 2; ctx.strokeRect(x0, y0, L, L);
  ctx.beginPath(); ctx.arc(cx, cy, r, 0, 2 * Math.PI); ctx.strokeStyle = "#2163ff"; ctx.stroke();
  let hit = 0;
  for (let i = 0; i < importanceN; ++i) {
    const x = 0.68 + importanceSigma * normal();
    const y = 0.42 + importanceSigma * normal();
    if (x < 0 || x > 1 || y < 0 || y > 1) continue;
    const px = x0 + x * L, py = y0 + y * L;
    const inside = (px-cx)**2 + (py-cy)**2 < r**2;
    if (inside) hit += 1;
    ctx.beginPath(); ctx.arc(px, py, inside ? 2.5 : 1.5, 0, 2*Math.PI);
    ctx.fillStyle = inside ? "#2163ff" : "rgba(100,100,100,0.28)"; ctx.fill();
  }
  ctx.fillStyle = "#222"; ctx.font = "20px system-ui";
  ctx.fillText(`по важности: в области ${hit} точек`, 395, 180);
  return canvas;
}

Равномерная выборка и выборка по важности

Рисунок 8.2: Выборка по важности сосредоточивает точки вблизи области, которая даёт основной вклад в интеграл

Рисунок 8.3: Умеренные веса дают много сопоставимых вкладов; редкие события с очень большими весами увеличивают дисперсию оценки

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

Для ручного выбора $q(\mathbf x)$ надо заранее знать важную область. Алгоритм VEGAS строит приближение к такой плотности по первым итерациям расчёта [2]. В каждой координате он сгущает сетку там, где интегранд вносит большой вклад. Полученная плотность имеет приближённо разделяющийся вид

\[ q(\mathbf x)=\prod_{i=1}^{d}q_i(x_i). \tag{8.35}\]

Рисунок 8.4: После нескольких итераций VEGAS располагает больше интервалов в области большого вклада

VEGAS сохраняет асимптотический закон $1/\sqrt N$. Выигрыш получается за счёт меньшей дисперсии. Разделяющаяся плотность хорошо описывает пики, положение которых видно по отдельным координатам. Узкий диагональный гребень или несколько далеко разнесённых пиков представляют более трудную задачу.

8.10.1 Контрольный восьмимерный интеграл

Рассмотрим произведение нормированных гауссовых плотностей в единичном гиперкубе:

\[ f(\mathbf x)= \prod_{i=1}^{8} \frac{1}{\sqrt{2\pi}\sigma} \exp\!\left[-\frac{(x_i-\mu_i)^2}{2\sigma^2}\right], \qquad 0\leq x_i\leq1, \tag{8.36}\]

где $\sigma=0.10$, а центры равны

\[ \boldsymbol\mu= (0.63,0.37,0.58,0.42,0.69,0.31,0.54,0.46). \tag{8.37}\]

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

\[ I_{\mathrm{exact}}=0.9978196359. \tag{8.38}\]

Простой Монте-Карло с одним и тем же начальным числом даёт результаты из таблицы 8.2.

Таблица 8.2: Равномерный метод Монте-Карло для восьмимерного интеграла

$N$	$\widehat I_N$
$10^4$	$0.946\pm0.48$
$10^5$	$0.864\pm0.13$
$10^6$	$1.001\pm0.061$

Миллион точек всё ещё даёт ошибку около шести процентов: почти вся равномерная выборка проходит мимо узкого пика. Тензорная квадратура Гаусса—Лежандра сходится, но число вызовов функции растёт как $m^8$.

Таблица 8.3: Тензорная квадратура для восьмимерного интеграла

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

Числа в обеих таблицах получены скриптом shared/notebooks/mc_vegas_demo.py с начальным состоянием 12345. После восьми адаптационных итераций по 20 000 точек и десяти рабочих итераций по 50 000 точек VEGAS даёт

\[ \widehat I_{\mathrm{VEGAS}}=0.99787\pm0.00053. \tag{8.39}\]

Этот пример специально удобен для VEGAS: интегранд раскладывается в произведение функций отдельных координат. Для нового интеграла результат проверяют независимыми запусками, увеличением числа точек и сравнением с более грубым методом. Очень узкий пик алгоритм может не встретить на первых итерациях, а значит, и не научиться генерировать точки в его окрестности.

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

Псевдоэксперимент — одна искусственная реализация измерения при заданной модели. В простом случае достаточно сгенерировать одно число. В полном анализе последовательно выполняются генерация физического процесса, моделирование детектора, реконструкция, отбор и та же статистическая процедура, которая применяется к данным.

Ансамбль псевдоэкспериментов позволяет найти:

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

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

8.11.1 Возвращаемся к счётному эксперименту

Для модели из формулы 8.1 один псевдоэксперимент состоит из генерации $N\sim\operatorname{Pois}(s+b)$ и вычисления $\widehat s=N-b$. Повторив эти действия несколько тысяч раз, получим распределение оценки.

viewof toySignal = Inputs.range([0, 20], {
  value: 5,
  step: 0.5,
  label: "ожидаемый сигнал s"
})

viewof toyBackground = Inputs.range([0.5, 30], {
  value: 10,
  step: 0.5,
  label: "ожидаемый фон b"
})

viewof toyN = Inputs.range([500, 20000], {
  value: 6000,
  step: 500,
  label: "число псевдоэкспериментов"
})

viewof regenerateToys = Inputs.button("Новый ансамбль", {
  value: 0,
  reduce: value => value + 1
})

seededUniformToys = seed => {
  let state = seed >>> 0;
  return () => {
    state = (state + 0x6D2B79F5) >>> 0;
    let value = state;
    value = Math.imul(value ^ (value >>> 15), value | 1);
    value ^= value + Math.imul(value ^ (value >>> 7), value | 61);
    return ((value ^ (value >>> 14)) >>> 0) / 4294967296;
  };
}

poissonToyCh8 = (lambda, random) => {
  const limit = Math.exp(-lambda);
  let product = 1;
  let count = 0;
  do {
    count += 1;
    product *= random();
  } while (product > limit);
  return count - 1;
}

toySeed = (
  20260908
  + 101 * Math.round(2 * toySignal)
  + 1009 * Math.round(2 * toyBackground)
  + 17 * toyN
  + 104729 * regenerateToys
) >>> 0

toyEstimates = {
  const random = seededUniformToys(toySeed);
  const lambda = toySignal + toyBackground;
  return Array.from(
    {length: toyN},
    () => poissonToyCh8(lambda, random) - toyBackground
  );
}

toyBins = {
  const counts = new Map();
  for (const value of toyEstimates) {
    counts.set(value, (counts.get(value) || 0) + 1);
  }
  return Array.from(counts, ([value, count]) => ({
    value,
    x1: value - 0.45,
    x2: value + 0.45,
    probability: count / toyEstimates.length
  })).sort((a, b) => a.value - b.value);
}

toyExact = {
  const lambda = toySignal + toyBackground;
  const maxCount = Math.max(
    Math.ceil(lambda + 8 * Math.sqrt(lambda)),
    Math.ceil(Math.max(...toyEstimates) + toyBackground)
  );
  const values = [];
  let probability = Math.exp(-lambda);
  for (let count = 0; count <= maxCount; ++count) {
    if (count > 0) probability *= lambda / count;
    values.push({
      value: count - toyBackground,
      probability
    });
  }
  return values;
}

toyMean = toyEstimates.reduce((sum, value) => sum + value, 0) / toyEstimates.length
toySigma = Math.sqrt(
  toyEstimates.reduce((sum, value) => sum + (value - toyMean) ** 2, 0)
  / (toyEstimates.length - 1)
)
toyNegative = toyEstimates.filter(value => value < 0).length / toyEstimates.length
toyNegativeExact = toyExact
  .filter(point => point.value < 0)
  .reduce((sum, point) => sum + point.probability, 0)

toyMin = Math.min(
  Math.floor(Math.min(...toyEstimates)) - 0.5,
  -toyBackground - 0.5
)
toyMax = Math.max(
  Math.ceil(Math.max(...toyEstimates)) + 0.5,
  toySignal + 7 * Math.sqrt(toySignal + toyBackground)
)
toyYMax = 1.18 * Math.max(
  ...toyBins.map(point => point.probability),
  ...toyExact.map(point => point.probability)
)

toyPlot = Plot.plot({
  width: 900,
  height: 390,
  marginLeft: 72,
  x: {
    label: "оценка сигнала ŝ = N − b",
    domain: [toyMin, toyMax],
    grid: true
  },
  y: {
    label: "вероятность",
    domain: [0, toyYMax],
    grid: true
  },
  marks: [
    Plot.rectY(toyBins, {
      x1: "x1",
      x2: "x2",
      y: "probability",
      fill: "#4ea1ff",
      fillOpacity: 0.72
    }),
    Plot.line(toyExact, {
      x: "value",
      y: "probability",
      stroke: "#ff8c42",
      strokeWidth: 2.2
    }),
    Plot.dot(toyExact, {
      x: "value",
      y: "probability",
      fill: "#ff8c42",
      r: 3.3
    }),
    Plot.ruleX([toySignal], {stroke: "#216b73", strokeWidth: 3}),
    Plot.ruleX([0], {stroke: "#777", strokeDasharray: "4,4"})
  ]
})

html`
<div class="book-applet-panel">
  ${toyPlot}
  <div class="book-applet-summary">
    <span>Монте-Карло: среднее ${toyMean.toFixed(3)}, σ ${toySigma.toFixed(3)}, P(ŝ&lt;0) ${(100 * toyNegative).toFixed(2)}%;</span>
    <span>точно: среднее ${toySignal.toFixed(3)}, σ ${Math.sqrt(toySignal + toyBackground).toFixed(3)}, P(ŝ&lt;0) ${(100 * toyNegativeExact).toFixed(2)}%;</span>
    <span>синие столбцы — ансамбль, оранжевые точки — распределение Пуассона.</span>
  </div>
</div>`

Отрицательное $\widehat s$ не означает отрицательного числа физических событий. Так проявляется флуктуация оценки $N-b$. Статистическая процедура должна уметь работать и с такими реализациями; удаление их из ансамбля изменит среднее и нарушит свойства метода.

Реальный эксперимент даёт одну точку из показанного распределения. Увеличение числа псевдоэкспериментов точнее определяет это распределение, но не уменьшает разброс результата одного реального измерения.

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

В физике частиц и нейтрино один псевдоэксперимент может опираться на длинную цепочку моделирования:

генератор задаёт процесс, типы частиц и их кинематику;
моделируется прохождение частиц через вещество, потери энергии и вторичные частицы;
вычисляется отклик сенсоров и электроники;
программа реконструкции восстанавливает энергию, координаты и тип события;
применяются триггер, критерии качества и отбор анализа.

Каждый этап содержит приближения. Их проверяют по контрольным данным, а оставшиеся расхождения входят в систематическую неопределённость.

В ролике Саймона Блита показана работа Opticks: распространение оптических фотонов рассчитывается на графическом процессоре для геометрий детекторов Daya Bay и JUNO. Описание проекта опубликовано на странице автора [3].

Моделирование оптических фотонов в Daya Bay и JUNO

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

Пусть истинный спектр разбит на бины. Элемент матрицы отклика

\[ R_{ji} =P(\text{событие восстановлено в бине }j \mid\text{истинный бин }i) \tag{8.40}\]

оценивается по моделированию. Если истинный спектр равен $\mathbf t$, то ожидаемый восстановленный спектр имеет вид

\[ \boxed{\mathbf r=R\mathbf t.} \tag{8.41}\]

Диагональные элементы описывают правильную реконструкцию бина, внедиагональные — миграции. Сумма элементов столбца может быть меньше единицы: часть событий теряется из-за геометрического аксептанса, порогов и отбора.

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

Если из $N_{\mathrm{gen}}$ сгенерированных событий прошло отбор $N_{\mathrm{sel}}$, оценка эффективности равна

\[ \widehat\varepsilon =\frac{N_{\mathrm{sel}}}{N_{\mathrm{gen}}}. \tag{8.42}\]

Для независимых событий с единичными весами

\[ N_{\mathrm{sel}} \sim\operatorname{Bin}(N_{\mathrm{gen}},\varepsilon), \tag{8.43}\]

поэтому стандартную ошибку обычно оценивают как

\[ \boxed{ \sigma_{\widehat\varepsilon} \simeq \sqrt{\frac{ \widehat\varepsilon(1-\widehat\varepsilon) }{N_{\mathrm{gen}}}}. } \tag{8.44}\]

У этой формулы хорошо видны границы применимости. Если все события прошли отбор, она даёт $\widehat\varepsilon=1$ и нулевую ошибку. Из конечной выборки следует лишь, что в ней не встретился ни один отказ. Следующее событие пройти отбор не обязано. Вблизи физических границ $0\leq\varepsilon\leq1$ нужны интервальные методы; они рассматриваются в следующих главах.

При наличии весов биномиальная модель уже неприменима в таком виде. Тогда проверяют суммы весов, суммы квадратов весов и эффективный размер выборки из формулы 8.34.

8.13 Как проверять расчёт

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

Если из псевдоэксперимента восстанавливается параметр $\theta$, на ансамбле с известным истинным значением вычисляют смещение

\[ \operatorname{bias}(\widehat\theta) =\mathbb E[\widehat\theta]-\theta \tag{8.45}\]

и дисперсию

\[ \operatorname{Var}(\widehat\theta) =\mathbb E\!\left[ (\widehat\theta-\mathbb E[\widehat\theta])^2 \right]. \tag{8.46}\]

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

Здесь участвуют два разных разброса. Физическая статистическая неопределённость описывает возможные результаты одного реального эксперимента. Численная ошибка Монте-Карло показывает, насколько точно конечная симуляция оценивает это распределение. Большее число псевдоэкспериментов уменьшает вторую ошибку. Первая определяется объёмом реальных данных.

Для взвешенной выборки отдельно проверяют распределение весов и $N_{\mathrm{eff}}$. Для адаптивного интегрирования сравнивают итерации и независимые начальные состояния. Для полного моделирования тестируют каждый этап цепочки и сопоставляют контрольные распределения с данными.

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

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

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

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

8.15 Итоги главы

Метод Монте-Карло строит численный ансамбль реализаций заданной модели.
Равномерные псевдослучайные числа преобразуются в нужное распределение обратной функцией, методом отбора или специальным алгоритмом.
Ошибка простого интегрирования убывает как $1/\sqrt N$; в большой размерности это часто выгоднее регулярной сетки.
Выборка по важности и VEGAS уменьшают дисперсию, направляя вычисления в области большого вклада.
Псевдоэксперименты позволяют проверить всю статистическую процедуру от генерации данных до оценки параметра.
Статистическая неопределённость реального эксперимента и численная ошибка конечной симуляции являются разными величинами.

8.16 Задачи

Задача 1. Обратная функция распределения

Пусть \[ U\sim\mathrm U(0,1). \]

Покажите, что \[ X=-\frac{1}{\lambda}\ln(1-U) \] имеет экспоненциальное распределение с параметром $\lambda$.
Сгенерируйте выборку $X$ и сравните гистограмму с аналитической плотностью.
Проверьте выборочные среднее и дисперсию.

Задача 2. Метод отбора

Требуется генерировать плотность \[ f(x)=6x(1-x), \qquad 0\le x\le1. \]

Найдите максимум $f(x)$.
Постройте алгоритм отбора из равномерного распределения.
Вычислите теоретическую эффективность алгоритма.
Проверьте её численно.

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

Вычислите \[ I=\int_0^1\frac{4}{1+x^2}\,dx=\pi \] методом Монте-Карло.

Запишите оценку $\widehat I_N$.
Оцените её статистическую ошибку по выборке.
Исследуйте зависимость ошибки от $N$.
Покажите численно, что ошибка убывает приблизительно как $N^{-1/2}$.

Задача 4. Редкое событие и выборка по важности

Требуется оценить вероятность \[ p=P(Z>4), \qquad Z\sim\mathrm N(0,1). \]

Оцените $p$ наивной генерацией из стандартного гаусса.
Объясните, почему для устойчивой оценки требуется очень большая выборка.
Предложите распределение $q(z)$, чаще генерирующее область $z>4$.
Запишите соответствующий вес $w(z)=f(z)/q(z)$.
Сравните дисперсии двух оценок.

Задача 5. Псевдоэксперименты для счётного эксперимента

Пусть \[ N\sim\operatorname{Pois}(s+b), \qquad \widehat s=N-b, \] где фон $b$ известен точно.

Найдите аналитические $\mathbb E[\widehat s]$ и $\operatorname{Var}(\widehat s)$.
Постройте ансамбль псевдоэкспериментов.
Сравните выборочные среднее и дисперсию с аналитическими значениями.
Оцените вероятность получить $\widehat s<0$.

Задача 6. Нелинейное распространение ошибки

Пусть коррелированные величины имеют распределение \[ \begin{pmatrix}X\\Y\end{pmatrix} \sim \mathrm N\!\left[ \begin{pmatrix}10\\5\end{pmatrix}, \begin{pmatrix} \sigma_X^2 & \rho\sigma_X\sigma_Y\\ \rho\sigma_X\sigma_Y & \sigma_Y^2 \end{pmatrix} \right]. \]

Для отношения \[ R=\frac{X}{Y} \] сравните:

линейную оценку ошибки;
стандартное отклонение выборки Монте-Карло;
центральный квантильный интервал $[q_{0.16},q_{0.84}]$.

Исследуйте зависимость результатов от $\sigma_Y$ и $\rho$.

Задача 7. Эффективность отбора

Из $N_{\mathrm{gen}}=10\,000$ сгенерированных событий отбор прошло $N_{\mathrm{sel}}=730$.

Оцените эффективность отбора.
Оцените её биномиальную статистическую ошибку.
Сколько событий требуется сгенерировать, чтобы относительная ошибка эффективности была меньше $1\%$?
Обсудите, как изменится задача при наличии весов событий.

Задача 8. Эффективный размер выборки

Для взвешенной выборки определено \[ N_{\mathrm{eff}} = \frac{(\sum_k w_k)^2}{\sum_k w_k^2}. \]

Покажите, что при одинаковых весах $N_{\mathrm{eff}}=N$.
Найдите $N_{\mathrm{eff}}$ для весов \[ (1,1,1,1,6). \]
Сравните результат с фактическим числом событий.
Объясните, почему большая вариация весов ухудшает статистическую точность.

Литература

[1] Metropolis, N. и Ulam, S. The Monte Carlo Method. Journal of the American Statistical Association, 44, 335–341. 1949. doi:10.1080/01621459.1949.10483310.

[2] Lepage, G. P. A New Algorithm for Adaptive Multidimensional Integration. Journal of Computational Physics, 27, 192–203. 1978. doi:10.1016/0021-9991(78)90004-9.

[3] Blyth, S. C. Opticks: GPU Optical Photon Simulation. https://juno.ihep.ac.cn/~blyth/.

\(N\)	\(\widehat I_N\)
\(10^4\)	\(0.946\pm0.48\)
\(10^5\)	\(0.864\pm0.13\)
\(10^6\)	\(1.001\pm0.061\)