Алгоритм Баума-Велша

Алгоритм Баума-Велша (англ. Baum–Welch algorithm) — алгоритм для нахождения неизвестных параметров скрытой Марковской модели. Использует алгоритм прямого-обратного хода.

История

Скрытые Марковские модели (HMMs) и алгоритм Баума-Велша впервые были описаны в заметках Леонарда Баума и его сверстников в конце 1960х. Одно из первых основных приложений на основе HMMs было использовано в области обработки речи. В 1980х HMMs стало эффективным инструментом в анализе биологических систем и информации, особенно в генном анализе.

Описание алгоритма

Пусть [math]Q_t[/math] — это дискретная случайная переменная, принимающая одно из [math]N[/math] значений [math](1..N)[/math]. Будем полагать, что данная модель Маркова, определенная как [math]P(Q_t | Q_{t - 1})[/math] однородна по времени, то есть независима от [math]t[/math]. Тогда можно задать [math]P(Q_t | Q_{t - 1}) [/math] как независящую от времени стохастическую матрицу перемещений [math]A = \{a_{ij}\} = p(Q_t = j | Q_{t - 1} = i)[/math]. Особый случай для времени [math]t = 1[/math] определяется начальным распределением [math]\pi_i = P(Q_1 = i)[/math].

Будем считать, что мы в состоянии [math]j[/math] в момент времени [math]t[/math], если [math]Q_t = j[/math]. Последовательность заданных состояний определяется как [math]q = (q_1, ..., q_T)[/math], где [math]q_t \in \{ 1..N\}[/math] является состоянием в момент времени [math]t[/math].

Наблюдение может иметь одно из [math]L[/math] возможных значений, [math]Q_t \in \{o_1, ..., o_L\}[/math]. Вероятность заданного вектора наблюдений в момент времени [math]t[/math] для состояния [math]j[/math] определяется как [math]b_j(o_t) = P(O_t = o_t | Q_t = j)( B = \{ b_{ij}\}[/math] — это матрица [math]L[/math] на [math]N)[/math]. Заданная последовательность наблюдений [math]O[/math] выражается как [math] O = (O_1 = o_1, ..., O_T = o_T)[/math].

Следовательно, мы можем описать скрытую модель Маркова с помощью [math] \lambda = (A, B, \pi)[/math]. При заданном векторе наблюдений [math]O[/math] алгоритм Баума-Велша находит [math] \lambda^*=\max_\lambda P(O\mid\lambda)[/math]. [math]\lambda[/math] максимизирует вероятность наблюдений [math]O[/math].

Исходные данные: [math] \lambda = (A, B, \pi)[/math] со случайными начальными условиями. Алгоритм итеративно обновляет параметр [math]\lambda[/math] до схождения в одной точке.

Прямая процедура

[math]a_i(t) = p(O_1 = o_1, ..., O_t = o_t, Q_t = _i | \lambda[/math], что является вероятностью получения заданной последовательности [math]\{ o_1, ..., o_t\}[/math] для состояния [math]i[/math] в момент времени [math]t[/math].

[math]a_i(t)[/math] можно вычислить рекурсивно:

1.[math]a_i(1) = \pi_i \cdot b_i(O_1) [/math];

2.[math]a_j(t + 1) = b_j(O_{t + 1})\displaystyle\sum^N_{i=1}a_i(t) \cdot a_{ij}[/math].

Обратная процедура

Данная процедура позволяет вычислить вероятность конечной заданной последовательности [math]o_{t + 1}, ..., o_T[/math] при условии, что мы начали из исходного состояния [math]i[/math], в момент времени [math]t[/math].

[math]\beta_i(t)[/math] можно вычислить рекурсивно:

1.[math]\beta_i(T) = 1[/math];

2. [math]\beta_i(t) = \displaystyle\sum^N_{j = 1}\beta_j(t + 1)a_{ij}b_j(O_{t + 1})[/math].

Обновление переменных

Определим временные переменные:

[math]\gamma_i(t) = p(Q_t=i\mid O,\;\lambda)=[/math] [math]\frac{\alpha_i(t)\beta_i(t)}{\displaystyle\sum^N_{j=1}\alpha_j(t)\beta_j(t)}[/math]

[math]\xi_{ij}(t) = p(Q_t=i,\;Q_{t+1}=j\mid O,\;\lambda)=[/math] [math]\frac{\alpha_i(t)a_{ij}\beta_j(t+1)b_j(o_{t+1})}{\displaystyle\sum^N_{i=1}\displaystyle\sum^N_{j=1}\alpha_i(t)a_{ij}\beta_j(t+1)b_j(O_{t+1})}[/math].

Имея [math]\gamma[/math] и [math]\xi[/math], можно определить:

[math]\bar\pi_i=\gamma_i(1)[/math],

[math]\bar{a}_{ij}=\frac{\displaystyle\sum^{T-1}_{t=1}\xi_{ij}(t)}{\displaystyle\sum^{T-1}_{t=1}\gamma_i(t)}[/math],

[math]\bar{b}_i(k)=\frac{\displaystyle\sum^T_{t=1}\delta_{O_t,\;o_k}\gamma_i(t)}{\displaystyle\sum^T_{t=1}\gamma_i(t)}[/math].

Используя новые переменные [math] A, B, \pi[/math] итерации продолжаются до схождения.

Пример

Предположим, у нас есть курица, с которой мы собираем яйца. Снесла ли курица яйца — зависит от некоторых неизвестных факторов. Для простоты предположим, что существуют лишь два состояния, которые определяют куриные ли это яйца. В начальный момент нам неизвестно текущее состояние, также нам неизвестна вероятность перехода из одного состояния в другое. Для начала возьмем произвольные матрицы переходов и состояний.

Переходы Состояние 1 Состояние 2 Состояние 1 0.5 0.5 Состояние 2 0.3 0.7

Состояния Яйца не отложены Яйца отложены Состояние 1 0.3 0.7 Состояние 2 0.8 0.2

Начальное состояние Состояние 1 0.2 Состояние 2 0.8

Рассмотрим набор наблюдений (E — яйца отложены, N — яйца не отложены): NN, NN, NN, NN, NE, EE, EN, NN, NN.

Следующим шагом оценим новую матрицу переходов:

Последовательность	Вероятность последовательности и состояний	Наибольшая вероятность наблюдения
NN	0.024	0.3584 S2,S2
NN	0.024	0.3584 S2,S2
NN	0.024	0.3584 S2,S2
NN	0.024	0.3584 S2,S2
NE	0.006	0.1344 S2,S1
EE	0.014	0.0490 S1,S1
EN	0.056	0.0896 S2,S2
NN	0.024	0.3584 S2,S2
NN	0.024	0.3584 S2,S2
Итог	0.22	2.4234

Таким образом получаем новую оценку перехода из [math]S1[/math] в [math]S2[/math], которая составляет [math]\frac{0.22}{2.4234}[/math][math] = 0.0908[/math]. После этого можно подсчитать вероятность переходов из [math]S2[/math] в [math]S1[/math], [math]S2[/math] в [math]S2[/math], [math]S1[/math] в [math]S1[/math] и изменим их так, чтобы в суммы вероятностей давали 1. В итоге получаем новую матрицу переходов:

Старая матрица Состояние 1 Состояние 2 Состояние 1 0.5 0.5 Состояние 2 0.3 0.7

Новая матрица (Псевдовероятности) Состояние 1 Состояние 2 Состояние 1 0.0598 0.0908 Состояние 2 0.2179 0.9705

Новая матрица (После изменения) Состояние 1 Состояние 2 Состояние 1 0.3973 0.6027 Состояние 2 0.1833 0.8167

Далее оценим новую матрицу состояний:

Последовательности	Наибольшая вероятность наблюдения Если допустимо, что E получено из [math]S1[/math]	Наибольшая вероятность наблюдения
NE	0.1344 S2,S1	0.1344 S2,S1
EE	0.0490 S1,S1	0.0490 S1,S1
EN	0.0560 S1,S2	0.0896 S1,S2
Итог	0.2394	0.2730

Новая оценка для E, полученная из [math]S1[/math], составляет [math]\frac{0.2394}{0.2730}[/math] [math] = 0.8769[/math].

Благодаря этому, возможно рассчитать матрицу состояний:

Старая матрица Яйца не отложены Яйца отложены Состояние 1 0.3 0.7 Состояние 2 0.8 0.2

Новая матрица (Оценка) Яйца не отложены Яйца отложены Состояние 1 0.0876 0.8769 Состояние 2 1.0000 0.7385

Новая матрица (После изменения) Яйца не отложены Яйца отложены Состояние 1 0.0908 0.9092 Состояние 2 0.5752 0.4248

Для оценки начальной вероятности, мы предполагаем, что все последовательности начаты со скрытого состояния [math]S1[/math] и рассчитаны с высокой вероятностью, а затем повторяем для [math]S2[/math]. После нормализации получаем обновленный исходный вектор.

Повторяем эти шаги до тех пор, пока вероятности не сойдутся.

Псевдокод

 // T — конечный момент времени
 int[] DynamicOptionalStateSequance([math]\lambda[/math], d):
     double [math]\gamma[/math][1, i] = [math]\pi[/math][i] * b[i, d[1]]
     int [math]\psi[/math][1, i] = []
     int ans[]
     for t = 2 to T
         for i = 1 to n
             if [math]\gamma[/math][t, j] < [math]\gamma[/math][t - 1, i] * a[i, j] * b[j, d[t]]
                 [math]\gamma[/math][t, j] = [math]\gamma[/math][t - 1, i] * a[i, j] * b[j, d[t]]
                 [math]\psi[/math][t, j] = i
     ans[T] = 1 
     for i = 2 to n
         if [math]\gamma[/math][T, i] > [math]\gamma[/math][T, i - 1]
             ans[T] = i    
     for t = T - 1 downto 1
         ans[t] = [math]\psi[/math][t + 1, ans[t + 1]]
 return ans

Источники информации

• Википедия — Алгоритм Баума-Велша

• Лекция "Скрытые Марковские Модели" Сергея Николенко

• Wikipedia — Baum–Welch algorithm