8.7: Додаток - Алгоритм обрізки Фельзенштейна
- Page ID
- 4702
Алгоритм обрізки Фельзенштейна (1973) є прикладом динамічного програмування, типу алгоритму, який має багато застосувань у порівняльній біології. У динамічному програмуванні ми розбиваємо складну задачу на ряд простих кроків, які мають вкладену структуру. Це дозволяє нам ефективно використовувати обчислення та прискорює час, необхідний для проведення розрахунків.
Найкращий спосіб проілюструвати алгоритм Фельзенштейна - на прикладі, який представлений на панелах нижче. Ми намагаємося обчислити ймовірність тридержавного характеру на філогенетичному дереві, що включає шість видів.
Малюнок 8.4A. Кожен наконечник і внутрішній вузол у дереві мають три поля, які міститимуть ймовірності для трьох станів символів у цій точці дерева. Перше поле представляє стан 0, другий стан 1 і третій стан 2. Зображення автора, може бути використано повторно за ліцензією CC-BY-4.0.
- Першим кроком в алгоритмі є заповнення ймовірностей для підказок. В даному випадку нам відомі стани на кінчиках дерева. Математично ми стверджуємо, що ми точно знаємо стан характеру на кінчиках; ймовірність того, що цей вид має стан, який ми спостерігаємо, становить 1, а всі інші стани мають нуль ймовірності:
- Далі ми визначаємо вузол, де всі його безпосередні нащадки є підказками. Завжди буде хоча б один такий вузол; часто, їх буде більше одного, в такому випадку ми довільно виберемо один. Для цього прикладу ми виберемо вузол, який є останнім загальним предком видів A і B, позначений як вузол 1 на малюнку 8.2B.
- Потім ми використовуємо рівняння 7.6 для обчислення умовної ймовірності для кожного стану символів для піддерева, що включає вершину, яку ми вибрали на кроці 2, та його нащадків. Для кожного стану символу умовною ймовірністю є ймовірність, враховуючи дані та модель, отримання станів символу підказки, якщо ви починаєте з цього стану символу в корені. Іншими словами, ми відстежуємо ймовірність для верхніх частин дерева, включаючи наші дані, якщо розглянутий нами вузол мав кожен з можливих станів символів. Цей розрахунок є:
\[ L_P(i) = (\sum\limits_{x \in k}Pr(x|i,t_L)L_L(x)) \cdot (\sum\limits_{x \in k}Pr(x|i,t_R)L_R(x)) \label{8.7}\]
де i та x - це обидва індекси для k символьних станів, з сумами, взятими у всіх можливих станах на кінчиках гілок (x), та термінами, розрахованими для кожного можливого стану у вузлі (i). Два шматки рівняння є лівим і правим нащадком цікавить вузла. Гілки можуть бути призначені як ліві або праві довільно, не впливаючи на кінцевий результат, і підхід також працює для політомій (але рівняння трохи інше). Крім того, кожен з цих двох частин сам має дві частини: ймовірність початку і закінчення кожного стану вздовж двох розглянутих гілок, і поточні умовні ймовірності, які входять до рівняння на кінчиках піддерева (L (x) і Л Р (х)). Довжини гілок позначаються як t L і t R для лівого і правого відповідно.
Можна думати про ймовірність «стікання» по гілках дерева, а умовні ймовірності для лівої та правої гілок об'єднуються шляхом множення на кожному вузлі, генеруючи умовну ймовірність для батьківського вузла для кожного символьного стану (L P ( i)).
Розглянемо піддерево, що веде до видів A і B у наведеному прикладі. Два стани характеру наконечника 0 (для виду A) і 1 (для видів B). Ми можемо обчислити умовну ймовірність для стану символу 0 у вузлі 1 як:
\[ L_P(0) = \left(\sum\limits_{x \in k}Pr(x|0,t_L=1.0)L_L(x)\right) \cdot \left(\sum\limits_{x \in k}Pr(x|0,t_R=1.0)L_R(x)\right) \label{8.8}\]
Далі ми можемо обчислити частки ймовірності з матриці ймовірностей P. У цьому випадку t L = t R = 1.0, тому і для лівої, і для правої гілки:
\[ \mathbf{Q}t = \begin{bmatrix} -2 & 1 & 1 \\ 1 & -2 & 1 \\ 1 & 1 & -2 \\ \end{bmatrix} \cdot 1.0 = \begin{bmatrix} -2 & 1 & 1 \\ 1 & -2 & 1 \\ 1 & 1 & -2 \\ \end{bmatrix} \label{8.9}\]
Так що:
\[ \mathbf{P} = e^{Qt} = \begin{bmatrix} 0.37 & 0.32 & 0.32 \\ 0.32 & 0.37 & 0.32 \\ 0.32 & 0.32 & 0.37 \\ \end{bmatrix} \label{8.10}\]
Тепер зверніть увагу, що, оскільки стан лівого символу, як відомо, дорівнює нулю, L (0) =1 і L L (1) = L L (2) =0. Аналогічно правий стан одне, тому L R (1) =1 і L R (0) = L R (2) =0.
Тепер ми можемо заповнити дві частини рівняння 8.2:
\[ \sum\limits_{x \in k}Pr(x|0,t_L=1.0)L_L(x) = 0.37 \cdot 1 + 0.32 \cdot 0 + 0.32 \cdot 0 = 0.37 \label{8.11} \]
і:
\[ \sum\limits_{x \in k}Pr(x|0,t_R=1.0)L_R(x) = 0.37 \cdot 0 + 0.32 \cdot 1 + 0.32 \cdot 0 = 0.32 \label{8.11B} \]
Отже:
Л Р (0) = 0,37 ⋅ 0,32 = 0,12.
Це означає, що за моделлю, якби стан у вузлі 1 дорівнював 0, ми мали б ймовірність 0,12 для цієї невеликої ділянки дерева. Ми можемо використовувати подібний підхід, щоб виявити, що:
(ур. 8,13)
Л Р (1) = 0,32 ⋅ 0,37 = 0,12.
Л Р (2) = 0,32 ⋅ 0,32 = 0,10.
Тепер у нас є ймовірність для всіх трьох можливих родових держав. Ці цифри можна ввести у відповідні графи:
- Потім ми повторюємо вищевказаний розрахунок для кожного вузла в дереві. Для вузлів 3-5 не всі члени L L (x) і L R (x) дорівнюють нулю; їх значення можна прочитати з квадратів на дереві. Результат всіх цих розрахунків:
- Тепер ми можемо обчислити ймовірність по всьому дереву, використовуючи умовні ймовірності для трьох станів в корені дерева.
\[ L = \sum\limits_{x \in k} \pi_x L_{root} (x) \label{8.14}\]
Де π x - попередня ймовірність цього символьного стану в корені дерева. Для цього прикладу ми візьмемо ці попередні ймовірності рівномірними, рівними для кожного стану (π х = 1/ k = 1/3). Імовірність для нашого прикладу, значить, така:
\[L = 1/3 ⋅ 0.00150 + 1/3 ⋅ 0.00151 + 1/3 ⋅ 0.00150 = 0.00150 \label{8.15}\]
Зауважте, що якщо ви спробуєте цей приклад в іншому програмному пакеті, наприклад GEIGER або PAUP*, програма обчислить ln-ймовірність -6.5, що є саме природним журналом значення, розрахованого тут.