Гауссова смесь распределений (GMM)

Описание гауссовой смеси распределений и как её обучать

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

Следующая тема, которую мы обсудим, – это гауссовы смеси распределений.

Гауссовы смеси распределений – это одна из форм оценки плотности распределения. Они дают нам приближённое значение распределения вероятностей наших данных. Гауссовы смеси распределений используются в случае, когда обнаруживается, что наши данные имеют мультимодальное распределение, то есть имеют несколько «горбов». Как вы помните из теории вероятностей, мода – это наиболее часто встречаемое значение.

Гауссова смесь распределений (GMM)

Гауссова смесь – это просто сумма всех взвешенных гауссиан. Для представления весовых коэффициентов вводится новый символ π. Например, πk означает вероятность того, что x принадлежит k-той гауссиане:

p(x)=\pi_1 N(\mu_1\sum_1)+\pi_2N(\mu_2\sum_2)+\pi_3N(\mu_3\sum_3) + ...

Обратите внимание, что существует ограничение – сумма всех π должна быть равной 1. Чтобы доказать это, используем тот факт, что все распределения вероятностей должны интегрироваться до 1. P(x) интегрируется до 1, и каждая гауссиана (поскольку мы знаем, что это действительные распределения) также интегрируется до 1. Таким образом, сумма всех π равна 1:

1 = \int p(x)dx = \int \pi_1 N((x|\mu_1,\sum_1)dx+\pi_2N((x|\mu_2,\sum_2)dx

1=\pi_1 * 1+\pi_2 * 1.

Другой способ представления – это введение новой скрытой переменной Z. Z показывает, какие данные представляет собой данная гауссиана. Можно записать:

\pi_k=P(Z=k).

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

Обучение гауссовой смеси распределений очень похоже на алгоритм метода k-средних. Оно проходит в два этапа, отражающих аналогичные этапы в методе k-средних.

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

\gamma_k^{(n)}=p(z^{(n)|}x)=\frac{\pi_k N (x^{(n)}|\mu_k,\sum_k)}{\sum^K_{j=1}\pi_j N(x^{(n)}|\mu_j,\sum_j).

Как вы можете видеть, если πk будет большим, то данная гауссиана перекроет остальные и станет примерно равной единице.

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

\mu_k=\frac{1}{N_k} \sum^N_{n=1} \gamma_k^{(n)}x^{(n)},

\sum_k=\frac{1}{N_k} \sum^N_{n=1} \gamma_k^{(n)}(x^{(n)}-\mu_k)(x^{(n)}-\mu_k)^T,

\pi_k = \frac {N_k}{N},  где N_k = \sum^N_{n=1} \gamma_k^{(n)}.

Сравнение гауссовой смеси распределений и метода k-средних

В этом разделе мы поэтапно сравним гауссовы смеси распределений с «мягким» методом k-средних. Это поможет нам не только оценить сходство между ними, но и понять ограничения метода k-средних и причину их существования.

Итак, давайте сравним два этапа каждого из алгоритмов:

Гауссова смесь распределений (GMM)

Как видите, у обоих первым этапом является расчёт принадлежностей, а вторым – расчёт параметров модели.

Теперь становится понятным, почему метод k-средних подразумевает поиск кластеров одинакового веса – потому что в нём нет переменной π, что эквивалентно равности всех значений π. Обратите внимание, что при этом гауссовы смеси распределений имеют то же ограничение, что и метод k-средних – вам и в этом случае придётся подбирать значение K.

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

Заметьте, что выражения для вычисления среднего значения при этом совершенно одинаковы.

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

Код для гауссовой смеси распределений на языке Python

В этой части мы рассмотрим код для гауссовой смеси распределений. Если вы не хотите сами писать код, а лишь посмотреть на уже написанный, перейдите по адресу https://github.com/lazyprogrammer/machine_learning_examples и в папке unsupervised_class найдите файл gmm.py.

Мы вновь начнём с импорта основных библиотек и функции main. Мы также используем те же самые данные, которые использовались перед этим для иерархической кластеризации и кластеризации методом k-средних. Таким образом, у нас будет три гауссовых распределения. Сразу же нарисуем диаграмму рассеяния, чтобы посмотреть, как выглядят данные, а также используем функцию gmm, которую опишем позже. Кроме того, нам необходимо импортировать функцию multivariate_normal из библиотеки SciPy.

Для лучшего отличия от метода k-средних параметры данных будут несколько изменены – например, количество примеров у нас станет равным 2 000, а сами данные будут иметь разную ковариации и плотность.

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import multivariate_normal

def main():

    D = 2

    s = 4

    mu1 = np.array([0, 0])

    mu2 = np.array([s, s])

    mu3 = np.array([0, s])

    N = 2000

    X = np.zeros((N, D))

    X[:1200, :] = np.random.randn(1200, D)*2 + mu1

    X[1200:1800, :] = np.random.randn(600, D) + mu2

    X[1800:, :] = np.random.randn(200, D)*0.5 + mu3

    plt.scatter(X[:,0], X[:,1])

    plt.show()

    K = 3

    gmm(X, K)

Теперь напишем код для работы с многоразмерными гауссианами. Первая часть очень похожа на код для метода k-средних, разве что у нас появляется новая переменная C для ковариации – это трёхмерная матрица, и одномерная переменная pi, представляющее наше π.

И сразу установим начальные средние расстояние случайным образом, а ковариацию – равной единице.

def gmm(X, K, max_iter=20):

    N, D = X.shape

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

    R = np.zeros((N, K))

    C = np.zeros((K, D, D))

    pi = np.ones(K) / K

    for k in xrange(K):

        M[k] = X[np.random.choice(N)]

        C[k] = np.diag(np.ones(D))

Теперь переходим к главному циклу нашей функции gmm. Первый этап очень похож на метод k-средних.

    costs = np.zeros(max_iter)

    weighted_pdfs = np.zeros((N, K))

    for i in xrange(max_iter):

        for k in xrange(K):

            for n in xrange(N):

                weighted_pdfs[n,k] = pi[k]*multivariate_normal.pdf(X[n], M[k], C[k])

        for k in xrange(K):

            for n in xrange(N):

                R[n,k] = weighted_pdfs[n,k] / weighted_pdfs[n,:].sum()

Второй же этап – это просто перерасчёт всех параметров – π, μ и ковариации. Тут тоже всё довольно похоже на метод k-средних, но есть и новый фактор – ковариация.

        for k in xrange(K):

            Nk = R[:,k].sum()

            pi[k] = Nk / N

            M[k] = R[:,k].dot(X) / Nk

            X_minus_mu = X – M[k]

            C[k] = np.sum(R[n,k]*np.outer(X_minus_mu, X_minus_mu) for n in xrange(N)) / Nk + np.diag(np.ones(D))*0.001

Теперь обновим значение функции затрат (ошибок). Заметьте, что это вовсе не та функция, что была в методе k-средних.

        costs[i] = np.log(weighted_pdfs.sum(axis=1)).sum()

        if i > 0:

            if np.abs(costs[i] – costs[i-1]) < 0.1:

                break

И последнее – вывести всё на экран. Мы, как и ранее, используем вывод в случайных цветах, так что эту часть можно скопировать из кода для метода k-средних.

    plt.plot(costs)

    plt.title(“Costs”)

    plt.show()

    random_colors = np.random.random((K, 3))

    colors = R.dot(random_colors)

    plt.scatter(X[:,0], X[:,1], c=colors)

    plt.show()

Возможно, вам будет интересно, сможет ли гауссова смесь распределений найти корректные значения для π, μ и ковариации, поэтому выведем и их значения на экран.

    print “pi:”, pi

    print “means:”, M

    print “covariances:”, C

Запустим нашу программу.

Первым видим наши гауссовы облака, несколько трудноразличимые.

Гауссова смесь распределений (GMM)

Далее – значение функции затрат. Она почти достигает нуля.

Гауссова смесь распределений (GMM)

И три наших кластера. Судя по виду, отмеченные правильно.

Гауссова смесь распределений (GMM)

Значения π, μ и ковариации тоже рассчитаны правильно.

Вопросы практического использования гауссовой смеси распределений с вырожденной ковариацией

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

Гауссова смесь распределений (GMM)

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

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

Одно из решений – использовать так называемую диагональную ковариационную матрицу. Она полезна не только для обхода проблемы вырожденной ковариации, но и для ускорения вычислений – имея диагональную ковариационную матрицу, очень легко найти обратную ей, потому она используется для проведения многомерных вычислений. Мы просто делим единицу на все элементы, не равные нулю (они располагаются на диагонали), и получаем обратные величины:

Гауссова смесь распределений (GMM)

Можете попробовать сделать это самостоятельно, чтобы убедиться, что в результате получится единичная матрица.

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

\sum_{ij} = E [(x_j-\mu_j) (x_j - \mu_j)].

Поскольку размерности независимы, можем записать:

\sum_{ij} = E [(x_i-\mu_i) (x_j - \mu_j)]= E[(x_i -\mu_i)] E [(x_j-\mu_j)] = (E(x_i)-\mu_i) (E(x_j) - \mu_j)] = (\mu_i - \mu_i)(\mu_j-\mu_j) = 0.

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

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

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

C = \frac{(X - mu)^T (X - mu)}{N},

где X – полная матрица размерностью NxD, muD-размерное среднее значение.

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

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