Методи Монте-Карло — випадкова вибірка для симуляції та інтегрування

1. Основна ідея: вибірка → оцінка

Кожен метод Монте-Карло дотримується одного й того ж трикрокового шаблону:

1. Визначити випадкову величину X, чиє математичне сподівання E[X] дорівнює потрібній вам величині.
2. Зробити N незалежних вибірок x₁, x₂, …, x_N з розподілу X.
3. Оцінити величину як вибіркове середнє: μ̂ = (x₁ + x₂ + … + xN) / N.

За законом великих чисел μ̂ збігається до E[X] при N → ∞. За центральною граничною теоремою похибка приблизно гаусова зі стандартним відхиленням σ/√N, де σ² = Var(X).

2. Оцінка π випадковими дротиками

Класичний вступ: кинути N випадкових дротиків рівномірно в одиничний квадрат [0,1]². Підрахувати частку, що приземлилася всередині чверті кола радіуса 1 з центром у початку координат. Оскільки чверть кола має площу π/4, частка збігається до π/4.

function estimatePi(N) {
  let inside = 0;
  for (let i = 0; i < N; i++) {
    const x = Math.random();
    const y = Math.random();
    if (x * x + y * y <= 1) inside++;
  }
  return 4 * inside / N;
}

console.log(estimatePi(1_000_000));  // ≈ 3.1415…  (похибка ~ 0.001)

За N = 10⁶ типова похибка близько 0.002 — приблизно 3 значущі цифри. Щоб отримати ще одну значущу цифру, потрібно у 100× більше вибірок: N = 10⁸. Ця збіжність 1/√N характерна для всіх методів Монте-Карло.

3. Інтегрування методом Монте-Карло

Загальна мета — обчислити інтеграл:

Якщо ми можемо рівномірно вибирати точки з Ω об’ємом |Ω|, оцінка така:

// Інтегруємо f(x) = x² на [0, 1] — точна відповідь: 1/3
function mcIntegrate(f, a, b, N) {
  let sum = 0;
  for (let i = 0; i < N; i++) {
    const x = a + Math.random() * (b - a);
    sum += f(x);
  }
  return (b - a) * sum / N;
}

const result = mcIntegrate(x => x * x, 0, 1, 1e6);
console.log(result);   // ≈ 0.3333…

// Багатовимірний: об'єм 4D одиничної гіперсфери (точно: π²/2 ≈ 4.935)
function hypersphereMC(d, N) {
  let inside = 0;
  for (let i = 0; i < N; i++) {
    let r2 = 0;
    for (let j = 0; j < d; j++) { const x = 2*Math.random()-1; r2 += x*x; }
    if (r2 <= 1) inside++;
  }
  return Math.pow(2, d) * inside / N;  // об'єм охоплюючого гіперкуба = 2^d
}
console.log(hypersphereMC(4, 1e6));   // ≈ 4.93

4. Дисперсія та швидкість збіжності

Ключова думка: вдвічі зменшити похибку коштує у 4× більше вибірок (O(1/ε²) за похибки ε). Це гірше за детерміновані методи для гладких 1D-інтегралів (які збігаються як 1/N^k для деякого k ≥ 2), але краще за детерміновані методи у високих вимірах (де розмір сітки вибухає).

Зменшення дисперсії без додаткових вибірок

Кілька технік зменшують σ(f) без збільшення N:

• Контрольні змінні: відняти корельовану функцію g з відомим інтегралом, потім додати її назад аналітично.
• Антитетичні змінні: поєднати кожну вибірку x з її дзеркальним відображенням 1−x; їхні значення функції від’ємно корельовані.
• Стратифікована вибірка: поділити область на підобласті й вибирати з кожної пропорційно.
• Квазі-Монте-Карло: замінити псевдовипадкові послідовності послідовностями з низькою розбіжністю (Halton, Sobol) для збіжності O(log(N)^d/N).

5. Вибірка за значущістю

Якщо f(x) велика лише в малій області, більшість рівномірних вибірок майже нічого не дають оцінці. Вибірка за значущістю концентрує вибірки там, де f велика, вибираючи з альтернативного розподілу p(x) та коригуючи відношенням правдоподібності:

Відношення w(x) = f(x)/p(x) називається вагою значущості. Якщо p(x) точно повторює форму |f(x)|, ваги майже сталі, і дисперсія різко падає. Підступ: нам потрібна пропозиція p, з якої ми можемо вибирати й яка покриває всі області, де f не є нехтовно малою (інакше ваги можуть вибухнути — сумнозвісна проблема важких хвостів).

// Оцінюємо ∫₀^∞ x·exp(−x) dx = 1, використовуючи Exponential(1) як пропозицію
// p(x) = exp(−x),  f(x)/p(x) = x·exp(-x)/exp(-x) = x
function sampleExponential() {
  return -Math.log(1 - Math.random());   // метод оберненої функції розподілу
}

function importanceSamplingEstimate(N) {
  let sum = 0;
  for (let i = 0; i < N; i++) {
    const x = sampleExponential();   // x ~ Exp(1)
    sum += x;                         // вага f(x)/p(x) = x
  }
  return sum / N;   // збігається до 1 з нульовою дисперсією для цієї підінтегральної функції!
}
console.log(importanceSamplingEstimate(100));   // ≈ 1.000 (точно!)

6. Монте-Карло на ланцюгах Маркова (MCMC)

Коли цільовий розподіл p(x) складний (наприклад, багатовимірний апостеріорний розподіл у баєсівському висновуванні), ми не можемо вибирати з нього напряму. MCMC будує ланцюг Маркова, чий стаціонарний розподіл є p(x), а потім використовує траєкторію ланцюга як корельовані вибірки.

Алгоритм Метрополіса-Гастінгса

1. Почати з будь-якого стану x₀.
2. Запропонувати новий стан x′ із симетричної пропозиції q(x′|xᵢ) (наприклад, гаусів крок).
3. Обчислити коефіцієнт прийняття: α = min(1, p(x′)/p(xᵢ)).
4. Прийняти x′ з імовірністю α; інакше лишитися в xᵢ.
5. Повторювати, вибудовуючи траєкторію (x₀, x₁, x₂, …).

// Метрополіс-Гастінгс: вибірка з p(x) ∝ exp(−x²/2) (стандартний нормальний)
function metropolis(logP, start, stepSize, N) {
  const samples = [start];
  let current = start;
  let logPcur  = logP(current);

  for (let i = 1; i < N; i++) {
    const proposal = current + (Math.random() * 2 - 1) * stepSize;
    const logPprop  = logP(proposal);
    const logAlpha  = logPprop - logPcur;       // лог коефіцієнта прийняття

    if (Math.log(Math.random()) < logAlpha) {   // приймаємо
      current = proposal;
      logPcur = logPprop;
    }
    samples.push(current);
  }
  return samples;
}

// Ціль: стандартний нормальний   log p(x) = −x²/2  (стала відкинута)
const draws = metropolis(x => -x * x / 2, 0, 1, 50000);
// Після розігріву вибірки слідують N(0,1) — обчисліть середнє й дисперсію для перевірки
const mean = draws.reduce((a, b) => a + b) / draws.length;   // ≈ 0

MCMC потребує періоду розігріву (burn-in) (ранні вибірки автокорельовані й зміщені початковою точкою) та має ефективний розмір вибірки менший за N через автокореляцію. Практична діагностика включає статистику Гелмана-Рубіна R̂ та огляд трас-графіків.

7. JavaScript: практичні поради

Якість генератора псевдовипадкових чисел

Math.random() у V8 використовує алгоритм xorshift128+ — швидкий, 128-бітний стан, період 2¹²⁸ − 1. Він проходить усі стандартні статистичні тести й годиться для симуляцій. Для криптографічних застосунків використовуйте натомість crypto.getRandomValues().

Сідування для відтворюваності

// Простий LCG PRNG (не крипто-безпечний — лише для відтворюваних симуляцій)
function makePRNG(seed) {
  let s = seed & 0xffffffff;
  return function() {
    s = (Math.imul(1664525, s) + 1013904223) & 0xffffffff;
    return (s >>> 0) / 4294967296;
  };
}
const rng = makePRNG(42);
console.log(rng(), rng());   // та сама послідовність щоразу

Векторизовані пакети з TypedArrays

// Заповнюємо Float64Array N випадковими точками — швидше, ніж цикл Math.random()
function mcPiBatch(N) {
  const data = new Float64Array(N * 2);
  for (let i = 0; i < data.length; i++) data[i] = Math.random();
  let inside = 0;
  for (let i = 0; i < N; i++) {
    const x = data[2*i], y = data[2*i+1];
    if (x*x + y*y <= 1) inside++;
  }
  return 4 * inside / N;
}

8. Застосування в різних дисциплінах

Фізика частинок та симуляції високих енергій

Інструментарій GEANT4 симулює взаємодії частинок у детекторах за допомогою вибірки методом Монте-Карло перерізів, мод розпаду й кутів розсіювання. Одне зіткнення на LHC може породити мільйони вторинних частинок, усі з яких ймовірнісно відстежуються крізь матеріал детектора.

Оцінка вартості фінансових опціонів (Блек-Шоулз MC)

Ціна європейського опціону — це дисконтована очікувана виплата за ризик-нейтральною мірою. Оскільки траєкторія акції слідує геометричному броунівському руху, ми симулюємо N незалежних траєкторій, обчислюємо виплату для кожної й усереднюємо:

Баєсівське висновування

MCMC (зокрема гамільтонів Монте-Карло / NUTS, як у Stan та PyMC) дозволяє науковцям виводити апостеріорні розподіли параметрів за даними, навіть із сотнями корельованих параметрів. Це лежить в основі всього — від калібрування кліматичних моделей до геномних асоціативних досліджень.

Комп’ютерна графіка: трасування шляхів

Рендеринг фізично коректного зображення означає обчислення інтеграла рендерингу — вхідної яскравості, проінтегрованої по всіх напрямках і всіх світлових шляхах. Трасування шляхів оцінює цей інтеграл, випадково вибираючи відбиття світла. Сучасні рендерери (Cycles, Arnold, V-Ray) — усі є трасувальниками шляхів Монте-Карло.