Скрытые марковские модели для непрерывных наблюдений

Гауссовы смеси распределений и скрытые марковские модели

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

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

Рассмотрим аудиосигнал. Аудиосигнал – это просто звук, а звук – это вибрирующий поток воздуха, причем под вибрацией подразумевается, что происходит процесс колебаний во времени. Аудио может храниться на компьютере. Очевидно, что компьютер не может хранить бесконечное число значений, а потому поток воздуха нужно квантовать, то есть представить в виде чисел, которые компьютер может воспроизвести. Квантованию подвергаются две вещи: амплитуда звуковых вибраций и момент, когда эта амплитуда возникла. Таким образом, и время само по себе нужно квантовать. Как указывалось ранее, звук обычно дискретизируется с частотой 44,1 кГц, что, разумеется, слишком много для нашего мозга, чтобы понять, что он не непрерывен.

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

При использовании скрытых марковских моделей и гауссовых смесей распределения появляются следующие изменения. Как вы помните, скрытая марковская модель определяется тремя параметрами: π, A и B. Именно B должно поменяться, поскольку раньше это была матрица размерностью MxV, но у нас больше нет V возможных символов. Напомним, что для гауссовых смесей распределений также нужны три параметра. Нам нужны соответствия, показывающие вероятность конкретной гауссианы, а затем для каждой отдельной гауссианы нам нужны среднее значение и ковариация. В связи с этим мы вместо B подставляем три новых параметра – R, μ и Σ. Кроме того, буквой K мы обозначаем количество гауссиан. Если вы хотите посмотреть, как работают гауссовы смеси распределений и каковы их параметры, то материалы по этому вопросу находятся в приложении.

R представляется в виде матрицы размерностью MxK; таким образом, для каждого состояния выделяется по одной строке, и у каждого состояния есть K различных вероятностей. Поскольку это матрица вероятностей, то сумма значений в каждой строке должна равняться 1.

Каждое конкретное μ имеет размерность D, а всего их K штук; кроме того, такое множество μ нужно для каждого из M состояний. Поэтому μ будет матрицей размерностью MxKxD – трёхмерной матрицей.

Аналогично Σ будет иметь размерность MxKxDxD – то есть будет четырёхмерной матрицей.

Выбрав некоторое состояние j, мы получим вероятность наблюдения, равную

где N – функция распределения вероятностей для многомерной гауссианы.

Хотя формально у нас больше нет матрицы B, мы некоторым образом можем её создать, вычисляя отдельно для каждой последовательности.

Вычислим также отдельные компоненты, не проводя суммирование:

Теперь мы можем вычислить новое значение γ, являющимся частью первого этапа алгоритма максимизации ожиданий. Оно нам понадобится для обновления значений R, μ и Σ:

Теперь вы понимаете, почему перед этим мы специально вычисляли матрицу компонентов?

Имея значение переменной γ, мы можем определить обновления для R, μ и Σ. Они выглядят очень похоже на переменные обычной гауссовой смеси распределений:

Эти же принципы мы можем применить и для масштабирования, и множественных наблюдений, обсуждавшихся ранее.

Генерация данных из реальной скрытой марковской модели

Хотите узнать, как можно сгенерировать непрерывные данные из марковской модели?

Мы создадим две функции. Первую назовём simple_init; она будет представлять более простую модель с одним состоянием, одной гауссианой и одной размерностью:

import numpy as np

import matplotlib.pyplot as plt

def simple_init():

    M = 1

    K = 1

    D = 1

    pi = np.array([1])

    A = np.array([[1]])

    R = np.array([[1]])

    mu = np.array([[[0]]])

    sigma = np.array([[[[1]]]])

    return M, K, D, pi, A, R, mu, sigma

Теперь создадим функцию big_init, которая будет создавать более сложный сигнал. У неё пять состояний, три гауссианы и две размерности. Значения А я подготовил заранее. Это будет 0,9 в диагонали и 0,025 во всех остальных местах. R будет иметь равномерное распределение, значения μ я также подготовил заранее

def big_init():

    M = 5

    K = 3

    D = 2

    pi = np.array([1, 0, 0, 0, 0])

    A = np.array([

        [0.9, 0.025, 0.025, 0.025, 0.025],

        [0.025, 0.9, 0.025, 0.025, 0.025],

        [0.025, 0.025, 0.9, 0.025, 0.025],

        [0.025, 0.025, 0.025, 0.9, 0.025],

        [0.025, 0.025, 0.025, 0.025, 0.9],

    ])

    R = np.ones((M, K)) / K

    mu = np.array([

        [[0, 0], [1, 1], [2, 2]],

        [[5, 5], [6, 6], [7, 7]],

        [[10, 10], [11, 11], [12, 12]],

        [[15, 15], [16, 16], [17, 17]],

        [[20, 20], [21, 21], [22, 22]],

    ])

    sigma = np.zeros((M, K, D, D))

    for m in xrange(M):

        for k in xrange(K):

            sigma[m,k] = np.eye(D)

    return M, K, D, pi, A, R, mu, sigma

Следующее – функция get_signals, которая будет делать семплирование. При этом N – количество последовательностей, T – длина последовательности, init – какую именно из ранее написанных функций использовать.

def get_signals(N=20, T=100, init=big_init):

    M, K, D, pi, A, R, mu, sigma = init()

    X = []

    for n in xrange(N):

        x = np.zeros((T, D))

        s = 0

        r = np.random.choice(K, p=R[s])

        x[0] = np.random.multivariate_normal(mu[s][r], sigma[s][r])

        for t in xrange(1, T):

            s = np.random.choice(M, p=A[s])

            r = np.random.choice(K, p=R[s])

            x[t] = np.random.multivariate_normal(mu[s][r], sigma[s][r])

        X.append(x)

    return X

if __name__ == ‘__main__’:

    T = 500

    x = get_signals(1, T)[0]

    axis = range(T)

    plt.plot(axis, x[:, 0], axis, x[:, 1])

plt.show()

Запустим программу и посмотрим, что получится.

Видим сгенерированный сигнал.

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

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