Построить смесь гауссиан с фиксированной ковариацией в Python

У меня есть некоторые 2D-данные (данные GPS) с кластерами (места остановки), которые, как я знаю, напоминают гауссианцы с характерным стандартным отклонением (пропорциональным присущему шуму образцов GPS). На рисунке ниже представлен пример, который, как я полагаю, имеет два таких кластера. Изображение имеет ширину 25 метров и высоту 13 метров.

Модуль sklearn имеет функцию sklearn.mixture.GaussianMixture, которая позволяет сопоставлять смесь гауссианов с данными. Функция имеет параметр covariance_type, который позволяет вам принимать разные вещи о форме гауссовых. Например, вы можете считать их равномерными с помощью аргумента 'tied'.

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

Я хочу предположить, что SD должен оставаться постоянным на расстоянии около 3 метров для каждого компонента, так как это примерно соответствует уровню шума моих образцов GPS.

Ответ 1

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

Класс будет выглядеть так (в Python 3):

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

class FixedCovMixture:
    """ The model to estimate gaussian mixture with fixed covariance matrix. """
    def __init__(self, n_components, cov, max_iter=100, random_state=None, tol=1e-10):
        self.n_components = n_components
        self.cov = cov
        self.random_state = random_state
        self.max_iter = max_iter
        self.tol=tol

    def fit(self, X):
        # initialize the process:
        np.random.seed(self.random_state)
        n_obs, n_features = X.shape
        self.mean_ = X[np.random.choice(n_obs, size=self.n_components)]
        # make EM loop until convergence
        i = 0
        for i in range(self.max_iter):
            new_centers = self.updated_centers(X)
            if np.sum(np.abs(new_centers-self.mean_)) < self.tol:
                break
            else:
                self.mean_ = new_centers
        self.n_iter_ = i

    def updated_centers(self, X):
        """ A single iteration """
        # E-step: estimate probability of each cluster given cluster centers
        cluster_posterior = self.predict_proba(X)
        # M-step: update cluster centers as weighted average of observations
        weights = (cluster_posterior.T / cluster_posterior.sum(axis=1)).T
        new_centers = np.dot(weights, X)
        return new_centers


    def predict_proba(self, X):
        likelihood = np.stack([multivariate_normal.pdf(X, mean=center, cov=self.cov) 
                               for center in self.mean_])
        cluster_posterior = (likelihood / likelihood.sum(axis=0))
        return cluster_posterior

    def predict(self, X):
        return np.argmax(self.predict_proba(X), axis=0)

В данных, подобных вашей, модель будет сходиться быстро:

np.random.seed(1)
X = np.random.normal(size=(100,2), scale=3)
X[50:] += (10, 5)

model = FixedCovMixture(2, cov=[[3,0],[0,3]], random_state=1)
model.fit(X)
print(model.n_iter_, 'iterations')
print(model.mean_)

plt.scatter(X[:,0], X[:,1], s=10, c=model.predict(X))
plt.scatter(model.mean_[:,0], model.mean_[:,1], s=100, c='k')
plt.axis('equal')
plt.show();

и вывод

11 iterations
[[9.92301067 4.62282807]
 [0.09413883 0.03527411]]

Вы можете видеть, что расчетные центры ((9.9, 4.6) и (0.09, 0.03)) близки к истинным центрам ((10, 5) и (0, 0)).

Ответ 2

Я думаю, что лучшим вариантом будет "roll your own" Модель GMM, определяя новый класс scikit-learn, который наследуется от GaussianMixture и перезаписывает методы для получения желаемого поведения. Таким образом, у вас просто есть реализация самостоятельно, и вам не нужно менять код scikit-learn (и создавать запрос на pull).

Другим вариантом, который может работать, является просмотр байесовской версии GMM в scikit-learn. Возможно, вы сможете установить предшествующую для ковариационной матрицы так, чтобы ковариация была исправлена. Похоже, что использовать дистрибутив Wishart в качестве предварительного ковариации. Однако я недостаточно осведомлен об этом дистрибутиве, чтобы помочь вам больше.

Ответ 3

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

В случае, когда вы хотите настроить дисперсию, проблема ухудшается при поиске только лучших центров для ваших компонентов. Вы можете сделать это, используя k-means, например. Если вы не знаете количество компонентов, вы можете перебирать все логические значения (например, от 1 до 20) и оценивать декремент в ошибке установки. Или вы можете оптимизировать свою собственную функцию EM, найти центры и количество компонентов одновременно.