MCMC Metropolis-Hastings у 50 рядках JavaScript

1. Що таке MCMC і навіщо він потрібен

Припустимо, ви знаєте густину ймовірності з точністю до константи, p(x) = f(x)/Z, де f(x) легко обчислити, а нормувальна константа Z = ∫ f(x) dx — нерозв'язна (типова ситуація в баєсівському висновуванні, де f — це апріорі × правдоподібність, а Z — маргінальна правдоподібність). Ви не можете вибрати з p(x) напряму, але все ще можете обчислити f(x) у будь-якій точці.

Ланцюг Маркова Монте-Карло (MCMC) будує ланцюг Маркова, стаціонарний розподіл якого дорівнює p(x) — випадкове блукання, яке після достатньої кількості кроків відвідує кожну область простору станів з точно правильною довгостроковою частотою, ніколи не потребуючи знати Z. Metropolis-Hastings (1953/1970) — класичний, універсальний рецепт побудови такого ланцюга для практично будь-якої цільової густини.

2. Цільовий та пропозиційний розподіли

Цільовий розподіл π(x) ∝ f(x) — це те, з чого ми хочемо отримати вибірки. Пропозиційний розподіл q(x'|x) генерує кандидата наступного стану x' з поточного стану x. Найпростіший і найпоширеніший вибір — симетрична пропозиція випадкового блукання:

Оскільки ця пропозиція симетрична (q(x'|x) = q(x|x')), вона спрощує формулу прийняття на наступному кроці — цей окремий випадок називається алгоритмом Metropolis. Metropolis- Hastings узагальнює його на асиметричні пропозиції з додатковим коригувальним членом.

3. Правило прийняття Metropolis-Hastings

На кожному кроці пропонуємо кандидата x' з q(x'|x), потім приймаємо його з імовірністю α, інакше залишаємось у x:

Для симетричної пропозиції q(x|x') = q(x'|x) скорочується, залишаючи простіше відношення Metropolis:

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

4. Реалізація у 50 рядках

Нижче logTarget(x) повертає логарифм ненормованої цільової густини (робота в логарифмічному просторі уникає underflow для мізерних імовірностей і перетворює відношення на віднімання). Приклад цілі — суміш двох гаусіанів — бімодальний розподіл, з якого було б незручно вибирати напряму методами оберненої функції розподілу на базі Math.random():

// Ненормований лог-густина: суміш N(-2, 1) та N(3, 0.7), ваги 0.4/0.6
function logTarget(x) {
  const logNormal = (x, mu, sigma) =>
    -0.5 * ((x - mu) / sigma) ** 2 - Math.log(sigma);
  const a = Math.log(0.4) + logNormal(x, -2, 1);
  const b = Math.log(0.6) + logNormal(x, 3, 0.7);
  // log-sum-exp для числової стабільності: log(e^a + e^b)
  const m = Math.max(a, b);
  return m + Math.log(Math.exp(a - m) + Math.exp(b - m));
}

// Вибірка стандартного нормального розподілу через перетворення Бокса-Мюллера
function gaussianStep(sigma) {
  const u1 = Math.random(), u2 = Math.random();
  return sigma * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
}

// Семплер Metropolis-Hastings (симетрична пропозиція випадкового блукання)
function metropolisHastings(logTarget, {
  nSamples = 20000, sigma = 1.0, x0 = 0,
} = {}) {
  const samples = new Array(nSamples);
  let x = x0;
  let logFx = logTarget(x);
  let accepted = 0;

  for (let i = 0; i < nSamples; i++) {
    const xProp   = x + gaussianStep(sigma);
    const logFxp  = logTarget(xProp);
    const logAlpha = logFxp - logFx;               // логарифм f(x')/f(x)

    if (Math.log(Math.random()) < logAlpha) {   // приймаємо якщо log(U) < log(alpha), alpha обмежене 1
      x = xProp;
      logFx = logFxp;
      accepted++;
    }
    // інакше: відхиляємо, залишаємось на поточному x (це повторення Є вибіркою)
    samples[i] = x;
  }
  return { samples, acceptanceRate: accepted / nSamples };
}

const { samples, acceptanceRate } = metropolisHastings(logTarget, { nSamples: 20000, sigma: 1.2 });
console.log('коефіцієнт прийняття:', acceptanceRate.toFixed(3)); // цільте на ~0.2-0.5

Це весь алгоритм: пропонуй, обчисли, прийми-або-відхили, повтори. Усе інше в цьому уроці — burn-in, проріджування, діагностика, багатовимірні цілі — це уточнення навколо цього ~30-рядкового ядра.

5. Burn-in і проріджування

Ланцюг починається з довільного x0, який може перебувати в області низької ймовірності, далеко від того, де насправді зосереджена маса цілі. Тому початкові вибірки не є репрезентативними для π(x) і мають бути відкинуті — це період burn-in.

function applyBurnInAndThin(samples, { burnIn = 2000, thin = 5 } = {}) {
  return samples
    .slice(burnIn)               // відкидаємо перші `burnIn` вибірок
    .filter((_, i) => i % thin === 0); // залишаємо кожну `thin`-ту вибірку
}

Проріджування — залишення лише кожної k-ї вибірки — зменшує сильну автокореляцію між послідовними вибірками ланцюга випадкового блукання (кожна вибірка близька до попередньої) ціною викидання більшості сирих вибірок. Це не збільшує справжній обсяг інформації на залишену вибірку, але суттєво спрощує подальший код, що припускає приблизно незалежні вибірки (наприклад, наївні гістограмні оцінки густини чи прості зведені статистики).

6. Діагностика: trace-графіки та автокореляція

Ніколи не довіряйте запуску MCMC наосліп — завжди перевіряйте, що він справді зійшовся і добре перемішався:

• Trace-графік: побудуйте значення вибірки проти номера ітерації. Добре перемішаний ланцюг виглядає як шумне «пухнасте гусеня», що коливається навколо носія цілі; погано перемішаний показує довгі плоскі ділянки (ланцюг застряг) або повільний дрейф (усе ще в burn-in).
• Коефіцієнт прийняття: занадто високий (> 70%) зазвичай означає, що крок пропозиції σ занадто малий — ланцюг повзе й повільно перемішується. Занадто низький (< 10%) означає, що σ занадто великий — більшість пропозицій відхиляється. Класична ціль для 1D Metropolis із випадковим блуканням — приблизно 20–50%.
• Функція автокореляції (ФАК): вимірює, наскільки корельована вибірка i з вибіркою i+k. Повільно спадна ФАК означає великий штраф для ефективного розміру вибірки.

// Автокореляція лагу k сирого (непроріджeного) ланцюга
function autocorrelation(samples, maxLag = 50) {
  const n = samples.length;
  const mean = samples.reduce((a,b) => a+b, 0) / n;
  const variance = samples.reduce((s,x) => s + (x-mean)**2, 0) / n;
  const acf = [];
  for (let k = 0; k <= maxLag; k++) {
    let cov = 0;
    for (let i = 0; i < n - k; i++) cov += (samples[i]-mean) * (samples[i+k]-mean);
    acf.push(cov / (n - k) / variance);
  }
  return acf; // acf[0] === 1 завжди; дивіться, як швидко спадає до 0
}

// Ефективний розмір вибірки (приблизна оцінка через суму ФАК)
function effectiveSampleSize(samples) {
  const acf = autocorrelation(samples, 100);
  let sum = 1;
  for (let k = 1; k < acf.length; k++) {
    if (acf[k] < 0.05) break;   // зупиняємось, коли автокореляція незначна
    sum += 2 * acf[k];
  }
  return Math.round(samples.length / sum);
}

7. Розширення на кілька вимірів і налаштування кроку

Базовий алгоритм узагальнюється напряму на вектор стану x ∈ ℝᵈ — замінюємо скалярний гаусів крок багатовимірним нормальним кроком, а logTarget — функцією від вектора:

function gaussianStepVec(dim, sigma) {
  return Array.from({ length: dim }, () => gaussianStep(sigma));
}

function metropolisHastingsND(logTarget, dim, { nSamples = 20000, sigma = 0.5 } = {}) {
  let x = new Array(dim).fill(0);
  let logFx = logTarget(x);
  const samples = [];
  for (let i = 0; i < nSamples; i++) {
    const step = gaussianStepVec(dim, sigma);
    const xProp = x.map((v, j) => v + step[j]);
    const logFxp = logTarget(xProp);
    if (Math.log(Math.random()) < logFxp - logFx) { x = xProp; logFx = logFxp; }
    samples.push([...x]);
  }
  return samples;
}

Зі зростанням розмірності d наївний ізотропний крок випадкового блукання стає дедалі неефективнішим — коефіцієнт прийняття падає, якщо σ не зменшується приблизно як 1/√d (класичний результат оптимального масштабування Metropolis із випадковим блуканням). У високих розмірностях просунутіші семплери, як-от Hamiltonian Monte Carlo (HMC) чи No-U-Turn Sampler (NUTS) — що використовуються в Stan і PyMC — застосовують градієнтну інформацію для набагато ефективніших ходів. Metropolis-Hastings залишається базовим будівельним блоком, що пояснює, чому вони всі працюють.