2.3: Максимальна вірогідність

Розділ 2.3a: Що таке ймовірність?

Оскільки всі підходи, описані в решті частини цієї глави, передбачають обчислення ймовірностей, я спочатку коротко опишу це поняття. Хорошим загальним оглядом ймовірності є Edwards (1992). Імовірність визначається як ймовірність, задана моделлю і набором значень параметрів, отримання певного набору даних. Тобто, враховуючи математичний опис світу, яка ймовірність того, що ми побачимо фактичні дані, які ми зібрали?

Щоб розрахувати ймовірність, ми повинні розглянути конкретну модель, яка, можливо, породила дані. Ця модель майже завжди матиме значення параметрів, які потрібно вказати. Ми можемо посилатися на цю задану модель (з конкретними значеннями параметрів) як гіпотезу, H. Імовірність тоді:

\[L(H|D)=Pr(D|H) \label{2.1}\]

Тут L і Pr означають ймовірність і ймовірність, D для даних і H для гіпотези, яка знову включає як розглянуту модель, так і набір значень параметрів. Символ | розшифровується як «дано», тому рівняння 2.1 можна прочитати як «ймовірність гіпотези з урахуванням даних дорівнює ймовірності даних, наведених гіпотезою». Іншими словами, ймовірність представляє ймовірність при заданій моделі та значеннях параметрів, що ми отримаємо дані, які ми насправді бачимо.

Для будь-якої заданої моделі використання різних значень параметрів, як правило, змінює ймовірність. Як ви могли здогадатися, ми віддаємо перевагу значенням параметрів, які дають нам найвищу ймовірність отримання даних, які ми бачимо. Таким чином, одним із способів оцінки параметрів за даними є знаходження значень параметрів, які максимізують ймовірність; тобто значення параметрів, які дають найбільшу ймовірність, і найбільшу ймовірність отримання даних. Ці оцінки потім називаються оцінками максимальної ймовірності (ML). У рамках ML ми припускаємо, що гіпотеза, яка найкраще підходить до даних, - це та, яка має найвищу ймовірність генерувати ці дані.

Для наведеного вище прикладу нам потрібно обчислити ймовірність як ймовірність отримання голови 63 з 100 перевертань ящірки, враховуючи якусь модель гортання ящірки. Загалом, ми можемо написати ймовірність для будь-якої комбінації H «успіхів» (сальто, що дають голови) з n випробувань. У нас також буде один параметр, p _H, який буде представляти ймовірність «успіху», тобто ймовірність того, що будь-який фліп спливе голови. Ми можемо обчислити ймовірність наших даних, використовуючи біноміальну теорему:

\[ L(H|D)=Pr(D|p)= {n \choose H} p_H^H (1-p_H)^{n-H} \label{2.2} \]

У наведеному прикладі n = 100 і H = 63, так:

\[ L(H|D)= {100 \choose 63} p_H^{63} (1-p_H)^{37} \label{2.3} \]

Малюнок 2.2. Поверхня ймовірності для параметра p _H, дана монета, яка була перевернута як голови 63 рази з 100. Зображення автора, може бути використано повторно за ліцензією CC-BY-4.0.

Ми можемо скласти графік ймовірності, L, як функція p _H (рис. 2.2). Коли ми це робимо, ми бачимо, що максимальне значення ймовірності p _H, яке ми можемо назвати $\ hat {p} _H$, знаходиться на $\ hat {p} _H = 0.63$. Це підхід «грубої сили» до пошуку максимальної ймовірності: спробуйте безліч різних значень параметрів і підберіть той з найбільшою ймовірністю. Ми можемо зробити це набагато ефективніше, використовуючи числові методи, як описано в наступних розділах цієї книги.

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

\[ \ln{L} = \ln{n \choose H} + H \ln{p_H}+ (n-H)\ln{(1-p_H)} \label{2.4} \]

Зверніть увагу, що перетворення природного журналу (ln) змінює наше рівняння з силової функції на лінійну функцію, яку легко вирішити. Ми можемо диференціювати:

\[ \frac{d \ln{L}}{dp_H} = \frac{H}{p_H} - \frac{(n-H)}{(1-p_H)}\label{2.5} \]

Максимум ймовірності являє собою пік, який ми можемо знайти, встановивши похідну $\ frac {d\ ln {L}} {DP_h} $ в нуль. Потім ми знаходимо значення p _H, яке вирішує це рівняння, яке буде нашою оцінкою $\ hat {p} _H$. Отже, у нас є:

\[ \begin{array}{lcl} \frac{H}{\hat{p}_H} - \frac{n-H}{1-\hat{p}_H} & = & 0\\ \frac{H}{\hat{p}_H} & = & \frac{n-H}{1-\hat{p}_H}\\ H (1-\hat{p}_H) & = & \hat{p}_H (n-H)\\ H-H\hat{p}_H & = & n\hat{p}_H-H\hat{p}_H\\ H & = & n\hat{p}_H\\ \hat{p}_H &=& H / n\\ \end{array} \label{2.6}\]

Зверніть увагу, що для нашого простого прикладу H/n = 63/100 = 0.63, що точно дорівнює максимальній ймовірності від цифри 2.2.

Оцінки максимальної ймовірності мають багато бажаних статистичних властивостей. Варто, однак, відзначити, що вони не завжди будуть повертати точні оцінки параметрів, навіть коли дані генеруються за фактичною моделлю, яку ми розглядаємо. Насправді параметри ML іноді можуть бути упередженими. Щоб зрозуміти, що це означає, нам потрібно формально ввести два нових поняття: упередженість і точність. Уявіть, що ми повинні були імітувати набори даних за деякою моделлю A з параметром a. Для кожного моделювання ми потім використовували ML для оцінки параметра $\ hat {a} $ для змодельованих даних. Точність нашої оцінки ML говорить нам, наскільки відрізняється, в середньому, кожен з наших оцінюваних параметрів $\ hat {a} _i$ один від одного. Точні оцінки оцінюються з меншою невизначеністю. Упередженість, з іншого боку, вимірює, наскільки близькі наші оцінки $\ hat {a} _i $ до істинного значення a. Якщо наша оцінка параметрів ML є упередженою, то середнє значення $\ hat {a} _i$ буде відрізнятися від істинного значення a. Не рідкість оцінки ML бути упередженими таким чином, що залежить від розміру вибірки, так що оцінки наближаються до істини в міру збільшення розміру вибірки, але можуть бути досить далекими в певному напрямку, коли кількість точок даних невелика порівняно з кількістю параметрів, що оцінюються.

У нашому прикладі перегортання ящірки ми оцінили значення параметра $\ hat {p} _H = 0.63$. Для конкретного випадку оцінки параметра біноміального розподілу наша оцінка ML, як відомо, є неупередженою. І ця оцінка відрізняється від 0,5 - що було нашим очікуванням при нульовій гіпотезі. Так чи справедлива ця ящірка? Або, як варіант, ми можемо відхилити нульову гіпотезу про те, що р _Н = 0,5? Щоб оцінити це, нам потрібно скористатися вибором моделі.

Розділ 2.3b: Тест на коефіцієнт ймовірності

Вибір моделі передбачає порівняння набору потенційних моделей та використання якогось критерію для вибору тієї, яка забезпечує «найкраще» пояснення даних. Різні підходи по-різному визначають «кращий». Спочатку я обговорю найпростіший, але і найбільш обмежений, з цих методик, тест на коефіцієнт ймовірності. Тести на коефіцієнт правдоподібності можна використовувати лише в одній конкретній ситуації: порівняти дві моделі, де одна з моделей є окремим випадком іншої. Це означає, що модель А точно еквівалентна більш складній моделі B з параметрами, обмеженими певними значеннями. Ми завжди можемо визначити простішу модель як модель з меншою кількістю параметрів. Наприклад, можливо, модель B має параметри x, y та z, які можуть приймати будь-які значення. Модель A така ж, як модель B, але з параметром z фіксованим на 0. Тобто A - це особливий випадок B, коли параметр z = 0. Це іноді описується як модель A вкладена в модель B, оскільки кожна можлива версія моделі A дорівнює певному випадку моделі B, але модель B також включає більше можливостей.

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

Для прикладу розглянемо ще раз наш приклад гортання ящірки. Одна з моделей полягає в тому, що ящірка «справедлива»: тобто ймовірність голів дорівнює 1/2. Інша модель може полягати в тому, що ймовірність голів є деяким іншим значенням p, яке може становити 1/2, 1/3 або будь-яке інше значення між 0 і 1. Тут остання (складна) модель має один додатковий параметр, p _H, в порівнянні з колишньою (простою) моделлю; проста модель є окремим випадком складної моделі, коли p _H = 1/2.

Для таких вкладених моделей можна обчислити статистику тесту коефіцієнта ймовірності як

\[ \Delta = 2 \cdot \ln{\frac{L_1}{L_2}} = 2 \cdot (\ln{L_1}-\ln{L_2}) \label{2.7}\]

Тут Δ - це тестова статистика коефіцієнта ймовірності, L ₂ - ймовірність більш складної (багатопараметричної) моделі, а L ₁ - ймовірність більш простої моделі. Так як моделі вкладені, ймовірність складної моделі завжди буде більше або дорівнює ймовірності простої моделі. Це прямий наслідок того, що моделі вкладені. Якщо ми знайдемо особливу ймовірність для простішої моделі, ми завжди можемо знайти ймовірність, рівну такій для складної моделі, встановивши параметри так, щоб складна модель була еквівалентною простій моделі. Таким чином, максимальна ймовірність для складної моделі буде або це значення, або деяке більш високе значення, що ми можемо знайти за допомогою пошуку простору параметрів. Це означає, що тестова статистика Δ ніколи не буде негативною. Насправді, якщо ви коли-небудь отримаєте статистику тесту з негативним коефіцієнтом ймовірності, щось пішло не так - або ваші розрахунки неправильні, або ви насправді не знайшли рішень ML, або моделі насправді не вкладені.

Для проведення статистичного тесту порівняння двох моделей ми порівняємо тестову статистику Δ з її очікуванням за нульовою гіпотезою. Коли розміри вибірки великі, нульовий розподіл тестової статистики коефіцієнта ймовірності слідує за ші-квадратом (²) розподілу зі ступенями свободи, рівною різниці в кількості параметрів між двома моделями. Це означає, що якби простіша гіпотеза була правдою, і один проводив цей тест багато разів на великих незалежних наборах даних, тестова статистика приблизно слідувала б цьому розподілу ^{σ 2}. Щоб відхилити простішу (нульову) модель, порівнюють статистику тесту з критичним значенням, отриманим із відповідного розподілу ^{σ 2}. Якщо тестова статистика більша за критичне значення, відхиляється нульова гіпотеза. В іншому випадку нам не вдасться відкинути нульову гіпотезу. У цьому випадку нам потрібно розглянути лише один хвіст тесту ², оскільки кожне відхилення від нульової моделі підштовхне нас до вищих значень Δ та до правого хвоста розподілу.

Для наведеного вище прикладу, наведеного вище, ми можемо обчислити ln-ймовірність за гіпотезою p _H = 0.5 як:

\[ \begin{array}{lcl} \ln{L_1} &=& \ln{\left(\frac{100}{63}\right)} + 63 \cdot \ln{0.5} + (100-63) \cdot \ln{(1-0.5)} \nonumber \\ \ln{L_1} &=& -5.92\nonumber\\ \end{array} \label{2.8}\]

Ми можемо порівняти це з ймовірністю нашої оцінки максимальної ймовірності:

\[ \begin{array}{lcl} \ln{L_2} &=& \ln{\left(\frac{100}{63}\right)} + 63 \cdot \ln{0.63} + (100-63) \cdot \ln{(1-0.63)} \nonumber \\ \ln{L_2} &=& -2.50\nonumber \end{array} \label{2.9}\]

Потім ми обчислюємо статистику тесту коефіцієнта ймовірності:

\[ \begin{array}{lcl} \Delta &=& 2 \cdot (\ln{L_2}-\ln{L_1}) \nonumber \\ \Delta &=& 2 \cdot (-2.50 - -5.92) \nonumber \\ \Delta &=& 6.84\nonumber \end{array} \label{2.10}\]

Якщо порівняти це з розподілом φ ² з одним d.f., ми виявимо, що P = 0,009. Оскільки це P-значення менше порогу 0,05, ми відкидаємо нульову гіпотезу і підтримуємо альтернативу. Робимо висновок, що це не справедлива ящірка. Як і слід було очікувати, цей результат узгоджується з нашою відповіддю з біноміального тесту в попередньому розділі. Однак підходи математично різні, тому два P-значення не ідентичні.

Хоча описані вище з точки зору двох конкуруючих гіпотез, тести на коефіцієнт ймовірності можуть бути застосовані до більш складних ситуацій з більш ніж двома конкуруючими моделями. Наприклад, якщо всі моделі утворюють послідовність зростаючої складності, з кожною моделлю окремий випадок наступної більш складної моделі, можна порівняти кожну пару гіпотез послідовно, зупиняючи перший раз статистику тесту несуттєво. Крім того, в деяких випадках гіпотези можуть бути розміщені в роздвоєному дереві вибору, і можна переходити від простих до складних моделей по певному шляху парних порівнянь вкладених моделей. Цей підхід зазвичай використовується для вибору моделей еволюції послідовності ДНК (Posada and Crandall 1998).

Розділ 2.3c: Інформаційний критерій Akaike (AIC)

Можливо, ви помітили, що описаний вище тест на коефіцієнт ймовірності має деякі обмеження. Спеціально для моделей, що включають більше одного параметра, підходи, засновані на тестах на коефіцієнт ймовірності, можуть зробити лише стільки. Наприклад, можна порівняти серію моделей, деякі з яких вкладені в інші, використовуючи впорядковану серію тестів коефіцієнта ймовірності. Однак результати часто сильно залежатимуть від того, в якому порядку проводяться тести. Крім того, часто ми хочемо порівняти моделі, які не є вкладеними, як того вимагають тести коефіцієнта ймовірності. З цих причин може бути корисним інший підхід, заснований на інформаційному критерії Akaike (AIC).

Значення АІК для конкретної моделі є простою функцією ймовірності L і кількості параметрів k:

\[AIC = 2k − 2\ln L \label{2.11}\]

Ця функція врівноважує ймовірність моделі і кількість параметрів, оцінених в процесі підгонки моделі до даних. Критерій АПК можна розглядати як визначення моделі, яка забезпечує найбільш ефективний спосіб опису закономірностей у даних з невеликою кількістю параметрів. Однак цей стенографічний опис АПК не відображає фактичного математичного та філософського обґрунтування рівняння (2.11). Насправді це рівняння не є довільним; натомість його точний компроміс між числами параметрів та різницею ймовірності журналу походить від теорії інформації (для отримання додаткової інформації див. Burnham and Anderson 2003, Akaike (1998)).

Вищенаведене вище рівняння AIC (2.11) справедливо лише для досить великих розмірів вибірки щодо кількості оцінюваних параметрів (для n зразків та k параметрів, n/k > 40). Більшість емпіричних наборів даних включають менше 40 незалежних точок даних на один параметр, тому слід використовувати невелику корекцію розміру вибірки:

\[ AIC_C = AIC + \frac{2k(k+1)}{n-k-1} \label{2.12}\]

Ця корекція штрафує моделі, які мають невеликі розміри вибірки щодо кількості параметрів; тобто моделі, де параметрів майже стільки ж, скільки точок даних. Як зазначають Burnham and Anderson (2003), ця корекція мало впливає, якщо розміри вибірки великі, і тому забезпечує надійний спосіб виправлення можливого зміщення в наборах даних будь-якого розміру. Я рекомендую завжди використовувати невелику корекцію розміру вибірки при розрахунку значень AIC.

Щоб вибрати серед моделей, можна порівняти їх оцінки A I _C, і вибрати модель з найменшим значенням. Простіше проводити порівняння в A I C c _{оцінками} між моделями, обчислюючи різницю, Δ A I C _c. Наприклад, якщо ви порівнюєте набір моделей, ви можете обчислити Δ A I _C для моделі i як:

\[ΔAIC_{c_i} = AIC_{c_i} − AIC_{c_{min}} \label{2.13}\]

де A I C _{c i} є A I C _c оцінка для моделі i і A I C _{c m i n} - це мінімальна оцінка A I _C для всіх моделей.

Як широке правило для порівняння значень A I C, будь-яка модель з Δ A, I _{C c i} менше чотирьох, приблизно еквівалентна моделі з найнижчим A I C _c значення. Моделі з Δ A I C _{c i} між 4 і 8 мають невелику підтримку в даних, тоді як будь-яка модель з Δ A I C _c i більше 10 можна сміливо ігнорувати.

Крім того, можна обчислити відносну підтримку для кожної моделі за допомогою ваг Akaike. Вага для моделі i в порівнянні з набором конкуруючих моделей розраховується як:

\[ w_i = \frac{e^{-\Delta AIC_{c_i}/2}}{\sum_i{e^{-\Delta AIC_{c_i}/2}}} \label{2.14} \]

Ваги для всіх розглянутих моделей складають 1, тому w _i для кожної моделі можна розглядати як оцінку рівня підтримки цієї моделі в даних порівняно з іншими моделями, що розглядаються.

Повертаючись до нашого прикладу гортання ящірки, ми можемо обчислити оцінки A I _C для наших двох моделей наступним чином:

\[ \begin{array}{lcl} AIC_1 &=& 2 k_1 - 2 ln{L_1} = 2 \cdot 0 - 2 \cdot -5.92 \\\ AIC_1 &=& 11.8 \\\ AIC_2 &=& 2 k_2 - 2 ln{L_2} = 2 \cdot 1 - 2 \cdot -2.50 \\\ AIC_2 &=& 7.0 \\\ \end{array} \label{2.15} \]

Наш приклад трохи незвичний у тому, що модель не має оцінених параметрів; це трапляється іноді, але не характерно для біологічних застосувань. Ми можемо виправити ці значення для нашого розміру вибірки, який у цьому випадку дорівнює n = 100 відворотів ящірки:

\[ \begin{array}{lcl} AIC_{c_1} &=& AIC_1 + \frac{2 k_1 (k_1 + 1)}{n - k_1 - 1} \\\ AIC_{c_1} &=& 11.8 + \frac{2 \cdot 0 (0 + 1)}{100-0-1} \\\ AIC_{c_1} &=& 11.8 \\\ AIC_{c_2} &=& AIC_2 + \frac{2 k_2 (k_2 + 1)}{n - k_2 - 1} \\\ AIC_{c_2} &=& 7.0 + \frac{2 \cdot 1 (1 + 1)}{100-1-1} \\\ AIC_{c_2} &=& 7.0 \\\ \end{array} \label{2.16} \]

Зверніть увагу, що в даному конкретному випадку виправлення не вплинуло на наші значення A I C, принаймні до одного знака після коми. Це пояснюється тим, що розмір вибірки великий щодо кількості параметрів. Зверніть увагу, що модель 2 має найменшу оцінку A I _C і, таким чином, є моделлю, яка найкраще підтримується даними. Відзначивши це, ми тепер можемо перетворити ці A I _C оцінки у відносну шкалу:

\[ \begin{array}{lcl} \Delta AIC_{c_1} &=& AIC_{c_1}-AIC{c_{min}} \\\ &=& 11.8-7.0 \\\ &=& 4.8 \\\ \end{array} \label{2.17} \]

\[ \begin{array}{lcl} \Delta AIC_{c_2} &=& AIC_{c_2}-AIC{c_{min}} \\\ &=& 7.0-7.0 \\\ &=& 0 \\\ \end{array} \]

Зверніть увагу, що Δ A I C _{c i} для моделі 1 більше чотирьох, що свідчить про те, що ця модель («справедлива» ящірка) має мало підтримки в даних. Це знову узгоджується з усіма результатами, які ми отримали до цього часу, використовуючи як біноміальний тест, так і тест на коефіцієнт ймовірності. Нарешті, ми можемо використовувати відносні оцінки AiCC для обчислення ваг Akaike:

\[ \begin{array}{lcl} \sum_i{e^{-\Delta_i/2}} &=& e^{-\Delta_1/2} + e^{-\Delta_2/2} \\\ &=& e^{-4.8/2} + e^{-0/2} \\\ &=& 0.09 + 1 \\\ &=& 1.09 \\\ \end{array} \label{2.18}\]

\[ \begin{array}{lcl} w_1 &=& \frac{e^{-\Delta AIC_{c_1}/2}}{\sum_i{e^{-\Delta AIC_{c_i}/2}}} \\\ &=& \frac{0.09}{1.09} \\\ &=& 0.08 \end{array} \]

\[ \begin{array}{lcl} w_2 &=& \frac{e^{-\Delta AIC_{c_2}/2}}{\sum_i{e^{-\Delta AIC_{c_i}/2}}} \\\ &=& \frac{1.00}{1.09} \\\ &=& 0.92 \end{array} \]

Наші результати знову узгоджуються з результатами тесту на коефіцієнт ймовірності. Відносна ймовірність несправедливої ящірки становить 0,92, і ми можемо бути цілком впевнені, що наша ящірка не чесний ласт.

Ваги AIC також корисні для іншої мети: ми можемо використовувати їх для отримання усереднених за моделлю оцінок параметрів. Це оцінки параметрів, які поєднуються між різними моделями, пропорційними підтримці цих моделей. Як приклад думки уявіть, що ми розглядаємо дві моделі, A і B, для конкретного набору даних. І модель A, і модель B мають однаковий параметр p, і саме цей параметр нас особливо цікавить. Іншими словами, ми не знаємо, яка модель є найкращою моделлю для наших даних, але те, що нам дійсно потрібно, це хороша оцінка p. Ми можемо зробити це за допомогою модельного усереднення. Якщо модель А має велику вагу АПК, то усереднена за моделлю оцінка параметрів для p буде дуже близька до нашої оцінки p за моделлю А; однак, якщо обидві моделі мають приблизно рівну підтримку, то оцінка параметрів буде близька до середнього показника двох різних оцінок. Усереднення моделі може бути дуже корисним у випадках, коли існує велика невизначеність у виборі моделі для моделей, які поділяють цікаві параметри. Іноді самі моделі не представляють інтересу, але їх потрібно розглядати як можливості; в цьому випадку усереднення моделі дозволяє оцінити параметри таким чином, який не так сильно залежить від нашого вибору моделей.