ML: Баєсівські методи

Вступ

У задачах машинного навчання та алгоритмах штучного інтелекту важливу роль відіграє формула Баєса. Документ є продовженням вступу до теорії ймовірностей і присвячений різним аспектам застосування цієї формули.

Формула Баєса

Нехай є $n$ несумісних (неперетинних) подій $\{H_1,...,H_n\}$, $H_i\,\cap\,H_j=\varnothing$, які далі називатимуться гіпотезами. Передбачається, що одна з них істинна (обов'язково відбувається у випробуванні): $$ \sum_i P(H_i) = 1. $$ У найпростішому випадку може бути дві несумісні гіпотези, які полягають у тому, що деяке твердження $D$ або істинне, або хибне:

$H_1: $ «людина хвора на хворобу $D$,    $H_2: $ «людина не хвора на хворобу $D$».

Будемо вважати, що відомі ймовірності гіпотез $P(H_i)$. Їх можна отримати зі спостережень (вище — це частка хворих у популяції — objectivist interpretation) або це може бути ступінь «переконаності» або «віри» в істинність гіпотези тим, хто приймає деяке рішення (subjectivist interpretation). Ці ймовірності називають апріорними («до досвідні», лат. «a priori»).

Коли настає деяка подія $E$ (що несе нову інформацію), вона може змінити ймовірності гіпотез на апостериорні ймовірності $P(E\to H_i)$ (отримані з досвіду, лат. «a posteriori»). Вони визначаються формулою Баєса: $$ P(E\to H_i) ~=~ \frac{P(E,H_i)}{P(E)} ~=~ \frac{P(H_i)\,P(H_i\to E)}{P(E)}. $$ Ймовірність події $P(E)$ для повних (що покривають усі випадки) і несумісних гіпотез обчислюється за формулою повної ймовірності: $$ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~P(E) ~=~ \sum_j P(H_j)\,P(H_j\to E) ~=~ \sum_j P(E,\,H_j). $$ У результаті сума $P(E\to H_i)$ по всіх $H_i$ дорівнює одиниці (істинна одна з гіпотез).

За формулою Баєса апостериорна ймовірність $P(E\to H_i)$ пропорційна апріорній ймовірності $P(H_i)$, помноженій на «функцію правдоподібності» (likelihood) $P(H_i\to E)$ події $E$ для гіпотези $H_i$. Зокрема, якщо ймовірність події не залежить від істинності гіпотези $P(H_i\to E)=P(E)$, то апріорна ймовірність не змінюється: $P(E\to H_i)=P(H_i)$.

Часто накопичити статистику для обчислення ймовірностей $P(H_i\to E)$ простіше, ніж для $P(E\to H_i)$. Наприклад, у випадку $H:$ «людина хвора на дану хворобу», можна взяти групу хворих і обчислити по ній імовірність $P(H\to E)$ для симптому $E$. Водночас симптом $E$ може виникати з безлічі найрізноманітніших причин (хвороб), відмінних від хвороби $H$.

Крім цього ймовірності $P(H\to E)$ часто стабільніші, ніж $P(E\to H)$. Наприклад, нехай почалася епідемія. Тоді ймовірність $P(H)$ зростає, а $P(H\to E)$ залишаться незмінними. За формулою Баєса $P(E\to H)$, звичайно, також збільшиться.

Приклад: Тестування хвороби

$p$ ${n}$ Tot
$d$ 9 1 10
$\bar{d}$ 90 9900 9990
Tot 99 9901 10000
Нехай у людини провели тестування з виявлення деякої хвороби $D$. Відомо, що цей тест у $90\%$ випадків дає позитивний результат, якщо людина дійсно хвора. Якщо ж людина здорова, тест помилково видає позитивний результат у 0.9% випадків. Нехай також відомо, що на 10'000 осіб цією хворобою хворіє 10 осіб.

Будемо замість подій говорити про випадкові величини, що приймають два значення
Гіпотеза $H=\{d,\bar{d}\}$, де $d$ — «людина хвора на хворобу $D$» і $\bar{d}$ — «людина не хвора на цю хворобу»;
Подія $~E=\{p,n\}$, де $p$ — тест позитивний, $n$ — тест негативний: $$ P(d\to p) = 0.9,~~~~~~~P(\bar{d}\to p) = 0.009,~~~~~~P(d) = 0.001. $$ Праворуч у таблиці наведено статистику (в особах), на підставі якої отримано оцінки для цих імовірностей. Значення $P(H\to E)$ обчислюються по рядках, а $P(E\to H)$ — по стовпцях: $P(d\to p)=9/10$, $P(p \to \bar{d})=90/99$.

Ймовірність виявити позитивний тест (незалежно від стану людини) дорівнює: $$ P(p ) = P(d)\,P(d\to p ) ~+~ P(\bar{d})\,P(\bar{d}\to p ) = 0.9*0.001 + 0.009*(1-0.001)= 0.0099 $$ Звичайно це ж значення можна отримати і з таблиці (tot у першій колонці 99/10000).

За формулою Баєса апостериорна ймовірність виявитися хворим дорівнює: $$ P(p \to d) = \frac{P(d)\,P(d\to p )}{P( p )} = \frac{0.001\cdot 0.9}{0.00989} = 0.091. $$ З таблиці цей результат отримується по першій колонці: 9/99.

Цікаво, що, хоча тест досить точно працює на здорових людях (помиляється $9$ разів на $1000$ тестів), ймовірність виявитися хворим при одноразовому позитивному тесті дорівнює лише $9\%$. Цей, на перший погляд, дивний результат пов'язаний з тим, що, хоча помилок у тесту мало, частка хворих серед усієї популяції ще менша. Варто проаналізувати таблицю, наведену вище (якщо тест провести у всій популяції позитивний результат буде у 99 осіб, хоча серед них лише 9 хворих).

☝ Зазвичай зручніше нормувати апостериорні ймовірності гіпотез наприкінці, опускаючи при обчисленнях однакові множники, які від гіпотез не залежать. Для прикладу з тестами це виглядатиме так (знак $\sim$ означає, що величини пропорційні, тобто рівні з точністю до загального множника, незалежного від $d$, $\bar{d}$): $$ \begin{array}{lclcl} P( p \to d) &\sim& P(d) \, P(d\to p ) &=& 0.9\cdot 0.001, \\ P( p \to \bar{d}) &\sim& P(\bar{d})\, P(\bar{d}\to p ) &=& 0.009\cdot 0.999. \end{array} $$ Запишемо $P( p \to H)$ як вектор, компоненти якого відповідають $H=d$ і $H=\bar{d}$: $$ P( p \to H) ~\sim~ \bigr\{ 0.9\cdot 0.001,~~~~ 0.009\cdot 0.999\bigr\} ~\sim~ \bigr\{ 0.1,~~~~ 0.999\bigr\}, $$ де з компонент знову винесені та опущені загальні множники. Тепер вектор необхідно пронормувати на одиницю: $P( p \to H)=\bigr\{ 0.1,~~~~ 0.999\bigr\}/(0.1+0.999)=\bigr\{ 0.091,~~~~ 0.909\bigr\} $.

Послідовність спостережень

При прийнятті рішень формулу Баєса використовують для перерахунку ймовірностей у міру надходження нової інформації. Нехай було проведено два незалежні тести (послідовно або одночасно, але в різних лабораторіях). Якими стануть апостериорні ймовірності в цьому випадку?

Позначимо через $E_1$, $E_2$ — дві умовно незалежні події (наприклад «pp» або «pn»). Ймовірність істинності гіпотези $H_i$ після виникнення подій за формулою Баєса дорівнює (вважаємо $\{E_1,\,E_2\}$ однією подією): $$ P(E_1,E_2 \to H_i) = \frac{P(H_i)\,P( H_i\to E_1,E_2 )}{P(E_1,E_2)} = \frac{P(H_i)\,P(H_i\to E_1)\,P(H_i\to E_2)}{P(E_1,E_2)}. $$ У другій рівності враховано, що події $E_1$ і $E_2$ умовно (при гіпотезі $H_i$) незалежні. Ймовірність у знаменнику, як і раніше, знаходиться з умови нормування $P(E_1,E_2 \to H_i)$: $$ P(E_1,E_2) ~=~ \sum_j P(H_j)\,P(H_j\to E_1)\,P(H_j\to E_2). $$

Звернімо увагу, що, хоча проведені тести умовно незалежні, але при цьому не безумовно незалежні: $P(E_1,E_2) \neq P(E_1)\,P(E_2)$, Чому так відбувається?

Відповідь Простіше зрозуміти це на прикладі тестування хвороби. Співвідношення $P(E_1,E_2) = P(E_1)\,P(E_2)$ було б вірним, якби результати тестів $E_1$, $E_2$ були б отримані у двох випадково вибраних людей. У нашому випадку тест проводиться з однією і тією ж людиною, яка або хвора, або здорова. Ми не знаємо, яка гіпотеза вірна, але в будь-якому випадку виконується одна з них. Тому $P(E_1,E_2) ~=~ P(h)\,P( h\to E_1,E_2)+P(\bar{h})\,P(\bar{h}\to E_1,E_2 )$.

Формулу для апостериорної ймовірності можна переписати таким чином: $$ P(E_1\to H_i) = \frac{P(H_i)\,P(H_i\to E_1 )}{P(E_1)},~~~~~~~~~~~~~~ P(E_1,E_2 \to H_i) = \frac{P(E_1\to H_i)\,P(H_i \to E_2)}{P(E_1\to E_2)}. $$ Таким чином, подія $E_1$ перетворює апріорну ймовірність $P(H_1)$ на $P(E\to H_1)$, після чого вона (за допомогою тієї ж формули Баєса) перетворюється на ще більш уточнену ймовірність гіпотези $P(E_1,E_2\to H_i)$ тощо. При цьому послідовність умовно незалежних подій не відіграє ролі.

Приклад обчислень на numpy

Наведемо приклад обчислень з тестуванням хвороби на Python за допомогою бібліотеки NumPy: 
import numpy as np
np.set_printoptions(precision=5, suppress=True)
Спільні ймовірності будемо починати з префікса P_, а умовні з C_.
Нумерація подій: E=[p,n] і гіпотез: H=[б,з]. Вихідні дані мають вигляд: 

NH, NE = 2, 2                          # число гіпотез і подій

P_H  = np.array([0.001, 1-0.001])      # [P(б), P(з)]

C_EH = np.array([[  0.9,   0.009],     # P(б -> p) P(з -> p)
                 [1-0.9, 1-0.009] ])   # P(б -> n) P(з -> n)
Результати одноразового тесту: 
C_HE = ( C_EH*P_H.reshape(1,NH) ).T    # [[P(p -> б), P(n -> б)],      [[0.09099 0.0001]
C_HE /= C_HE.sum(0)                    #  [P(p -> з), P(n -> з)]]   =   [0.90901 0.9999]]
Результати двох тестів: 
C_HEE = (C_EH.reshape(NE,1,NH)*C_EH.reshape(1,NE,NH)*P_H.reshape(1,1,NH)).transpose(2,0,1)
C_HEE /= C_HEE.sum(0)

# P(pp -> б)  P(pn -> б)  [[[0.90917 0.01   ]
# P(np -> б)  P(nn -> б)    [0.01    0.00001]]
#
# P(pp -> з)  P(pn -> з)   [[0.09083 0.99   ]
# P(np -> з)  P(nn -> з)    [0.99    0.99999]]]

Таким чином, імовірність виявитися хворим при двох позитивних тестах $P(pp\to \text{б})=0.909$ суттєво вища, ніж при одному: $P(p\to \text{б})=0.091$.

Задача класифікації

У задачі класифікації з неперетинними класами об'єкт з ознаками $\mathbf{x}=\{x_1,...,x_n\}$ необхідно віднести до одного з $C$ класів. Нехай $P(\mathbf{x} \to c)$ — умовна ймовірність того, що об'єкт з ознаками $\mathbf{x}$ належить $c$-тому класу: $c\in [0...C-1]$. Тоді формулу Баєса можна записати таким чином: $$ P(\mathbf{x} \to c) ~=~ \frac{p(c,\mathbf{x})}{p(\mathbf{x})} = \frac{P(c)\,p(c\to \mathbf{x})}{p(\mathbf{x})}, ~~~~~~~~~~~~~ p(\mathbf{x}) ~=~ \sum_\alpha P(\alpha)\,p(\alpha\to \mathbf{x}), $$ де $P(c)$ — ймовірності належності випадково вибраного об'єкта до класу $c$ (незалежно від значень його ознак), а $p(c\to \mathbf{x}) = p(c,\mathbf{x})/P(c)$ — розподіл щільності ймовірностей значень ознак $\mathbf{x}$ для об'єктів з даного класу $c$. Нормування в знаменнику $p(\mathbf{x})$ — це розподіл об'єктів у просторі ознак. Якщо є $N$ об'єктів, то $p(\mathbf{x})\cdot N$ — це щільність числа об'єктів на одиницю об'єму.

Таким чином, щоб отримати ймовірність $P(\mathbf{x}\to c)$, необхідно знати ймовірності класів $P(c)$ і щільність розподілу $p(c\to \mathbf{x})$. Якщо навчальна вибірка не зміщена (класи в ній представлені з такою ж частотою, як і в тестовій), то $P(c)$ отримати нескладно. А ось значення $p(c\to \mathbf{x})$ у багатовимірному просторі ознак часто знайти непросто. Якщо об'єкти класу утворюють компактні «кластери», то $p(c\to \mathbf{x})$ можна апроксимувати гладкими функціями (наприклад, багатовимірними розподілами Гаусса). Параметри цих функцій отримують по множині навчальних прикладів. Інший підхід — це пошук найближчих сусідів до точки $\mathbf{x}$ у кожному класі і побудова по них апроксимуючих поверхонь $p(c\to \mathbf{x})$.

Іноді непогані результати дає наївний Баєс, у якому передбачається умовна незалежність ознак: $$ p(c\to \mathbf{x}) ~=~p(c\to x_1)\cdot...\cdot p(c\to x_n). $$ Ймовірності $p(c\to x_i)$ особливо легко обчислюються для категоріальних або бінарних ознак. У цьому випадку такий класифікатор може виявитися ефективнішим за метричні моделі. Якщо ознаки неперервні, то для $p(c\to x_i)$ за навчальними даними можна будувати гістограми або обчислювати статистики (наприклад, коли вважається, що $p(c\to x_i)$ має нормальний розподіл).

При використанні наївного Баєса можуть виникати дві проблеми — втрата точності при перемноженні невеликих чисел і рівність нулю деяких імовірностей $p(c\to x_i)$ (якщо навчальних даних мало).
Перша проблема розв'язується переходом до логарифмів (замість добутку буде сума). Для боротьби з нульовими ймовірностями застосовують згладжування Лапласа.

Визначення частин мови

Баєсівську формулу можна використовувати для присвоєння кожному слову $w_i$ речення його частини мови $t_i$: іменник (N), прикметник (A), дієслово (V) тощо (Part-of-Speech Tagging). Нехай $t_i$ — один із маркерів (N, A, V...). Якщо для тексту слів $w_1...w_n$ відома умовна ймовірність $P(w_1,...,w_n\to t_1,...,t_n)$, то задача розв'язується пошуком максимуму цієї ймовірності: $$ t_1,...,t_n = \arg\max_{t_1,...,t_n}\,P(w_1...w_n\to t_1...t_n). $$ Для оцінки ймовірності запишемо формулу Баєса: $$ P(w_1...w_n\to t_1...t_n) ~=~ \frac{P(w_1...w_n,t_1...t_n)}{P(w_1...w_n)} ~=~ \frac{P(t_1...t_n)\,P( t_1...t_n \to w_1...w_n)}{P(w_1...w_n)} $$ Для максимізації по $t_i$ достатньо знайти лише чисельник. Ймовірність $P(t_1...t_n)$ характеризує ступінь правдоподібності даної послідовності частин мови (наприклад N V більш імовірне, ніж A V). У першому наближенні для її оцінки можна використовувати біграми. Для цього в ланцюговому правилі робиться марковське спрощення: $P(t_1...t_{k-1}\to t_k)=P(t_{k-1}\to t_k)$: $$ P( t_1...t_n) = P(t_1)\cdot P(t_1\to t_2)\cdot... \cdot P(t_{n-1}\to t_n). $$ Друга ймовірність у найпростішому випадку оцінюється за допомогою таких біграм: $$ P( t_1...t_n \to w_1...w_n) ~=~ P( t_1...t_n \to w_1)\cdot ... \cdot P( t_1...t_n \to w_1) ~=~ P(t_1\to w_1)\cdot...\cdot P(t_n\to w_n). $$ При цьому робиться припущення (достатньо грубе) про умовну незалежність слів у послідовності. Зауважимо, що за тією ж теоремою Баєса $P(t\to w)=P(w\to t)\,P(w)/P(t)$.

Вгадати задумане

Досить поширені інтелектуальні системи, завдання яких полягає у з'ясуванні намірів людини. Для цього система задає послідовність питань, на які людина відповідає Так або Ні. Відомим прикладом є гра Akinator, яка вгадує задуманого людиною персонажа.

Припустимо, людина задумала один з об'єктів $X=\{x_1,...,x_N\}$. Система задає питання $Q_i\Rightarrow \{y_i, n_i\}$, на яке отримує бінарну відповідь $y_i$ (yes=так) або $n_i$ (no=ні). Наприклад: «Ваш персонаж жінка?», «Ваш персонаж з реального життя?», «Ваш персонаж дружить з ведмедем?» тощо. При побудові подібної системи необхідно враховувати, що людина, відповідаючи, може навмисно або випадково давати невірні відповіді. Тому в загальному випадку від нуля відмінні обидві умовні ймовірності: $P(X\to y_i)$ і $P(X\to n_i)$ — отримати на питання $Q_i$ відповідь $y_i$ або $n_i$, якщо задумано об'єкт із множини $X$.

yes $y_1$ $y_2$ ... $y_M$ $P(x)$
$x_1$ 9 1 ... 2 0.001
$x_2$ 10 0 ... 25 0.008
... ... ... ... ... ...
$x_N$ 0 890 ... 100 0.02
Нехай є фіксований набір $N$ об'єктів і $M$ питань. У процесі роботи система повинна навчатися, запам'ятовуючи відповіді людей. У результаті виникає таблиця, подібна до наведеної праворуч (число відповідей $y_j$ на питання $Q_j$, коли задумано об'єкт $x_i$) і аналогічна таблиця для відповідей «no». На підставі цих даних для пар $(x_i,y_j)$ і $(x_i,n_j)$ можна отримати ймовірності $P(X\to Q)$ — по рядках таблиць і $P(Q\to X)$ — по їхніх стовпцях. Крім цього, таблиці дають імовірності $P(X)$ задуманих людьми персонажів.

Позначимо через $I_k$ послідовність $k$ питань-відповідей, отриманих системою: $I_k=\{Q_1,...,Q_k\}$.
На підставі цієї інформації формуються умовні ймовірності $P(I_k\to X)$ для кожного об'єкта з $X$.
За допомогою формули Баєса можна записати: $$ P(I_k\to X) = \frac{P(X)\,P(X\to I_k)}{P(I_k)} = \frac{P(X)\,P(X\to Q_1)...P(X\to Q_k)}{P(I_k)}, $$ де в другій рівності зроблено припущення (порівняно грубе) про «не вплив» пар питань-відповідей одна на одну. Нормувальна ймовірність $P(I_k)$, як зазвичай, знаходиться з умови нормування (сума по всіх $x_i$), за припущення, що людина задумала відомий системі об'єкт. Оскільки задача полягає в пошуку $x$ з максимальним значенням $P(I_k\to X)$, знаменник можна опустити: $$ x= \text{arg}\max_X ~\,P(X)\,P(X\to Q_1)...P(X\to Q_k). $$

Вибір чергового питання $Q$ після серії питань-відповідей $I_{k-1}$ проводиться аналогічно алгоритму ID3 побудови дерев рішень . Єдина відмінність полягає в тому, що «ознака» $Q$ може мати два значення $\{y,n\}$. Оптимальним є питання, яке робить розподіл ймовірностей $P=P(I_{k-1},\, Q_k\to X)$ максимально нерівномірним (підвищує ймовірності окремих об'єктів). Ентропія $H[P]$ такого розподілу повинна бути мінімальною (усереднюємо ентропії для відповідей «так» і «ні»): $$ Q = \text{arg}\min_{Q}~ H[ P(I_{k-1},\, y\to X)]\,P(y) + H[ P(I_{k-1},\, n\to X)]\,P(n), $$ де ймовірності $P(Q)$ знаходяться з урахуванням відомої на цей момент імовірності об'єкта $x$: $$ P(Q) = \sum_{x} P(X\to Q)\,P(I_{k-1}\to X). $$

Оцінка параметрів моделі

Часто необхідно оцінити деякий, наприклад, неперервний параметр $\theta$ за наявним спостереженням $x$. Зазвичай $\theta$ є параметром статистики (ймовірність події, параметри розподілу тощо).

Нехай відомий апріорний розподіл імовірностей $p(\theta)$ параметра $\theta$. При повній відсутності інформації про його значення розподіл є рівномірним $p(\theta)=\text{const}$. Якщо були попередні спостереження або існують деякі суб'єктивні міркування, функція розподілу $p(\theta)$ може бути більш визначеною.

Апостериорний розподіл параметра (після спостереження $x$) обчислюється, як зазвичай, за формулою Баєса: $$ p(x\to \theta) ~=~ \frac{p(\theta)\,p(\theta\to x)}{p(x)},~~~~~~~~~~~~~~~p(x) ~=~ \int p(\theta)\,p(\theta\to x)\,d\theta, $$ де інтегрування ведеться по всій області «дозволених» значень $\theta$. Для обчислення $p(x\to \theta)$, крім апріорного розподілу $p(\theta)$, необхідно знати функцію правдоподібності $p(\theta\to x)$. Найчастіше вона задається з теоретичних міркувань.

Коли параметр $\theta$ є багатовимірною векторною величиною $\{\theta_1,...,\theta_n\}$, обчислення $n$-вимірного інтеграла з потрібною точністю може бути утрудненим. Щоб уникнути подібних труднощів, вибирають таку $p(\theta)$, щоб інтеграл не лише обчислювався аналітично, але й отриманий розподіл $p(x\to \theta)$ при подальших спостереженнях також призводив до аналітичних результатів.

◊ Розглянемо як приклад визначення ймовірності події за спостереженням факту її настання.
У цьому випадку $x=1$ (подія в спостереженні настала) і $x=0$ (не настала). Виберемо як апріорний розподіл закон Бернуллі: $$ p(\theta) = \frac{\theta^{\alpha-1}\,(1-\theta)^{\beta-1}}{\mathrm{B}(\alpha,\beta)}. $$ Дійсні величини $\alpha,\beta$ є гіперпараметрами, задавання яких призводить до різних форм розподілу $p(\theta)$. Наприклад $\alpha=\beta=1$ — це рівномірний розподіл. При $\alpha=\beta=10$ — розподіл має симетричний максимум в околі $\theta=0.5$ тощо. Нормувальний множник $\mathrm{B}(\alpha,\beta)$ є бета-функцією, що виражається через гамма-функції: $$ \mathrm{B}(\alpha,\beta) = \int\limits_0^1 t^{\alpha-1}(1 - t)^{\beta-1}\,dt = \frac{\Gamma(\alpha) \Gamma(\beta)}{\Gamma(\alpha+\beta)}. $$

Якщо подія випадкова і володіє стабільною частотою, то ймовірність того, що вона відбудеться або не відбудеться, дорівнює: $$ p(\theta\to x) = \left\{ \begin{array}{lll} \theta & x=1\\ 1-\theta & x=0 \end{array} \right. $$

У загальному випадку, після $n+m$ спостережень, у яких подія відбулася $n$ разів і $m$ разів не відбулася, апостериорна ймовірність дорівнюватиме: $$ P(n,m\to \theta) = \frac{x^{n+\alpha-1}\,(1-x)^{m+\beta-1}}{\text{B}(n+\alpha,\,m+\beta)}. $$ Залежно від значення гіперпараметрів $\alpha,\beta$ і результатів спостережень отримується той чи інший розподіл для значення ймовірності $\theta$.