Python中拟合固定协方差高斯混合模型的方法（GPS数据场景）

Great question! I’ve dealt with this exact challenge when working on GPS trajectory clustering—fixing the covariance to match known sensor noise is a smart move to avoid overfitting to noise and get more reliable cluster centers. Here’s a clean, maintainable way to do this without touching sklearn’s source code:

Solution: Custom GaussianMixture Subclass #

Instead of modifying sklearn’s source, we can inherit from GaussianMixture and override the M-step of the EM algorithm to force fixed covariance matrices. The EM algorithm has two core steps:

• E-step: Calculate the posterior probability of each sample belonging to each cluster.
• M-step: Update cluster weights, means, and covariances using the posterior probabilities.

We’ll let sklearn handle the E-step and most of the M-step logic, but overwrite the covariance update to use our pre-defined fixed values instead of estimating them from raw data.

Step-by-Step Implementation #

First, define your fixed covariance matrix. For 2D GPS data with a 3m standard deviation (isotropic noise), the covariance matrix is [[9, 0], [0, 9]] (since variance = std²). We’ll create a subclass that uses this fixed value for all clusters:

import numpy as np
from sklearn.mixture import GaussianMixture

class FixedCovarianceGaussianMixture(GaussianMixture):
    def __init__(self, fixed_covariances, **kwargs):
        super().__init__(**kwargs)
        # Store fixed covariances: shape depends on your chosen covariance_type
        # For 'full': (n_components, n_features, n_features)
        # For 'diag': (n_components, n_features)
        self.fixed_covariances = fixed_covariances

    def _m_step(self, X, log_resp):
        # Run standard M-step for weights and means (reuse sklearn's tested logic)
        n_samples, n_features = X.shape
        self.weights_ = np.mean(log_resp, axis=0)
        self.means_ = np.dot(log_resp.T, X) / np.sum(log_resp, axis=0)[:, np.newaxis]

        # Override covariances with our fixed values
        if self.covariance_type == 'full':
            self.covariances_ = self.fixed_covariances
        elif self.covariance_type == 'diag':
            # Extract diagonal values if using diagonal covariance format
            self.covariances_ = np.array([np.diag(cov) for cov in self.fixed_covariances])
        elif self.covariance_type == 'spherical':
            # Use average variance for spherical covariance setup
            self.covariances_ = np.array([np.trace(cov)/n_features for cov in self.fixed_covariances])
        elif self.covariance_type == 'tied':
            # For tied covariance, use the first fixed matrix (assume all are identical)
            self.covariances_ = self.fixed_covariances[0]

        # Normalize weights to sum to 1 (critical for valid probability outputs)
        self.weights_ /= self.weights_.sum()

How to Use the Custom Model #

Here’s how to apply this to your GPS data (we’ll simulate sample data for demonstration):

# Simulate 2D GPS data with two clusters (replace with your actual dataset)
np.random.seed(42)
n_samples = 200
cluster1 = np.random.multivariate_normal([5, 5], [[9, 0], [0, 9]], n_samples//2)
cluster2 = np.random.multivariate_normal([20, 8], [[9, 0], [0, 9]], n_samples//2)
X = np.vstack([cluster1, cluster2])

# Define fixed covariance matrices for 2 components (2D each)
fixed_cov = np.array([[[9, 0], [0, 9]], [[9, 0], [0, 9]]])

# Initialize and fit the model
model = FixedCovarianceGaussianMixture(
    fixed_covariances=fixed_cov,
    n_components=2,
    covariance_type='full',  # Matches the shape of our fixed_cov array
    random_state=42
)
model.fit(X)

# Extract results
print("Estimated cluster means:\n", model.means_)
print("\nFixed covariances (unchanged):\n", model.covariances_)

# Predict cluster labels for your raw GPS data
cluster_labels = model.predict(X)

Key Notes #

• Flexibility: If your GPS noise is anisotropic (different standard deviations in x/y directions), just adjust the diagonal values of fixed_cov (e.g., [[16,0],[0,9]] for 4m x-std and 3m y-std).
• Compatibility: This subclass inherits all features of the original GaussianMixture, so you can use methods like predict_proba(), score_samples(), and aic() exactly as you would with the standard model.
• Stability: By reusing sklearn’s built-in E-step and mean/weight calculation, you avoid reinventing the wheel and benefit from sklearn’s numerical stability checks and optimizations.

Alternative: Manual Likelihood Maximization #

If you prefer not to subclass, you could use scipy.optimize to maximize the log-likelihood directly, fixing covariance parameters. However, this requires implementing the EM algorithm from scratch, which is more error-prone and misses out on sklearn’s optimized workflows. The subclass method is far cleaner for most real-world use cases.

内容的提问来源于stack exchange，提问作者Ulf Aslak