一个样本能同时像两个簇吗？高斯混合模型与 EM 软聚类

输入：X [N,D]，成分数 K
初始化：π [K]，μ [K,D]，Σ [K,D,D]

repeat:
    # E 步：必须在 log 空间计算
    log_joint[i,k] = log π[k] + log Normal(X[i] | μ[k], Σ[k])
    log_norm[i] = logsumexp(log_joint[i,:])
    Γ[i,k] = exp(log_joint[i,k] - log_norm[i])

    # M 步
    N_k[k] = sum_i Γ[i,k]
    π[k] = N_k[k] / N
    μ[k] = sum_i Γ[i,k] X[i] / N_k[k]
    Σ[k] = weighted covariance around new μ[k] + reg_covar · I

    lower_bound = mean_i log_norm[i]
until lower_bound improvement < tol or max_iter reached

输出：π、μ、Σ；推理时重新执行 E 步得到概率 [Q,K]

import numpy as np

X = np.array([[0.0], [1.0], [2.0]])               # [N=3,D=1]
means = np.array([[0.0], [2.0]])                  # [K=2,D=1]
variances = np.ones((2, 1))                       # [K,D]
weights = np.array([0.5, 0.5])                    # [K]

for _ in range(20):
    diff = X[:, None, :] - means[None, :, :]      # [N,K,D]
    log_gaussian = -0.5 * (
        np.log(2 * np.pi * variances)[None, :, :]
        + diff**2 / variances[None, :, :]
    ).sum(axis=2)                                 # [N,K]
    log_joint = np.log(weights)[None, :] + log_gaussian

    row_max = log_joint.max(axis=1, keepdims=True)
    log_norm = row_max + np.log(
        np.exp(log_joint - row_max).sum(axis=1, keepdims=True)
    )                                             # [N,1]
    responsibilities = np.exp(log_joint - log_norm)  # [N,K]

    effective_count = responsibilities.sum(axis=0)    # [K]
    weights = effective_count / X.shape[0]
    means = responsibilities.T @ X / effective_count[:, None]
    diff = X[:, None, :] - means[None, :, :]
    variances = (
        responsibilities[:, :, None] * diff**2
    ).sum(axis=0) / effective_count[:, None]
    variances = np.maximum(variances, 1e-6)

assert responsibilities.shape == (3, 2)
assert np.allclose(responsibilities.sum(axis=1), 1.0)
assert np.isclose(weights.sum(), 1.0)
assert np.all(variances > 0)

import numpy as np
from sklearn.mixture import GaussianMixture
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

model = make_pipeline(
    StandardScaler(),
    GaussianMixture(
        n_components=2,
        covariance_type='full',
        init_params='k-means++',
        n_init=10,
        reg_covar=1e-6,
        tol=1e-3,
        max_iter=300,
        random_state=42,
    ),
)
model.fit(X_train)                                # X_train [N,D]

gmm = model.named_steps['gaussianmixture']
prob = model.predict_proba(X_new)                  # [Q,K]
label = model.predict(X_new)                       # [Q]，argmax(prob)
log_density = model.score_samples(X_new)           # [Q]

assert prob.shape == (X_new.shape[0], 2)
assert np.allclose(prob.sum(axis=1), 1.0)
assert gmm.means_.shape == (2, X_train.shape[1])
assert gmm.covariances_.shape == (2, X_train.shape[1], X_train.shape[1])
print(gmm.converged_, gmm.n_iter_, gmm.lower_bound_)

from sklearn.mixture import GaussianMixture

candidates = []
for k in range(1, 7):
    gm = GaussianMixture(
        n_components=k,
        covariance_type='full',
        n_init=10,
        reg_covar=1e-5,
        random_state=42,
    ).fit(X_train_scaled)
    candidates.append({
        'k': k,
        'bic_train': gm.bic(X_train_scaled),
        'valid_log_likelihood': gm.score(X_valid_scaled),
        'smallest_weight': gm.weights_.min(),
        'converged': gm.converged_,
    })