Алгоритмы Витерби и Баума-Велши в скрытых марковских моделях

Алгоритм Витерби

Здравствуйте и вновь добро пожаловать на занятия по теме «Машинное обучение без учителя: скрытые марковские модели на языке Python».

Мы продолжаем разбирать скрытые марковские модели для дискретных наблюдений.

В этой лекции мы поговорим об алгоритме Витерби.

У вас может возникнуть вопрос касательно скрытых марковских моделей, особенно если ваши скрытые состояния смоделированы на основе некоторых реальных физических процессов или действующих систем, а не просто произвольно взятых скрытых переменных, – какова последовательность скрытых состояний с учётом имеющихся наблюдений? К примеру, пусть дана последовательность звуков чьей-то речи, а нам нужно определить звучащие слова.

Это именно то, что вычисляет алгоритм Витерби, – наиболее вероятную последовательность скрытых состояний с учётом наблюдающейся последовательности в модели. Как вы увидите, алгоритм Витерби во многом похож на только что обсуждавшийся алгоритм прямого хода, с тем исключением, что вместо нахождения суммы используется нахождение максимума. Кроме того, придётся добавить возврат в прошлые моменты времени, поскольку нам необходимо отслеживать фактические состояния, а не просто определить конечную вероятность.

Для этого нам придётся ввести две новые переменные. Во-первых, это δ(t,i) с индексами времени и состояния, которая представляет максимальную вероятность окончания пребывания в состоянии i в момент t. Это совместное распределение вероятностей по последовательности состояний и наблюдаемой последовательности:

Во-вторых, мы также определим ψ(t,i), также с индексами времени и состояния, чтобы отслеживать фактические последовательности состояний, заканчивающихся в момент t в состоянии i.

И вновь у нас три этапа. Этап инициации очень похож на аналог из алгоритма прямого хода:

Далее идёт этап рекурсии, на протяжении которого вычисляются значения δ и ψ для всех моментов времени и всех состояний:

И, наконец, завершающий этап, когда мы находим максимальную вероятность для лучшей последовательности состояний и лучшего конечного состояния:

Чтобы определить остальную часть лучшей последовательности, мы возвращаемся назад и используем ψ, повторяя эту процедуру по нисходящей от момента T – 1 до 1:

Если бы мы писали код, он выглядел бы примерно так:

        delta = np.zeros((T, self.M))

        psi = np.zeros((T, self.M))

        delta[0] = self.pi*self.B[:,x[0]]

        for t in xrange(1, T):

            for j in xrange(self.M):

                delta[t,j] = np.max(delta[t-1]*self.A[:,j]) * self.B[j, x[t]]

                psi[t,j] = np.argmax(delta[t-1]*self.A[:,j])

Как вы можете видеть, δ и ψ являются массивами размерностью TxM. Мы инициируем δ[0] таким же образом, как делали это для α, после чего используем цикл для прохода от момента времени 1 до момента T, а также цикл для прохода по всем состояниям с одновременным обновлением δ и ψ. Это точная реализация предыдущих формул.

Следующее – возврат назад для получения всех состояний:

        # backtrack

        states = np.zeros(T, dtype=np.int32)

        states[T-1] = np.argmax(delta[T-1])

        for t in xrange(T-2, -1, -1):

            states[t] = psi[t+1, states[t+1]]

        return states

Мы определяем пустой массив размерностью T и названием states, инициируем последнее состояние как функцию argmax от последнего значения δ и используем обратный цикл про ψ.

Наглядное представление алгоритма Витерби

Основное отличие состоит в том, что вместо нахождения суммы по всем предыдущим состояниям, которые могли бы привести к состоянию текущему, мы находим максимум. Это связано с тем, что нам нужно отыскать последовательность, которая даст нам наибольшую вероятность наблюдения.

 

Итак, вначале рассмотрим только одно наблюдение x(1). Поскольку одно наблюдение соответствует лишь одному состоянию, то нам нужно найти наиболее вероятное единственное состояние. Оно равно максимуму произведения π и B, являющемся вероятностью перехода к данному состоянию и наблюдению того, что наблюдается из этого состояния:

Не забывайте, что мы определяем произведение π и B как δ с индексами времени и состояния, поэтому

Как и в случае алгоритма прямого хода, теперь представим себе двушаговую последовательность (на рисунке x(1) не показан для компактности). Предположим, мы рассматриваем только состояние 1. Как в него попасть? Вообще говоря, мы можем попасть в него из любого другого состояния. Таким образом, основной посыл таков: нужно вновь выбрать то, что даёт максимальную вероятность.

Как же получить эту вероятность? Это предыдущее значение δ (которое, если вы помните, является произведением π и B в момент t = 1) и последующий переход из того состояния в текущее, которое мы принимаем равным 1, а затем умножаем на вероятность наблюдения из этого состояния:

Это и есть определение следующего значения δ. Таким образом, δ также выражается рекурсивно через свои предыдущие значения.

В итоге мы получаем лучшую последовательность состояний, учитывающую все наблюдаемые переменные в целом. Это значит, что когда мы рассматриваем только первое наблюдение – предположим, пребывая в состоянии 1 – мы получаем просто лучшую общую вероятность. Но если мы рассматриваем уже два наблюдения, включая пребывание в состоянии 2, а потом переход в состояние 1, – мы получаем лучшую общую вероятность для всей последовательности. Поэтому мы должны указать, что состояние в момент t = 1 – это состояние 2.

Поскольку мы можем определить лучшую последовательность только при самом последнем значении δ, нам необходимо отслеживать все возможные переходы состояний, с которыми сталкивались. Тут нет экспоненциального роста, поскольку нет нужды учитывать все возможные переходы состояний в целом, а лишь те из них, которые дают наибольшую вероятность на всём пути.

Поскольку δ(T, i) является наибольшей вероятностью для наблюдения всей последовательности с окончанием в состоянии i, то её максимум даёт нам наилучшую вероятность в целом. Кроме того, функция argmax даёт нам лучшее последнее состояния. А поскольку мы отслеживали все переходы, мы знаем, как попали в последнее состояние из предпоследнего, как попали в предпоследнее из предпредпоследнего и так далее.

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

Настало время поговорим о последней и самой важной из трёх задач, решаемых в рамках скрытых марковских моделей, – о том, как их обучать.

Как и в случае, когда мы изучали гауссовы смеси распределений, у нас есть скрытые переменные, для которых надо суммировать все возможности, и мы не можем просто найти максимум правдоподобия. Нам придётся использовать итерационный алгоритм максимизации ожидания, названный алгоритмом Баума-Велша в честь математиков, которые его изобрели. Основывается он на алгоритме прямого-обратного хода, о котором мы узнали ранее.

Действительно, первым этапом алгоритма Баума-Велша является вычисление переменных прямого и обратного хода – α и β.

После этого вводится новая величина φ с тремя индексами t, i и j, которая определяется как вероятность пребывания в состоянии i в момент t с переходом в состояние j в момент t + 1 при данной последовательности наблюдений:

Диаграмма на слайде, надеюсь, даёт наглядное представление об этом процессе.

После φ вводится ещё одна переменная γ, зависящая только от t и i и потому суммирующаяся по j, что представляет собой все состояния: Таким образом, γ – это всего лишь вероятность φ с исключённым j:

Суть в том, что когда мы суммируем эти переменные по всем моментам времени, то γ будет представлять ожидаемое количество переходов из состояния i, а φ – ожидаемое количество переходов из состояния i в состояние j:

Это уже очень похоже на процесс обучения марковских моделей. В связи с этим можно переписать наши обновления следующим образом:

При этом следует иметь в виду, что сейчас рассматривается лишь одна последовательность. Но, разумеется, мы хотим приспособить нашу модель для всех обучающих последовательностей. Как это сделать? Об этом – в следующей лекции.

Популярное изложение алгоритма Баума-Велша

Мы должны более подробно обсудить алгоритм Баума-Велша, поскольку, как мне сообщали, понимание его действия может вызывать затруднения. Так что давайте подробнее разберём эту тему, чтобы её было легче усвоить.

Чтобы полностью понять алгоритм Баума-Велша, необходимо разобраться в нескольких его частях и этапах. Во-первых, необходимо понять уравнения алгоритма на чисто механическом уровне – все эти произведения, суммирования и деления. Мы вводим некоторые новые переменные, представляющие собой некоторые вероятности. Если вам непонятно, почему они представляют эти вероятности, то следует вернуться назад и вновь изучить определения α и β – пока они не станут понятными. Опять же, это просто следование правилам теории вероятностей. Ну, и не следует упускать из виду теорему Байеса.

Во-вторых, вы должны уметь воплотить эти уравнения в коде. Вновь-таки, это просто умножение и деление, которые мы уже знаем, как реализовать в Numpy. Вы можете вычислять суммы, используя цикл for, но куда лучше частично, а если это возможно, то и полностью, преобразовать обновления в векторную форму. Простая практика в работе с уравнениями поможет вам лучше их понимать.

И последний этап, самый сложный – понять, как обновляется алгоритм Баума-Велша. Его основополагающие принципы выходят за рамки данного курса, но основываются они на так называемом алгоритме максимизации ожиданий. Если вы хотите ознакомиться ближе с алгоритмом максимизации ожиданий, то вам следует обратиться к курсу по байесовому моделированию или графическим моделям. В нашем же случае лучший способ получить интуитивно понятное представление о том, как он работает, – это вернуться к кластеризации методом k-средних и гауссовым смесям распределений. Основная же идея в том, что это итерационная процедура, покоящаяся на двух этапах – этапе ожиданий и этапе максимизации.

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

Как легко видеть, в связи со структурой английского языка существует марковская вероятностная модель, лежащая в основе последовательности разметки по частям речи. При прогнозировании у нас будет новое предложение без разметки по частям речи. Мы должны отыскать наиболее вероятную последовательность скрытых состояний, которая и будет прогнозируемой последовательностью разметки по частям речи. Какой тут нужен алгоритм? Это обсуждавшийся ранее алгоритм Витерби, дающий наиболее вероятную последовательность скрытых состояний.

Интересным в этой задаче является то, что в обучающем наборе данных скрытые состояния на самом деле не являются скрытыми, мы точно знаем, что они из себя представляют. Вследствие этого мы можем использовать функцию максимума правдоподобия в аналитическом виде для определения всех переходов по скрытым состояниям. Так, например, вероятность перехода от существительного к глаголу – это количество переходов от существительного к глаголу, поделённая на количество существительных:

Аналогичным образом можно найти вероятность наблюдения. К примеру, слово «печь» может быть как существительным, так и глаголом. Легко найти и вероятность слова «печь» как существительного, и его вероятность как глагола, поскольку в обучающем наборе у нас уже есть все эти данные:

Таким образом, всё сводится к простому подсчёту.

В общем случае проблема скрытых марковских моделей состоит в том, что скрытые состояние неизвестны. Что же делать в таком случае?

Вернёмся на минуту к гауссовым смесям распределений. Напомним, что в них обновление происходит в два этапа: сначала находится лучшие пропорции смесей при данных текущих средних значениях и дисперсии, а затем находятся лучшие средние значения и дисперсия при данных текущих пропорциях смесей. Причина, по которой этой делается, заключается в том, что мы не знаем, от какой гауссианы пришли конкретные данные. Поэтому-то мы и определяем их с точки зрения вероятностей, являющихся пропорциями смесей, – какой из гауссиан является скрытой переменной в гауссовых смесях распределений. Если бы мы знали, какая гауссиана отвечает за данные, нам не потребовалось бы обращаться к максимизации ожиданий.

В скрытых марковских моделях идея та же. Мы не знаем, какие скрытые состояния являются истинными, потому определяем возможные состояние, вероятностным способом – это этап вычисления ожиданий:

Этап максимизации же состоит в нахождении лучших значений π, A и B с точки зрения этих вероятностно определённых скрытых состояний.

Другой способ представления обновлений алгоритма Баума-Велша – с точки зрения множителей Лагранжа. Метод множителей Лагранжа сам по себе довольно сложен, так что пригодится он только если вы уже знакомы с ним. Напомним, что этот метод работает при решении задачи оптимизации некоторой целевой функции при заданных ограничениях. В данном случае нам нужно оптимизировать вероятность наблюдения последовательности при некоторых заданных параметрах π, A и B:

 P=P(X|\theta), \theta = \{\pi,A,B\}.

Ограничениями тут являются требования, чтобы сумма всех π равнялась единице, а также чтобы суммы всех строк A и B тоже равнялись единице:

Если вы знаете метод множителей Лагранжа достаточно хорошо, чтобы решить данную задачу, то в итоге у вас получится вот такое решение:

Решив эти уравнения для производной, вы в конце концов получите то же решение, которое мы уже видели.

Резюмируем: независимо от используемого метода, он должен быть довольно изощрённым. Вам понадобится значительно больше знаний, чем предлагается в данном курсе. Но лучший способ достичь понимания – это сначала усвоить механику, реализацию и интуитивное представление.

Обновления алгоритма Баума-Велша для множественных наблюдений

Теперь мы можем узнать, как обучать скрытую марковскую модель в случае, если у нас есть несколько последовательностей наблюдений.

Вы легко можете представить ситуацию, когда наши обучающие данные состоят не из одного лишь конкретного случая, а нам нужно, чтобы модель точно соответствовала всем обучающим данным в целом.

Вначале подумаем о том, как можно представить эти данные. Обычно когда у нас есть N примеров размерностью D, мы можем поместить их в массив, или матрицу, Numpy размерностью NxD. Но когда речь заходит о последовательностях, возникает проблема. Почему так происходит? Предположим, мы рассматриваем примеры голосовых сообщений. Тогда чьё-то сообщение «Привет, мир» будет гораздо короче, чем сообщение «Я сэкономил 15% на автомобильной страховке, перейдя к компании Geico».

Простейший способ – просто сохранять каждый конкретный пример в качестве элемента в списке Python. Внутри такого списка у нас могут быть отдельные массивы Numpy любой длины. Мы можем их обозначить T(n), где n = 1..N. При этом каждая из N последовательностей будет иметь собственные значения всех обновляемых переменных, о которых мы говорили ранее. То есть, у нас будет P(n), представляющую вероятность последовательности, αn – переменную прямого хода, и βn, представляющую переменную обратного хода. Сейчас нам удобнее выражать обновления через α и β, поэтому пока что забудем о φ и γ.

Обновление для π будет равно

В качестве упражнения попытайтесь доказать, что их сумма будет равна 1.

Обновление для A(i,j) представляет собой весьма большое выражение, но включает в себя только α, A, B и β, то есть все уже обсуждавшиеся переменные:

Обновление для B(j,k) равны следующему:

Что интересно в этих уравнениях, так это то, что каждая из внутренних сумм умножается на величину, обратную вероятности наблюдения. О чём это говорит? Это значит, что если в рамках данной модели наблюдение, которое мы желаем промоделировать, является очень маловероятным, то обновления, основанные на этом наблюдении, будут более значимыми.

Итак, как это всё может выглядеть в коде?

Для π всё очень просто. Мы лишь берём первый элемент α, первый элемент β, суммируем по всем n и делим на N:

self.pi = np.sum((alphas[n][0] * betas[n][0])/P[n] for n in xrange(N)) / N

Для обновления A мы инициируем матрицу размерностью MxM, представляющую числитель, организовываем цикл по всем n образцам, внутри которого проходим по i, равному всем состояниям, j, равным всем состояниям, и t, равным всем моментам времени, кроме последнего:

       a_num = np.zeros((self.M, self.M))

       for n in xrange(N):

                x = X[n]

                T = len(x)

     den1 += (alphas[n][:-1] * betas[n][:-1]).sum(axis=0, keepdims=True).T / P[n]

                a_num_n = np.zeros((self.M, self.M))

                for i in xrange(self.M):

                    for j in xrange(self.M):

                        for t in xrange(T-1):

        a_num_n[i,j] += alphas[n][t,i] * betas[n][t+1,j] * self.A[i,j] * self.B[j, x[t+1]]

           a_num += a_num_n / P[n]

         self.A = a_num / den1

В качестве упражнения можете попробовать полностью преобразовать это уравнение в векторную форму. Обратите внимание, что знаменатель уже преобразован, но в качестве ещё одного упражнения попробуйте перевести его в исходную форму и посмотреть, всё ли правильно.

Следующее – обновления для B. Для этого мы для числителя инициируем матрицу размерностью MxV и суммируем в числителе, только если x(t) = j:

            b_num = np.zeros((self.M, V))

            for n in xrange(N):

                x = X[n]

                T = len(x)

                den2 += (alphas[n] * betas[n]).sum(axis=0, keepdims=True).T / P[n]

                b_num_n = np.zeros((self.M, V))

                for i in xrange(self.M):

                    for j in xrange(V):

                        for t in xrange(T):

                            if x[t] == j:

                                b_num_n[i,j] += alphas[n][t][i] * betas[n][t][i]

                b_num += b_num_n / P[n]

            self.B = b_num / den2

Вновь-таки, можете в качестве упражнения преобразовать знаменатель из векторной формы, чтобы посмотреть, всё ли правильно.  Можно также попробовать преобразовать в векторную форму числитель.

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

Итак, вначале мы просто определим переменные скрытой марковской модели:

symbol_map = [‘H’, ‘T’]

pi = np.array([0.5, 0.5])

A = np.array([[0.1, 0.9], [0.8, 0.2]])

B = np.array([[0.6, 0.4], [0.3, 0.7]])

M, V = B.shape

А затем создадим функцию для генерации последовательности длиной N:

def generate_sequence(N):

    s = np.random.choice(xrange(M), p=pi) # начальное состояние

    x = np.random.choice(xrange(V), p=B[s]) # начальное наблюдение

    sequence = [x]

    for n in xrange(N-1):

        s = np.random.choice(xrange(M), p=A[s]) # следующее состояние

        x = np.random.choice(xrange(V), p=B[s]) # следующее наблюдение

        sequence.append(x)

    return sequence

Как и в обычных марковских моделях, вначале мы случайным образом выбираем состояние, а выбрав его, случайным образом выбираем наблюдение. Заметьте, что это в точности то, что предполагается при моделировании физической системы или чего-то подобного, когда данные поступают извне.

Файл данных называется coin_data.txt и представляет собой прямоугольник размером 50х30, заполненный символами H и T.

Андрей Никитенко
Андрей Никитенко
Задать вопрос эксперту
Понравилась статья? Поделить с друзьями:
Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!:

Share via
Copy link