Методи Монте Карло: Отримання Правильної Відповіді Загадуванням

Відповідь на запитання шляхом кидання кісток

Метод Монте-Карло замінює обчислення, яке неможливо виконати експериментом, який можна повторити. Запишіть бажане значення як очікуване, виберіть випадкові зразки із відповідного розподілу та обчисліть середнє арифметичне. Закон великих чисел гарантує, що середнє значення сходиться до справжнього значення; центральна теорема обмеженості розповідає вам, наскільки швидко це відбувається.

Ідея виникла в Лос-Аламос у 1940-х роках. Станіслав Улам, одужуючи від хвороби, намагався визначити ймовірність того, що гра в одинокі шашки закінчиться перемогою — і зрозумів, що легше прорахувати сто карт, ніж перераховувати комбінаторику. Джона фон Нейман побачив це одразу, і Ніколас Мертоф запропонував назву коду, натхненну казино в Монако, де у тітки Улама було велике захоплення.

π з квадрата та кола

Це класичний приклад демонстрації, і той, що працює в симуляції на цьому сайті. Його суть – оцінка числа π. Візьміть одиничне квадратне та заглинуте у нього чвертсере круга. Випадково кидайте точки в квадрат. Ймовірність того, що точка потрапить всередину чвертькруга, дорівнює відношенню площ, тобто π/4. Отже, порахуйте кількість попадань, поділіть на кількість кидків і помножте на чотири.

призначено змінну inside = 0; для (i = 0; i < n; i++) { const x = Math.random(), y = Math.random(); якщо x * x + y * y <= 1 то inside++; // не потрібно обчислювати квадратний корінь } постійна piEstimate = 4 * inside / n; Це нічого конкретного не стосується кіл. Будь-який інтеграл – це площа, а будь-яка площа – це ймовірність того, що випадкова точка потрапить під криву — тому цей самий трирядковий цикл є загальним-призначенням інтегратором. Єдина річ, яка змінюється від задачі до задачі – це те, що ви зраховуєте і те, що рахуєте.

let inside = 0;
for (let i = 0; i < n; i++) {
  const x = Math.random(), y = Math.random();
  if (x * x + y * y <= 1) inside++;      // no sqrt needed
}
const piEstimate = 4 * inside / n;

Закон 1/√N та чому він одночасно жахливий і чудовий

Кожен бросок є Бернуллівським експериментом з ймовірністю успіху p = π/4 ≈ 0.7854, отже його дисперсія становить p(1 − p) ≈ 0.1685. Оцінювач – це середнє значення N таких експериментів, тому його стандартна похибка дорівнює √(p(1−p)/N), а оцінка π — в чотири рази більша за цю, має типову похибку

σ(π̂) = 4 · sqrt(p(1 − p) / N) ≈ 1.64 / sqrt(N) N = 10⁴ → типова похибка ≈ 0.016 (приблизно 2 правильні знаки) N = 10⁶ → типова похибка ≈ 0.0016 (приблизно 3 правильні знаки) N = 10⁸ → типова похибка ≈ 0.00016 (приблизно 4 правильні знаки) Це визначальна властивість Монте-Карло: похибка падає як 1/√N, тому щоб отримати один додатковий десятковий розряд, потрібно набагато більше зразків. Як спосіб обчислення π це безглуздо — серії Машина отримують більше знаків за кілька термів, ніж цей цикл отримує за мільярд кидків.

Але подивіться, чого немає в цій формулі: виміру. Детерміноване правило квадрантури на сітці потребує mᵈ точок, щоб розмістити m точок вздовж кожного з d осей, і його точність руйнується зі збільшенням d — прокляття розмірності. 1/√N Монте-Карло не залежить від d взагалі. У двох вимірах це жарт; у двадцяти або двісті вимірах інтеграла транспорту світла, фінансової винагороди чи статистичної механіки розподілу частинок, воно часто є єдиним, що працює. Кожен трасер шляхів, кожен оцінювач опціонів і кожен бейсіанський зразник заднього перерозподілу будується на цьому спостереженні.

σ(π̂) = 4 · sqrt(p(1 − p) / N) ≈ 1.64 / sqrt(N)

N = 10⁴    → typical error ≈ 0.016     (about 2 correct digits)
N = 10⁶    → typical error ≈ 0.0016    (about 3 correct digits)
N = 10⁸    → typical error ≈ 0.00016   (about 4 correct digits)

Зменшення дисперсії: єдиний важіль

Оскільки показник у 1/√N фіксований, єдиний спосіб отримати кращий результат за ту ж саму вартість – зменшити константу перед ним – дисперсію величини, яку ми підсумовуємо. Чотири стандартні техніки:

важливість семплювання здійснюйте зразок там, де інтеграл великий, а потім діліть його на щільність семплінгу (вагу). Класичний випадок: безпосереднє відбиття світла від джерела замість направлення променів у переважно чорне полусфери. шарове семплювання розділіть область на K шарів і зразок кожен, заборонивши клапти та прогалини чистого випадкового. антитетичні варіанти для кожного зразка u також використовуйте 1 − u. Негативна кореляція скасовує частину помилки. контрольні варіанти відніміть корельовану величину, середнє значення якої ви вже знаєте, і додайте це середнє значення назад.

Також є квазі-Monte Carlo, який повністю відмовляється від випадковості для послідовності з низькою дисперсією (Sobol, Halton) – детерміновані точки, які були спроектовані для більш рівномірного заповнення простору, ніж будь-яка випадкова величина. Для гладких, помірних багатовимірних інтегралів її помилка може наблизитися до 1/N замість 1/√N, що є величезним успіхом; ціна полягає в тому, що класичні межі похибок більше не застосовуються, і перевага зникає, коли зростає розмірність. Випадкове квазі-Monte Carlo (перемішування послідовності) повертає оцінку помилки.

importance sampling   sample where the integrand is large, then divide
                      by the sampling density (the weight). The classic
                      case: sampling a light source directly instead of
                      firing rays into a mostly black hemisphere.

stratified sampling   split the domain into K strata and sample each,
                      forbidding the clumps and gaps of pure randomness.

antithetic variates   for every sample u, also use 1 − u. The negative
                      correlation cancels part of the error.

control variates      subtract a correlated quantity whose exact mean
                      you already know, and add that mean back.

Когда нельзя напрямую отбирать выборку из распределения

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

α = min(1, π(x_new) / π(x_old)) // симметричное предложение принять → цепь двигается отклонить → цепь остается на месте, и x_old подсчитывается снова Нормализующая константа отменяется в отношении, что является сутью вопроса. Стоимость заключается в том, что последовательные выборки коррелированы, поэтому N MCMC шагов несет меньше информации, чем N независимых выборок, и необходимо отбрасывать начальный период «сгорания», пока цепь не забудет, где она начала.

α = min(1, π(x_new) / π(x_old))       // symmetric proposal

accept  → the chain moves
reject  → the chain stays where it is, and you count x_old AGAIN

Випадкові числа самі по собі

Результат Монте-Карло є таким же хорошим, як і генератор випадкових чисел, який ви використовуєте. Три практичні правила. По-перше, використовуйте добре перевірений PRNG – сучасний на основі лічильника або сімейство xorshift, а не домашній лінійний конгруентальний генератор, низькі біти якого відомі як не випадкові та які здебільшого падають у гіперплощини. По-друге, явно підсиджуйте його і записуйте насіння: запуск Монте-Карло, який не можна відтворити, неможливо виправити. По-третє, стежте за періодом та кореляцією між потоками, якщо ви паралелізуєте – два працівники випадково діляться потоком є безшумною та катастрофічною помилкою, яку жоден тест не виявить, оскільки відповідь все ще здається правдоподібною.

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

Frequently asked questions

Спробуйте наживо

Усе, що вище, працює прямо у вашому браузері — відкрийте Monte Carlo π і змінюйте параметри під час роботи. Нічого не встановлюється, нічого не завантажується на сервер, уся модель живе в одній вкладці.