Гауссово моделирование смеси (GMM)
Дата публикации Mar 8, 2019
В предыдущемпослеЯ обсуждал кластеризацию k-средних как способ обобщения текстовых данных. Я также говорил о некоторых ограничениях k-средних и в каких ситуациях это может быть не самым подходящим решением. Вероятно, самым большим ограничением является то, что каждый кластер имеет одинаковую диагональную ковариационную матрицу. Это создает сферические кластеры, которые являются довольно негибкими с точки зрения типов распределений, которые они могут моделировать. В этом посте я хотел бы рассмотреть некоторые из этих ограничений и рассказать об одном конкретном методе, который может избежать этих проблем,Гауссово моделирование смесей (GMM).Формат этого поста будет очень похож на последний, где я объясняю теорию GMM и как она работает. Затем я хочу погрузиться в кодирование алгоритма на Python, и мы увидим, как результаты отличаются от k-средних и почему использование GMM может быть хорошей альтернативой.
GMM сделано просто (иш)
В самом простом случае, GMM также является типом алгоритма кластеризации. Как следует из названия, каждый кластер моделируется в соответствии с различным распределением Гаусса. Этот гибкий и вероятностный подход к моделированию данных означает, что вместо жестких назначений в кластеры, такие как k-средних, у нас есть мягкие назначения. Это означает, что каждая точка данных могла быть сгенерирована любым из распределений с соответствующей вероятностью. По сути, каждый дистрибутив несет некоторую «ответственность» за создание определенной точки данных.
Как мы можем оценить этот тип модели? Ну, одна вещь, которую мы могли бы сделать, это ввестискрытая переменнаяγ(гамма) для каждой точки данных. Это предполагает, что каждая точка данных была сгенерирована с использованием некоторой информации о скрытой переменнойγ, Другими словами, он говорит нам, какой гауссов генерировал конкретную точку данных. На практике, однако, мы не наблюдаем эти скрытые переменные, поэтому нам нужно оценить их. как нам это сделать? Что ж, к счастью для нас уже есть алгоритм, который можно использовать в таких случаях,Алгоритм максимизации ожидания (EM)и это то, что мы обсудим дальше.
EM алгоритм
Приведенное выше уравнение часто приводит к сложной функции, которую трудно максимизировать. Что мы можем сделать в этом случае, это использоватьНеравенство Дженсена построитьнижняя границафункция что гораздо проще оптимизировать. Если мы оптимизируем это путем минимизацииДивергенция КЛ (разрыв) между двумя распределениями мы можем аппроксимировать исходную функцию. Этот процесс показан на рисунке 1 ниже. Я также предоставил ссылку на видео выше, которая показывает происхождение дивергенции KL для тех из вас, кто хочет более строгое математическое объяснение.
Чтобы оценить нашу модель по существу, нам нужно всего лишь выполнить два шага. На первом шаге (E-шаг) мы хотим оценить апостериорное распределение наших скрытых переменных 𝛾 в зависимости от наших весов (π), средних (µ) и ковариации (Σ) наших гауссианов. Вектор параметров обозначен как θ на рисунке 1. Оценка E-шага требует сначала инициализации этих значений, и мы можем сделать это с помощью k-средних, что обычно является хорошей отправной точкой (подробнее об этом в коде ниже). Затем мы можем перейти ко второму шагу (M-шагу) и использовать 𝛾, чтобы максимизировать вероятность относительно наших параметров θ. Этот процесс повторяется до тех пор, пока алгоритм не сходится (функция потерь не изменяется).
Визуализация EM-алгоритма
Почему бы нам не попытаться визуализировать этот процесс, используя рисунок 1? Мы вычисляем апостериорное распределение 𝛾 на первом этапе, которое, как оказывается, эквивалентно значению, которое мы получили бы, минимизируя расхождение KL между этими двумя распределениями. Затем мы устанавливаем апостериор равным q (запутанная запись, которую я знаю, но это всего лишь 𝛾) и максимизируем эту функцию по отношению к параметрам θ. Мы можем видеть из графика, когда мы повторяем и выполняем эти вычисления, мы движемся к оптимуму (или, по крайней мере, к локальному оптимуму).
EM-алгоритм для GMM
E-Step
M-Step
После вычисления нашего апостериорного значения все, что нам нужно сделать, это получить оценку параметров каждого гауссиана, определенных уравнениями ниже, а затем оценить логарифмическую вероятность. Эти два шага затем повторяются до сходимости.
Код Python
Теперь, когда мы объяснили теорию моделирования, я хочу кодировать этот алгоритм с помощью Python. Как и мой предыдущийпослеЯ собираюсь использовать тот же набор данных, чтобы мы могли сравнить результаты между k-means и GMM. Шаги предварительной обработки точно такие же, как и в предыдущем посте, и я предоставляю ссылку на полный код в конце этого поста.
Оценка K-средних
Как я упоминал ранее, чтобы запустить алгоритм (выполнить 1-й шаг E), нам нужны начальные значения для наших параметров. Вместо того, чтобы просто устанавливать эти значения случайным образом, обычно хорошей идеей является их оценка с использованием k-средних. Это обычно дает нам хорошую отправную точку и может помочь нашей модели быстрее сходиться. Прежде чем мы оценим GMM, давайте кратко рассмотрим, какие кластеры нам дает k-means.
Используя наши оценки из sklearn, мы можем создать красивую визуализацию наших кластеров (рисунок 2). Обратите внимание, что кластеры имеют сферическую форму и одинаковый размер. Сферические кластеры, по-видимому, не очень хорошо моделируют распространение данных, что указывает на то, что k-среднее в этом конкретном случае может быть не лучшим подходом. Это иллюстрирует одно из ограничений k-средних, поскольку все ковариационные матрицы диагональны с единичной дисперсией. Это ограничение означает, что модель не особенно гибкая. Имея это в виду, давайте попробуем GMM и посмотрим, какие результаты это дает нам.
Оценка GMM
Рисунок 3 ниже иллюстрирует, что делает GMM. Это ясно показывает три кластера, смоделированных тремя различными гауссовыми распределениями. Я использовал игрушечный набор данных здесь только для того, чтобы проиллюстрировать это наглядно, поскольку он менее понятен с набором данных Enron. Как вы можете видеть, по сравнению с рисунком 2, смоделированным с использованием сферических кластеров, GMM гораздо более гибок, что позволяет нам генерировать гораздо более подходящие распределения.
Класс GMM Python
Хорошо, теперь мы приступим к написанию кода нашего класса GMM на Python. Как всегда, мы начинаем с метода init. Единственное, что я инициализирую здесь, это количество раз, которое мы хотим запустить наш алгоритм, и количество кластеров, которые мы хотим смоделировать. Самый интересный метод в этом фрагменте кодаcalculate_mean_covariance, Это помогает нам рассчитать значения для наших начальных параметров. Он принимает наши данные, а также наши прогнозы от k-средних и вычисляет весовые коэффициенты, средние значения и ковариационные матрицы для каждого кластера.
Следующий фрагмент кода реализует наш метод initialise_parameters, который использует k-средства из библиотеки sklearn для вычисления наших кластеров. Обратите внимание, что эта функция на самом деле вызывает наш метод Calculate_mean_covariance, определенный выше. Возможно, мы могли бы использовать один метод для расчета наших кластеров и исходных параметров, но обычно гораздо проще отлаживать и избегать ошибок, если каждый метод выполняет только одну конкретную задачу.
После того, как мы вычислили это значение для каждого гауссиана, нам просто нужно нормализовать гамму (𝞬), соответствующую знаменателю в уравнении 3. Это сделано для того, чтобы наши гамма были действительными вероятностями. Если мы суммируем значения по кластерам для каждой точки данных, они должны равняться 1.
После того, как мы вычислим значения для обязанностей (𝞬), мы можем подать их вМ-шагом.Снова цель M-шага состоит в том, чтобы вычислить наши новые значения параметров, используя результаты E-шага, соответствующие уравнениям 4, 5 и 6. Чтобы упростить отладку, я разделилm_stepметод иcompute_loss_functionМетод в моем коде ниже. Функция compute_loss_function делает именно то, что подразумевает ее имя. Он принимает на себя обязанности и параметры, возвращаемые E-шагом и M-шагом, и использует их для вычисления нашей функции потерь с нижней границей, определенной в уравнении 9.
Все наши самые важные методы теперь закодированы. В соответствии со sklearn я собираюсь определить метод подгонки, который будет вызывать методы, которые мы только что определили. В частности, мы начнем с инициализации значений наших параметров. После этого мы выполняем шаги, описанные в EM-алгоритме для выбранного количества итераций. Обратите внимание, что на самом деле не требуется большого количества итераций, чтобы сходиться, особенно когда вы используете k-средства для получения значений начальных параметров (я думаю, что мой алгоритм сходился примерно за 30 итераций).
Поскольку мы, вероятно, также заинтересованы в использовании этой модели для прогнозирования того, к чему могут принадлежать новые гауссовские данные, мы можем реализоватьпрогнозироватьа такжеpredict_probaметод. Метод предиката_проба будет принимать новые точки данных и прогнозировать ответственность для каждого гауссиана. Другими словами, вероятность того, что эта точка данных пришла от каждого распределения. В этом суть мягкого задания, о котором я упоминал в начале поста. Метод прогнозирования, по сути, делает то же самое, но назначает кластер с наибольшей вероятностью, используя np.argmax.
Примерка нашей модели
После этого объяснения я думаю, что пришло время оценить нашу модель и посмотреть, что мы получим. Надеемся, что представленная выше визуализация GMM дала хорошее представление о том, что делает модель. Мы собираемся сделать то же самое для нашего набора данных Enron. Код ниже просто оценивает нашу модель GMM в нашем наборе данных, используя 3 разных гауссиана. Для построения графика я также вычисляю точку наибольшей плотности каждого распределения, соответствующую центру, который полезен в качестве средства визуализации. Наконец, мы также используем параметры модели, чтобы нарисовать форму каждого распределения на рисунке 4.
Основным выводом на этом рисунке является то, что распределения явно больше не являются сферическими. GMM позволил нам ослабить наши ограничения на ковариационную матрицу, что позволило распределению лучше соответствовать данным. Это особенно полезно, учитывая, что форма наших данных была явно не сферической. Теперь это, вероятно, не идеальное решение, и есть некоторые точки данных, которые не очень хорошо подходят для любого распределения, но это улучшение по сравнению с k-средних.
Внедрение GMM sklearn
Теперь, просто чтобы убедиться, что мы не сделали ничего сумасшедшего в нашем коде, я собираюсь повторить эту оценку, используя sklearn, и посмотреть, будет ли мое решение таким же. Приведенный ниже код в значительной степени совпадает с приведенным выше кодом, поэтому я не буду подробно останавливаться на нем. Похоже, у нас очень похожий результат по сравнению со sklearn. Единственное отличие состоит в том, что один из наших кластерных центров кажется другим. В реализации sklearn центр ближе к 0.4, а в нашей реализации ближе к 0.6. Возможно, это связано с немного другой инициализацией в sklearn?
Хорошо, ребята, вот и все для этого поста. Я надеюсь, что это было полезным и довольно интуитивно понятным объяснением гауссового моделирования смесей. Если кто-то из вас хочет глубже понять материал, я рекомендую курс CourseraБайесовские методы машинного обучения,Я получил много материала из этого курса, и я думаю, что он дает действительно хорошие и подробные объяснения представленных здесь концепций. Я также рекомендовал бы книгу,Распознавание образов и машинное обучение епископом,Эта книга является отличным справочником по большинству классических алгоритмов, с которыми вы столкнетесь в машинном обучении. Ниже я приведу полный код для класса GMM, описанный в посте, а также ссылку на ядро Kaggle, где я провел весь анализ. Как всегда, обратная связь приветствуется
Полный код Python для класса GMM
Источник:Кристофер М. Бишоп 2006, Распознавание образов и машинное обучение
Gaussian Mixture Models Visually Explained
A normal distribution is widely known and used in statistics. Because of the central limit theorem the normal distribution explains the data well in many cases. However, in practice the data distribution might be multimodal, i.e. the probability density function has several ‘bumps’. For example, the distribution of hardcover and paperback books prices, as producing a hardcover book is more expensive.
One way to exploit the nice properties of normal distributions for such multimodal data is to use the Gaussian Mixture Model. This combines several normal distributions with different parameters.
The Gaussian Mixture Model is a generative model that assumes the data is distributed as a Gaussian mixture. It can be used for density estimation and clustering. But, first things first.
Gaussian Mixture
The Gaussian Mixture Model defines a probability distribution on the data of the specific form — the mixture of Gaussians. Though you may rightly ask, “What is the Gaussian mixture exactly?” Gaussian mixture distribution is a linear superposition of m Gaussian components. This is not to be confused with a sum of random variables distributed as Gaussians — this sum is distributed as Gaussian too. Here we sum up not the random variables, but probability density functions with weight coefficients which summing up to 1.
This way we obtain the distribution which comprises of multiple Gaussian components. If the means are far enough, then we might get a distribution with several ‘bumps’, in other words, multimodal distribution. The parameters of this distribution are the weights π, means μ, and covariance matrices Σ for each component. The weight coefficients, also called mixing parameters, show how much of the respective component is in the resulting distribution.
To get a better intuition of this, let’s consider the simplest example; when one coefficient is one and all others are zero. In this case, the resulting distribution would be just Gaussian with respective mean and covariance. In a non-degenerative case, the resulting distribution captures all components to some extent. In the figure below, you can see how the resulting distribution changes depending on the value of the mixture parameter.
We defined this distribution as a combination of multiple Gaussian probability density functions. Equivalently, we can define it as a result of the following sampling procedure.
We first choose a component, and then sample a random variable from the chosen component. In more detail this might look like this:
The sampling procedure is shown in the animation below. The right subplot shows iteratively how many samples are chosen from each component. The resulting histogram is close to the barplot of the weights of components. The left subplot shows the samples drawn from the respective component’s Gaussian distribution. The components contain a different number of points. For example, the upper one (in histogram component 1) contains the least points.
Gaussian Mixture Model
Now imagine we know (or at least assume) the data is generated from the Gaussian mixture. However, the parameters of the distribution remain unknown. How do we learn the parameters? As often in machine learning models, this is done via likelihood maximization (or equivalently, negative log-likelihood minimization).
If we had only one component (that would be just a Gaussian distribution) we could take a log. (a gentle reminder: logarithm of the product is just a sum log(ab) = log(a) + lob(b) for a, b > 0). When we do this the exponent would disappear. Then we take a derivative, set it to zero, and the problem is solved (of course, we noticed that the log-likelihood is concave; hence, the solution is a maximizer).
However, in more general cases we have more than one component, and we don’t know the mixing coefficients. Although we still can set derivatives to zero, we won’t obtain a nice closed-form solution.
Let’s note, that if we knew which component a point is from (the class labels) we could learn the parameters of each component independently. Just like how things work when performing the density estimation of normally distributed data. On the other hand, if we knew the parameters of all components, we could assign points to their most probable classes. So, these subproblems are relatively easy to solve. This could lead us to the idea of updating parameters or class labels iteratively based on the previous estimate of class labels or components parameters, respectively. The idea of the EM algorithm is similar, but instead of class labels, it uses responsibilities — probabilities of a point belonging to a cluster.
EM Algorithm
Before we deep dive into the EM algorithm details, I’d like to remind you that this is only one of the approaches to find a maximizer. It can be shown that the algorithm converges, although not always to the global maximum. Another alternative would be, for example, gradient-based approaches.
To get further intuition on this, let’s consider an example. Imagine that centers of K clusters are kindergarten teachers looking after points under their responsibility. The teacher’s responsibility for a point depends on the point location, their location (μ), the area they frequently check (zone of responsibility, Σ), their commitment (π), and the same parameters of other teachers. The responsibility may be shared between teachers (for example, if a point is between two centers), but the responsibilities for a point of all teachers should sum up to 1. In other words, every point needs to be fully taken care of.
The contours around the teachers are level lines indicating how well they look after points on that line. The further the contour, the more difficult is to check on that location.
The goal is to find locations, zones of responsibility, and teachers’ commitment to increase overall ‘childcare quality’, which is an informal interpretation of likelihood. Another goal is to assign every point to a teacher, whose responsibility for it is largest among other teachers to find groups of points.
The teachers can easily compute their responsibility for all points using the current parameters of all teachers (expectation step).
Then they change the location and their zone of responsibility to better look after the points under their responsibility, moving closer to the points they are mainly responsible for and accounting with lower weight for points they have shared responsibility for (maximization step). Also, their commitment grows with the overall current responsibility.
Important point: Even if their responsibility for a point is the lowest among all other teachers (for example 0.2 and 0.8 for another teacher) they still account for it. The objective that they maximize is a weighted version of the log-likelihood with responsibilities acting as weights. The lower the responsibility, the less the respective summand. Hence the lower effect on the parameters of the cluster with low responsibility.
After changing the parameters, they recalculate their responsibilities and move again. The process is repeated until convergence (the changes are so small that the teachers are too lazy to move). It can be shown that the overall ‘childcare quality’ is increasing after each expectation-maximization iteration.
Now let’s move to the formal description of the EM algorithm.
Expectation Step
On the expectation step, we compute responsibilities. The responsibility r_ik is the posterior probability that point i belongs to cluster k. It is calculated using the Bayes’ theorem.
Maximization Step
As mentioned before, the log-likelihood cannot be maximized in closed form. However, the problem becomes easy once we have class labels. So let’s introduce latent variables Z, i = 1…n, with Z_i being a class label of point i. If we observed both X and Z we could easily deduce parameters of every Gaussian component by maximizing the log-likelihood. However, in our case Z depends on the parameters. Another consideration is that GMM is a probabilistic model, and even if a point is close to the center of a cluster, it could still arise from another cluster with lower probability. When we make inference, we choose the most probable class of each point, but even under GMM assumptions, it doesn’t mean we chose the ‘true’ class.
So now we have the log-likelihood function we want to maximize but cannot compute because of unknown Z… Don’t worry; there’s a solution. One way to estimate a quantity containing randomness is to take its expectation. Although we don’t know Z, we have estimates of its distribution — that’s r_ik computed on the expectation step. We can take expectation of log-likelihood with respect to the posterior distribution of cluster assignments.
In the second line, we used that only one indicator z_i = k is one, and all others are zero. Therefore, the product contains only one multiplier referring to the true class. The reason to use a product is that we can take it out of the logarithm and obtain a sum.
The maximizers now can be found by setting derivatives to zero and have the following form:
The EM algorithm monotonically increases log-likelihood and converges to local maximum. However, this algorithm depends on initialization and is not guaranteed to find the global likelihood maximizer.
Practice
It is a good exercise to implement the EM algorithm for GMM from scratch to check and deepen your understanding. But here, I want to give a quick overview of sklearn implementation usage.
First, importing libraries and setting matplotlib parameters up:
I used simulated data for simplicity and visual clarity.
The GaussianMixture() function creates an object, which we then fit to the data to learn the parameters.
Furthermore, sklearn implementation also allows for regularization, e.g., using the same covariance matrices for all components. It can be done by setting covariance_type=’tied’ in the initialization of GMM.
GMM is a generative model, which means that it can generate new data.
Summary
All animations were made with Python (check my gmm_visualisation repository if you want to see the implementation). The drawings I created in Procreate.
