Bayes' theorem

Інтуїція

Тест, який виявляє 99% людей із хворобою, повертає позитивний результат. Більшість людей робить висновок, що вони, мабуть, хворі. Вони помиляються, і не трохи: відповідь тут — 7.4%.

Помилка в тому, що число про тест читають як число про пацієнта. «На жаль, p(A|B) не є тим самим, що p(B|A), але є спосіб перейти від одного до іншого» [1, §1.3], і цим способом є теорема Баєса. Те, що вона додає до точності тесту, є рівно тим, чого точність не каже: скільки людей мали хворобу ще до того, як когось протестували.

Теорема — це арифметика, якої ніхто не заперечує. Вартою сторінки її робить те, що її відповідь суперечить сильній інтуїції, у ситуаціях, які важать, і суперечність ця не є тонкістю. Це множник понад тринадцять.

Формальне визначення

Виведення — три рядки. Кон’юнкцію можна записати двома способами — p(A and B) = p(A) p(B|A) для будь-яких A і B [1, §1.2], і так само з переставленими літерами — тож два добутки рівні, а ділення на p(B) дає теорему [1, §1.4].

Версія, варта запам’ятовування, — це та, що у формулі вгорі, де H позначає гіпотезу, а D — дані. Це діахронічне тлумачення, назване так тому, що «щось відбувається в часі; у цьому випадку ймовірність гіпотез змінюється, з часом, у міру того як ми бачимо нові дані» [1, §1.5]. Кожен член має назву:

• p(H) — апріорна ймовірність, ймовірність гіпотези до даних.
• p(D|H) — правдоподібність, ймовірність даних, якщо гіпотеза істинна. «Правдоподібність зазвичай найлегше обчислити».
• p(D) — нормувальна константа, «ймовірність даних за будь-якої гіпотези взагалі», яку «у найзагальнішому випадку … важко прицвяхувати».
• p(H|D) — апостеріорна ймовірність, те, чого ви хотіли.

Усі чотири процитовано з [1, §1.5].

Нормувальна константа перестає бути важкою тієї миті, коли гіпотези є взаємовиключними — «щонайбільше одна гіпотеза в наборі може бути істинною» — і сукупно вичерпними — «інших можливостей немає» [1, §1.5]. Тоді закон повної ймовірності дає p(D) як суму по гіпотезах, і це просто підсумок стовпця.

Часто її можна пропустити цілком. Коли ви порівнюєте гіпотези проти тих самих даних, «знаменником є нормувальний член p(x), який не залежить від значення мітки y. Як наслідок, нам треба турбуватися лише про порівняння чисельника між різними значеннями y» [2, §22.9.3]. Кожен класифікатор, який оцінює prior × likelihood і бере найбільше, робить саме це.

Як це працює

Рисунок — це той самий скринінговий тест, намальований. Четверо людей із тисячі мають хворобу, і тест знаходить усіх, крім однієї сотої з них; п’ятдесят із решти 996 позначено однаково, бо 5% великого числа більші за 99% малого. П’ятдесят чотири позитивні результати, чотири з них справжні.

У самому тесті немає нічого несправного. Чутливість і специфічність є властивостями приладу; апріорна ймовірність є властивістю популяції, а апостеріорна є їхнім добутком. Змініть популяцію — тестуйте людей із симптомами, а не всіх — і той самий прилад дасть цілком іншу відповідь [3, §2.6.5].

Де застосовується

• Діагностика й скринінг, де базова частота зазвичай є панівним членом і зазвичай же тим, який у звіті опускають.
• Наївні баєсові класифікатори, які застосовують теорему до кожного класу, а тоді припускають, що ознаки умовно незалежні за відомої мітки. Без якогось такого припущення ланцюгове правило лишає «приблизно 2^d параметрів» для оцінювання [2, §22.9.3] — а модель, якій треба «побачити кожен можливий приклад, щоб передбачити відповідну мітку», «насправді не вивчає візерунок, а просто запам’ятовує набір даних» [2, §22.9.2].
• Фільтрування спаму, A/B-тести, злиття даних із датчиків — усюди, де свідчення надходять частинами, а вчорашня апостеріорна ймовірність стає сьогоднішньою апріорною.
• Баєсова регресія та ієрархічні моделі, де те саме оновлення прогонюють по неперервному параметру, а не по жменьці гіпотез. Див. баєсову лінійну регресію.
• Генеративні класифікатори, як-от GDA, які моделюють p(x|y) і обертають її.

Чому і коли

Беріть її щоразу, коли умовна ймовірність у вас не в той бік. У цьому й увесь трюк: «якщо p(A|B) важко обчислити або важко виміряти експериментально, перевірте, чи не легше обчислити інші члени теореми Баєса» [1, §1.8].

Не бійтеся апріорної ймовірності на тій підставі, що вона суб’єктивна. За достатньої кількості даних вона перестає важити — «за достатньої кількості даних люди, які починають із різних апріорних, схильні збігтися до однієї апостеріорної» [1, §4.3]. Практичний приклад міряє, як швидко.

А ось нульової апріорної бійтеся. «У баєсовому оновленні ми множимо кожну апріорну ймовірність на правдоподібність, тож якщо p(H) дорівнює 0, то p(H|D) теж дорівнює 0, хоч би яким було D» [1, §4.6]. Це правило Кромвеля — «уникайте давати апріорну ймовірність 0 будь-якій гіпотезі, яка є бодай віддалено можливою» [1, §4.6].

Беріть спряжену апріорну, коли оновлення має бути дешевим. Для біноміальної правдоподібності спряженим є бета-розподіл, тож «якщо апріорна — це бета-розподіл із параметрами alpha і beta, а ми бачимо дані з h орлами і t решками, то апостеріорна — це бета-розподіл із параметрами alpha+h і beta+t. Інакше кажучи, ми можемо зробити оновлення двома додаваннями» [1, §4.5].

Практичний приклад

Два питання, на обидва відповідають підрахунком, а не твердженням. Перше: наскільки позитивний результат є тестом, а наскільки — популяцією? Друге: як довго міцно тримана апріорна ймовірність переживає зіткнення з даними?

import numpy as np

rng = np.random.default_rng(0)

# --- 1. скринінговий тест, змодельований по одній людині за раз ------------
PREV, SENS, SPEC = 0.004, 0.99, 0.95
N = 2_000_000

sick = rng.random(N) < PREV
pos1 = np.where(sick, rng.random(N) < SENS, rng.random(N) >= SPEC)

print(f"  {N:,} people, prevalence {PREV:.1%}, sensitivity {SENS:.0%}, specificity {SPEC:.0%}")
print(f"    tested positive          {pos1.sum():>8,}")
print(f"    of those, actually sick  {(pos1 & sick).sum():>8,}")
print(f"    P(sick | positive)         {(pos1 & sick).sum() / pos1.sum():>8.4f}   counted")
bayes = PREV * SENS / (PREV * SENS + (1 - PREV) * (1 - SPEC))
print(f"                               {bayes:>8.4f}   Bayes's theorem")

# --- 2. ті самі люди проходять тест удруге ---------------------------------
pos2 = np.where(sick, rng.random(N) < SENS, rng.random(N) >= SPEC)
both = pos1 & pos2
print(f"\n    positive twice           {both.sum():>8,}")
print(f"    of those, actually sick  {(both & sick).sum():>8,}")
print(f"    P(sick | positive twice)   {(both & sick).sum() / both.sum():>8.4f}   counted")
b2 = bayes * SENS / (bayes * SENS + (1 - bayes) * (1 - SPEC))
print(f"                               {b2:>8.4f}   yesterday's posterior as today's prior")

# --- 3. троє аналітиків, які не згодні щодо монети -------------------------
TRUE_P = 0.35
priors = {
    "uniform      Beta(1, 1)": (1.0, 1.0),
    "sure it's ~0.9  B(30, 3)": (30.0, 3.0),
    "sure it's ~0.1  B(3, 30)": (3.0, 30.0),
}
flips = rng.random(10_000) < TRUE_P
print(f"\n  a coin with p = {TRUE_P}, three priors, posterior mean after n flips")
print("    " + "prior".ljust(26) + "".join(f"{n:>9}" for n in (0, 10, 100, 1_000, 10_000)))
for label, (a, b) in priors.items():
    row = []
    for n in (0, 10, 100, 1_000, 10_000):
        h = int(flips[:n].sum())
        row.append(f"{(a + h) / (a + b + n):>9.3f}")
    print("    " + label.ljust(26) + "".join(row))

# --- 4. аналітик, який виключив істину заздалегідь -------------------------
grid = np.linspace(0.0005, 0.9995, 1000)
flat = np.ones_like(grid)
closed = np.where(grid >= 0.5, 1.0, 0.0)  # "I am certain the coin favours heads"
print()
for prior, name in ((flat, "open mind    "), (closed, "p >= 0.5 only")):
    for n in (100, 10_000):
        h = int(flips[:n].sum())
        # у логарифмічному просторі: theta**3500 зникає в нуль у float64
        loglik = h * np.log(grid) + (n - h) * np.log(1 - grid)
        post = prior * np.exp(loglik - loglik.max())
        post = post / post.sum()
        print(f"    {name}  after {n:>6,} flips   posterior mean {(grid * post).sum():.3f}")

  2,000,000 people, prevalence 0.4%, sensitivity 99%, specificity 95%
    tested positive           107,359
    of those, actually sick     7,940
    P(sick | positive)           0.0740   counted
                                 0.0737   Bayes's theorem

    positive twice             12,776
    of those, actually sick     7,851
    P(sick | positive twice)     0.6145   counted
                                 0.6116   yesterday's posterior as today's prior

  a coin with p = 0.35, three priors, posterior mean after n flips
    prior                             0       10      100     1000    10000
    uniform      Beta(1, 1)       0.500    0.417    0.314    0.352    0.351
    sure it's ~0.9  B(30, 3)      0.909    0.791    0.459    0.370    0.353
    sure it's ~0.1  B(3, 30)      0.091    0.163    0.256    0.344    0.351

    open mind      after    100 flips   posterior mean 0.314
    open mind      after 10,000 flips   posterior mean 0.351
    p >= 0.5 only  after    100 flips   posterior mean 0.512
    p >= 0.5 only  after 10,000 flips   posterior mean 0.501

Арифметика й підрахунок збігаються, і саме заради цього його й прогонюють. 0.0737 проти 0.0740 на двох мільйонах людей: теорема Баєса тут не є ні наближенням, ні вибором моделі, вона є тим, що дає підрахунок.

Другий позитивний результат змінює все. Апостеріорна ймовірність іде з 0.0740 до 0.6145, бо оновлення прогнали з 0.0737 як апріорною, а не з 0.004. Це діахронічне прочитання, зроблене конкретним [1, §1.5], і саме тому повторне тестування є стандартною відповіддю на скринінг рідкісної хвороби.

Зауважте, що двічі-позитивна цифра трохи вища за ланцюгову арифметику (0.6145 проти 0.6116). Два тести змодельовано як незалежні за відомого справжнього статусу, а це саме те припущення, яке робить ланцюжок; розрив є вибірковим шумом на 12,776 людях, а не вадою котрогось із них.

Двоє аналітиків, які не могли б розходитися сильніше, опиняються в одному місці. Один починає з 0.909, другий з 0.091 — віра дев’ять до одного в протилежних напрямках. Після десяти тисяч підкидань вони на 0.353 і 0.351, проти істини 0.35. Жодного з них не переконали аргументом; їх переголосували дані. Це і є затоплення апріорних [1, §4.3], виміряне.

Єдина апріорна, яка ніколи не відновлюється, — це та, що виключила відповідь. Аналітик, упевнений, що монета схиляється до орла, усе ще на 0.501 після десяти тисяч підкидань монети, яка випадає орлом у 35% випадків — притиснутий до краю того, що він був готовий розглядати, і жодна кількість свідчень його не рухає, бо кожну гіпотезу, яка могла б пояснити дані, помножено на нуль [1, §4.6].

Типові помилки

Плутати p(D|H) з p(H|D). Тест має чутливість 99% і позитивний результат є інформативним на 7.4%, водночас, без жодної суперечності. Це найбільш наслідкова арифметична помилка в прикладній статистиці.

Викидати базову частоту. Це той член, якого ніхто не міряє і який усі припускають. Якщо поширеність невідома, то апостеріорна ймовірність невідома — точність тесту її не постачить.

Зчіплювати оновлення, які не є умовно незалежними. Практичний приклад прогонює два тести, незалежні за відомого справжнього статусу. Два прогони того самого аналізу на тій самій погано взятій пробі незалежними не є, і поводитися з ними так, ніби вони незалежні, роздуває апостеріорну ймовірність: якщо провести той самий тест двічі, «додаткового розуміння це не дає» [3, §2.6.5].

Нульова апріорна або набір гіпотез, який не є вичерпним. Обидва є тією самою помилкою: істина поза набором, а «якщо p(H) дорівнює 0, то p(H|D) теж дорівнює 0, хоч би яким було D» [1, §4.6]. Тож апостеріорна ймовірність скупчується на найменш хибній із дозволених гіпотез і повідомляє високу впевненість — як аналітик із практичного прикладу, що допускав лише p ≥ 0.5 і застряг на 0.501.

Читати апостеріорну ймовірність так, ніби модель певна. «Увесь баєсів аналіз спирається на рішення щодо моделювання. Якщо ми з вами оберемо не ту саму модель, ми можемо тлумачити дані по-різному» [1, §4.6]. Апріорні затоплює; хибну функцію правдоподібності — ні.

Множити сирі ймовірності на довгому наборі даних. Останній блок вище потребує свого переписування в лог-простір: theta ** 3500 є нулем у float64, і апостеріорна ймовірність виходить nan. Див. максимальну правдоподібність.