2.4: Байєсова статистика

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

\[ Pr(H|D) = \dfrac{Pr(D|H) \cdot Pr(H)}{Pr(D)} \label{2.19}\]

Перевага байєсівських підходів полягає в тому, що вони дозволяють оцінити ймовірність того, що гіпотеза вірна з урахуванням спостережуваних даних,\(Pr(H|D)\). Це дійсно така ймовірність, що більшість людей мають на увазі, коли вони думають про цілі свого дослідження. Однак теорема Байєса також виявляє вартість такого підходу. Поряд з ймовірністю\(Pr(D|H)\), потрібно також включити попередні знання про ймовірність того, що будь-яка дана гіпотеза вірна - P r (H). Це являє собою попереднє переконання, що гіпотеза вірна, ще до розгляду наявних даних. Ця попередня ймовірність повинна бути чітко кількісно визначена у всіх байєсівських статистичних аналізах. На практиці вчені часто прагнуть використовувати «неінформативні» пріори, які мало впливають на задній розподіл - хоча навіть термін «неінформативний» може бути заплутаним, оскільки пріор є невід'ємною частиною байєсового аналізу. Термін також\(Pr(D)\) є важливою частиною теореми Байєса і може бути розрахований як ймовірність отримання даних, інтегрованих над попередніми розподілами параметрів:

\[Pr(D)=∫_HPr(H|D)Pr(H)dH \label{2.20}\]

Однак\(Pr(D)\) є постійним при порівнянні відповідності різних моделей для заданого набору даних і, таким чином, не має впливу на вибір байєсової моделі за більшості обставин (і всі приклади в цій книзі).

У нашому прикладі гортання ящірки ми можемо зробити аналіз в байєсівських рамках. Для моделі 1 вільних параметрів немає. Через це\(Pr(H)=1\) і\(Pr(D|H)=P(D)\), так що\(Pr(H|D)=1\). Це може здатися дивним, але результат означає, що наші дані не впливають на структуру моделі. Ми нічого не дізнаємося про модель без вільних параметрів, збираючи дані!

Якщо розглядати модель 2 вище, то параметр p _H повинен бути оцінений. Ми можемо встановити рівномірний попередній між 0 і 1 для p _H, так що f (p _H) =1 для всіх p _H в інтервалі [0,1]. Ми також можемо написати це як «наш\(p_h\) попередній для U (0,1)». Потім:

\[ Pr(H|D) = \frac{Pr(D|H) \cdot Pr(H)}{Pr(D)} = \frac{P(H|p_H,N) f(p_H)}{\displaystyle \int_{0}^{1} P(H|p_H,N) f(p_h) dp_H} \label{2.21} \]

Далі відзначимо, що\(Pr(D|H)\) це ймовірність наших даних, заданих моделлю, яка вже зазначена вище як Equation\ ref {2.2}. Підключивши це до нашого рівняння, ми маємо:

\[ Pr(H|D) = \frac{\binom{N}{H} p_H^H (1-p_H)^{N-H}}{\displaystyle \int_{0}^{1} \binom{N}{H} p_H^H (1-p_H)^{N-H} dp_H} \label{2.22} \]

Це потворне рівняння насправді спрощує бета-розподіл, який можна виразити більш просто як:

\[ Pr(H|D) = \frac{(N+1)!}{H!(N-H)!} p_H^H (1-p_H)^{N-H} \label{2.23} \]

Ми можемо порівняти цей задній розподіл нашої оцінки параметрів, p _H, враховуючи дані, з нашим рівномірним попереднім (рис. 2.3). Якщо ви оглянете цей графік, ви побачите, що задній розподіл сильно відрізняється від попереднього — тобто дані змінили наше уявлення про значення, які повинні приймати параметри. Знову ж таки, цей результат якісно узгоджується як з описаними вище підходами частот, так і ML. У цьому випадку ми бачимо з заднього розподілу, що ми можемо бути цілком впевнені, що наш параметр p _H не дорівнює 0,5.

Як видно з цього прикладу, теорема Байєса дозволяє об'єднати наше попереднє переконання про значення параметрів з інформацією з даних, щоб отримати задній. Ці задні розподіли дуже легко інтерпретувати, оскільки вони виражають ймовірність параметрів моделі з урахуванням наших даних. Однак ця ясність приходить за ціною вимагати явного попереднього. Пізніше в книзі ми дізнаємося, як використовувати цю функцію Байєсівська статистика на нашу користь, коли ми насправді маємо деякі попередні знання про значення параметрів.

Малюнок 2.3. Байєсівські попередні (пунктирна лінія) і задній (суцільна лінія) розподіли для перегортання ящірки. Зображення автора, може бути використано повторно за ліцензією CC-BY-4.0.

Розділ 2.4b: Байєський MCMC

Іншим основним інструментом у наборі інструментів байєсівських порівняльних методів є використання інструментів Маркова Монте-Карло (MCMC) для обчислення задніх ймовірностей. Методи MCMC використовують алгоритм, який використовує «ланцюжок» обчислень для вибірки заднього розподілу. MCMC вимагає обчислення ймовірностей, але не складної математики (наприклад, інтеграції розподілу ймовірностей, як у Equation\ ref {2.22}, і тому являє собою більш гнучкий підхід до байєсовських обчислень. Часто інтеграли в Equation\ ref {2.21} є нерозв'язними, тому найефективнішим способом підходити до байєсовських моделей є використання MCMC. Крім того, налаштування MCMC, на мій досвід, простіше, ніж люди очікують!

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

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

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

Наступний алгоритм використовує алгоритм Метрополіса-Гастінгса для проведення байєсівського аналізу MCMC з одним вільним параметром.

Алгоритм Метрополіса-Гастінгса

1. Отримати початкове значення параметра.
   • Вибірка початкового значення параметра, p ₀, з попереднього розподілу.
2. Починаючи з i = 1, запропонуйте новий параметр для покоління i.
   • З огляду на поточне значення параметра, p, виберіть нове запропоноване значення параметра, p ′, використовуючи густину пропозиції Q (p ′| p).
3. Обчисліть три співвідношення.
   • а Коефіцієнт попереднього коефіцієнта. Це співвідношення ймовірності витягування значень параметрів p і p ′ від попереднього. \[R_{prior} = \frac{P(p')}{P(p)} \label{2.24}\]
   • б. співвідношення щільності пропозиції. Це співвідношення ймовірності пропозицій, що йдуть від p до p ′ і зворотного. Часто ми цілеспрямовано будуємо щільність пропозиції, яка є симетричною. Коли ми це робимо, Q (p ′| p) = Q (p | p ′) і a ₂ = 1, спрощуючи розрахунки. \[R_{proposal} = \frac{Q(p'|p)}{Q(p|p')} \label{2.25}\]
   • c Коефіцієнт ймовірності. Це співвідношення ймовірностей даних з урахуванням двох різних значень параметрів. \[R_{likelihood} = \frac{L(p'|D)}{L(p|D)} = \frac{P(D|p')}{P(D|p)} \label{2.26}\]
4. Помножити. Знайдіть добуток попередніх коефіцієнтів, коефіцієнт щільності пропозиції та коефіцієнт ймовірності. \[R_{accept} = R_{prior} ⋅ R_{proposal} ⋅ R_{likelihood} \label{2.27}\]
5. Прийняти або відхилити. Намалюйте випадкове число x з рівномірного розподілу між 0 і 1. Якщо x < R _{a c c e p t}, прийміть запропоноване значення p ′ (p _i = p ′); в іншому випадку відхилити і зберегти поточне значення р _{(р і = р).}
6. Повторити. Повторіть кроки 2-5 велику кількість разів.

Виконуючи ці дії, отримується набір значень параметрів p _i, де i - від 1 до загальної кількості поколінь в MCMC. Як правило, ланцюг має період «вигоряння» на початку. Це час до того, як ланцюг досяг стаціонарного розподілу, і його можна спостерігати, коли значення параметрів показують тенденції через час, а ймовірність для моделей ще не плато. Якщо усунути цей період «прогорання», то, як було розглянуто вище, кожен крок ланцюга - проба з заднього розподілу. Ми можемо узагальнити задні розподіли параметрів моделі різними способами; наприклад, шляхом обчислення засобів, 95% довірчих інтервалів або гістограм.

Ми можемо застосувати цей алгоритм до нашого прикладу перегортання монети. Ми будемо розглядати той же попередній розподіл, U (0, 1), для параметра p. Визначимо також щільність пропозиції, Q (p ′| p) U (p −, p +). Тобто ми додамо або віднімемо невелике число (≤ 0,01) для генерації запропонованих значень p ′ заданого p.

Для початку алгоритму малюємо значення p від попереднього. Скажімо для ілюстративних цілей, що значення, яке ми малюємо, дорівнює 0,60. Це стає нашим поточним параметром оцінки. Для другого кроку ми пропонуємо нове значення, p ′, спираючись на наш розподіл пропозицій. Ми можемо використовувати = 0,01, тому розподіл пропозицій стає U (0,59, 0,61). Припустимо, що наше нове запропоноване значення p ′=0.595.

Потім ми обчислюємо наші три співвідношення. Тут все простіше, ніж ви могли очікувати з двох причин. По-перше, нагадаємо, що наш попередній розподіл ймовірностей - U (0, 1). Щільність цього розподілу є постійною (1,0) для всіх значень p і p ′. Через це коефіцієнт попереднього коефіцієнта для цього прикладу завжди:

\[ R_{prior} = \frac{P(p')}{P(p)} = \frac{1}{1} = 1 \label{2.28}\]

Аналогічно, оскільки розподіл наших пропозицій симетричний, Q (p ′| p) = Q (p | p ′) і R _{p r o p o s a l} = 1. Це означає, що нам потрібно лише обчислити коефіцієнт ймовірності, R _{l i k e l i h o o d} для p та p ′. Ми можемо зробити це, підключивши наші значення для p (або p ′) до Equation\ ref {2.2}:

\[ P(D|p) = {N \choose H} p^H (1-p)^{N-H} = {100 \choose 63} 0.6^63 (1-0.6)^{100-63} = 0.068 \label{2.29} \]

Так само,

\[ P(D|p') = {N \choose H} p'^H (1-p')^{N-H} = {100 \choose 63} 0.595^63 (1-0.595)^{100-63} = 0.064 \label{2.30}\]

Коефіцієнт ймовірності тоді:

\[ R_{likelihood} = \frac{P(D|p')}{P(D|p)} = \frac{0.064}{0.068} = 0.94 \label{2.31}\]

Тепер ми можемо обчислити R _{а с с е р т} = R _{р р i о р} ⋅ R р _{р р о р о с а л} ⋅ R _{л i k e l i h o o d} = 1 ⋅ 1 ⋅ 0.94 = 0.94. Далі ми вибираємо випадкове число між 0 і 1 - скажімо, що ми малюємо х = 0,34. Потім ми помічаємо, що наше випадкове число х менше або дорівнює R _{a c c e p t}, тому ми приймаємо запропоноване значення p ′. Якби випадкове число, яке ми намалювали, було більше 0,94, ми б відхилили запропоноване значення і зберегли наше початкове значення параметра p = 0.60, що йде в наступне покоління.

Якщо повторити цю процедуру велику кількість разів, то отримаємо довгий ланцюжок значень p. Результати такого пробігу ви можете побачити на малюнку 2.4. На панелі A я побудував ймовірності для кожного наступного значення p. Ви можете бачити, що ймовірність збільшення ймовірності для перших ~ 1000 або близько того поколінь, а потім досягають плато навколо l n L = −3. На панелі B наведено графік значень p, які швидко сходяться до стабільного розподілу навколо p = 0,63. Ми також можемо побудувати гістограму цих задніх оцінок р. У панелі C я зробив це - але з поворотом. Оскільки алгоритм MCMC створює ряд оцінок параметрів, ці цифри показують автокореляцію - тобто кожна оцінка схожа на оцінки, які надходять безпосередньо до і відразу після. Ця автокореляція може спричинити проблеми для аналізу даних. Найпростішим рішенням є підвибірка цих значень, вибираючи лише, скажімо, одне значення кожні 100 поколінь. Це те, що я зробив на гістограмі в панелі C. Ця панель також включає аналітичний задній розподіл, який ми розрахували вище - зверніть увагу, наскільки добре наш алгоритм Метрополіса-Гастінгса зробив у реконструкції цього розподілу! Для порівняльних методів загалом аналітичні задні розподіли важко або неможливо побудувати, тому наближення за допомогою MCMC є дуже поширеним явищем.

Малюнок 2.4. Байєсійський MCMC з ящірки перегортання приклад. Зображення автора, може бути використано повторно за ліцензією CC-BY-4.0.

Цей простий приклад охоплює деякі деталі алгоритмів MCMC, але ми розглянемо ці деталі пізніше, і є багато інших книг, які розглядають цю тему дуже глибоко (наприклад, Christensen et al. 2010). Справа в тому, що ми можемо вирішити деякі проблеми, пов'язані з байєсівською статистикою, використовуючи числові «трюки», такі як MCMC, які використовують потужність сучасних комп'ютерів, щоб відповідати моделям та оцінювати параметри моделі.

Розділ 2.4c: Фактори Байєса

Тепер, коли ми знаємо, як використовувати дані та до обчислення заднього розподілу, ми можемо перейти до теми вибору байєсової моделі. Ми вже вивчили один загальний метод вибору моделі з використанням АПК. Ми також можемо зробити вибір моделі в байєсівському рамках. Найпростіший спосіб - обчислити, а потім порівняти задні ймовірності для набору розглянутих моделей. Зробити це можна, розрахувавши коефіцієнти Байєса:

\[ B_{12} = \frac{Pr(D|H_1)}{Pr(D|H_2)} \label{2.32}\]

Фактори Байєса - це співвідношення граничних ймовірностей P (D | H) двох конкуруючих моделей. Вони являють собою ймовірність усереднених даних над заднім розподілом оцінок параметрів. Важливо зазначити, що ці граничні ймовірності відрізняються від ймовірностей, що використовуються вище для порівняння моделей A I C важливим чином. За допомогою A I C та інших пов'язаних тестів ми обчислюємо ймовірності для заданої моделі та певного набору значень параметрів — у прикладі перегортання монети ймовірність для моделі 2, коли p _H = 0,63. На відміну від цього, граничні ймовірності факторів Байєса дають ймовірність усереднених даних над усіма можливими значеннями параметрів для моделі, зваженими за їх попередньою ймовірністю.

Через використання граничних ймовірностей фактор Байєса дозволяє нам робити вибір моделі таким чином, що враховує невизначеність в наших оцінках параметрів - знову ж таки, за ціною вимагати явних попередніх ймовірностей для всіх параметрів моделі. Такі порівняння можуть сильно відрізнятися від тестів на коефіцієнт ймовірності або порівнянь балів A I C _c. Коефіцієнти Байєса представляють порівняння моделей, які інтегруються над усіма можливими значеннями параметрів, а не порівнюють придатність моделей лише за значеннями параметрів, які найкраще відповідають даним. Іншими словами, A I C _c оцінки порівнюють придатність двох моделей, заданих конкретними оціночними значеннями для всіх параметрів у кожній з моделей. На відміну від цього, фактори Байєса проводять порівняння між двома моделями, які враховують невизначеність у їх оцінках параметрів. Це зробить найбільшу різницю, коли деякі параметри однієї або обох моделей мають відносно широку невизначеність. Якщо всі параметри можна оцінити з точністю, результати обох підходів повинні бути схожими.

Розрахунок коефіцієнтів Байєса може бути досить складним, вимагаючи інтеграції між розподілами ймовірностей. У випадку нашої задачі перегортання монети ми вже зробили це для отримання бета-розподілу в Equation\ ref {2.22. Потім ми можемо обчислити коефіцієнти Байєса, щоб порівняти відповідність двох конкуруючих моделей. Давайте порівняємо дві моделі для перегортання монет, розглянуті вище: модель 1, де р _Н = 0,5, і модель 2, де р _Н = 0,63. Потім:

\[ \begin{array}{lcl} Pr(D|H_1) &=& \binom{100}{63} 0.5^{0.63} (1-0.5)^{100-63} \\\ &=& 0.00270 \\\ Pr(D|H_2) &=& \int_{p=0}^{1} \binom{100}{63} p^{63} (1-p)^{100-63} \\\ &=& \binom{100}{63} \beta (38,64) \\\ &=& 0.0099 \\\ B_{12} &=& \frac{0.0099}{0.00270} \\\ &=& 3.67 \\\ \end{array} \label{2.33}\]

У наведеному вище прикладі β (x, y) є бета-функцією. Наші розрахунки показують, що коефіцієнт Байєса становить 3,67 на користь моделі 2 порівняно з моделлю 1. Це, як правило, трактується як суттєві (але не вирішальні) докази на користь моделі 2. Знову ж таки, ми можемо бути досить впевнені, що наша ящірка не є справедливим ластом.

У прикладі перегортання ящірки ми можемо обчислити коефіцієнти Байєса саме тому, що ми знаємо розв'язку інтеграла в Equation\ ref {2.33. Однак, якщо ми не знаємо, як вирішити це рівняння (типова ситуація в порівняльних методах), ми все ще можемо наблизити фактори Байєса з наших пробігів MCMC. Методи для цього, включаючи вибірку зарозумілості та моделі сходинок (Xie et al. 2011; Perrakis et al. 2014), є складними і виходять за рамки цієї книги. Однак один із поширених методів наближення коефіцієнтів Байєса передбачає обчислення гармонійного середнього ймовірності за ланцюгом MCMC для кожної моделі. Співвідношення цих двох ймовірностей потім використовується як наближення коефіцієнта Байєса (Newton and Raftery 1994). На жаль, цей метод вкрай ненадійний, і, ймовірно, ніколи не повинен використовуватися (див. Ніл 2008 для більш докладної інформації).