在 Python 中拟合具有固定协方差的高斯混合答案

【问题标题】：Fit mixture of Gaussians with fixed covariance in Python在 Python 中拟合具有固定协方差的高斯混合
【发布时间】：2018-07-08 04:49:48
【问题描述】：

我有一些带有集群（停止位置）的二维数据（GPS 数据），我知道这些数据类似于具有特征标准偏差（与 GPS 样本的固有噪声成正比）的高斯。下图可视化了一个我期望有两个这样的集群的样本。图像宽 25 米，高 13 米。

sklearn 模块有一个函数sklearn.mixture.GaussianMixture，它允许您将高斯混合拟合到数据中。该函数有一个参数covariance_type，它使您能够假设关于高斯形状的不同事物。例如，您可以使用 'tied' 参数假设它们是统一的。

但是，假设协方差矩阵保持不变似乎并不直接可能。从sklearn源代码来看，进行修改以启用此功能似乎微不足道，但通过允许此功能的更新发出拉取请求感觉有点过分（我也不想在sklearn中意外添加错误） .在每个高斯的协方差矩阵固定的情况下，是否有更好的方法可以将混合拟合到数据？

我想假设每个组件的 SD 应保持恒定在 3 米左右，因为这大致是我的 GPS 样本的噪声水平。

【问题讨论】：

StackOverflow 鼓励发布一个由 MCVE 制定的问题。（其中 MCVE 代表 Minimum Complete Verifiable Example）。您介意发布为此收集的完整数据集吗？如果没有 DataSET，MCVE 只是在-C-完整且不可-V-可验证的，是吗？
@user3666197 感谢您的意见！不幸的是，我不能分享完整的数据集，它基于人类的流动性并且是高度机密的。另外，我觉得提供一个解决方案示例会相当简单，该解决方案只是模拟一些具有适当特性的集群。

标签： python machine-learning scikit-learn gmm

【解决方案1】：

编写自己的EM algorithm 实现非常简单。它还可以让您对该过程有一个很好的直觉。我假设协方差是已知的，并且分量的先验概率是相等的，并且仅拟合均值。

这个类看起来像这样（在 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)）。

【讨论】：

我们去！超级好的解决方案。抱歉，我只是将赏金授予@Xochipilli，但您的正是我一直在寻找的解决方案。非常感谢！

【解决方案2】：

我认为最好的选择是 "roll your own" GMM 模型，方法是定义一个新的 scikit-learn 类，该类继承自 GaussianMixture 并覆盖方法以获得您想要的行为。这样您就可以自己实现，而不必更改 scikit-learn 代码（并创建拉取请求）。

另一个可行的方法是查看 scikit-learn 中的 Bayesian version of GMM。您也许可以设置协方差矩阵的先验，以便固定协方差。它似乎使用Wishart distribution 作为协方差的先验。但是我对这个发行版不够熟悉，无法为您提供更多帮助。

【讨论】：

【解决方案3】：

首先，您可以使用spherical 选项，这将为您提供每个组件的单一方差值。这样你就可以检查自己，如果收到的方差值相差太大，那么就出了问题。

在您想要预设方差的情况下，您的问题会退化为只为您的组件找到最佳中心。例如，您可以使用k-means 来实现。如果您不知道组件的数量，您可以扫描所有逻辑值（如 1 到 20）并评估拟合误差的减少量。或者您可以优化自己的 EM 函数，以同时找到中心和分量数。

【讨论】：