如何调整高斯混合模型(GMM)以准确拟合一维临床数据？

数据背景 #

我有一维临床分析数据（如血液白细胞数据），数据分布呈现双模态特征，推测服从双高斯分布——分别对应重症患者与其余患者，其中一个分布集中在中间区域，另一个在数值较低的右侧区域。

已尝试的方法与问题 #

1. GMM拟合分类 #

使用sklearn的GaussianMixture进行双组分拟合：

gmm = GaussianMixture(n_components=2, random_state=42)
gmm.fit(data['obs_value_nm'].values.reshape(-1,1))
data['target_class'] = gmm.predict(data['obs_value_nm'].values.reshape(-1,1))

绘制分类前后的KDE图后，发现分类结果不符合预期，未能准确区分两个对应患者群体的分布。

2. 分布参数可视化 #

尝试用scipy.stats将拟合的高斯分布曲线绘制到直方图上：

x = np.linspace(group['Min'].values[0], group['Max'].values[0], 400)
y = norm.pdf(x, loc = stat['media'].values[0], scale = stat['std'].values[0])

但最终拟合曲线与原始直方图的贴合度较差，未能准确呈现双高斯分布特征。

核心需求 #

• 准确获取两个高斯分布的均值与标准差
• 基于这些参数，通过均值距离和偏差比较识别数据中的异常值
• 将两个高斯分布准确绘制在原始数据的直方图上

解决方案 #

1. 优化GMM拟合效果 #

GMM拟合不准确多源于局部最优或初始化问题，可通过以下方式调整：

• 数据标准化：先对数据做标准化处理，帮助模型收敛到合理组分：

  from sklearn.preprocessing import StandardScaler
  scaler = StandardScaler()
  scaled_data = scaler.fit_transform(data['obs_value_nm'].values.reshape(-1,1))
  # 增加初始化次数避免局部最优
  gmm = GaussianMixture(n_components=2, random_state=42, n_init=10)
  gmm.fit(scaled_data)
  

• 手动指定初始均值：若从直方图观察到两个分布的大致峰值，可手动设置初始均值引导拟合：

  # 替换为你观察到的两个峰值近似值
  init_means = np.array([[2], [8]])
  gmm = GaussianMixture(n_components=2, random_state=42, means_init=init_means, n_init=10)
  gmm.fit(scaled_data)
  

2. 提取原始尺度的分布参数 #

拟合完成后，将标准化后的参数反变换回原始数据尺度：

# 提取两个组分的均值（反变换）
means = scaler.inverse_transform(gmm.means_)[:, 0]
# 提取两个组分的标准差（协方差矩阵开根号后乘以标准化缩放因子）
stds = scaler.scale_ * np.sqrt(gmm.covariances_)[:, 0, 0]
# 提取每个组分的样本占比权重
weights = gmm.weights_

3. 准确绘制直方图+双高斯曲线 #

结合原始数据直方图，叠加带权重的拟合高斯曲线（权重匹配样本占比，保证曲线高度与直方图密度一致）：

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

plt.figure(figsize=(12,6))
# 绘制原始数据直方图（density=True保证与PDF曲线刻度一致）
counts, bins, _ = plt.hist(data['obs_value_nm'], bins=30, density=True, alpha=0.5, label='原始数据')

# 生成覆盖数据范围的x轴
x = np.linspace(bins[0], bins[-1], 400)

# 绘制两个高斯分布曲线
for mean, std, weight in zip(means, stds, weights):
    # 乘以权重匹配直方图密度
    y = weight * norm.pdf(x, loc=mean, scale=std)
    plt.plot(x, y, linewidth=2, label=f'高斯组分 (μ={mean:.2f}, σ={std:.2f})')

plt.xlabel('观测值')
plt.ylabel('密度')
plt.title('原始数据直方图+双高斯拟合')
plt.legend()
plt.show()

4. 基于参数识别异常值 #

可通过两种方式标记异常值：

• 3σ原则：对每个组分，标记超出μ ± 3σ范围的数据点为异常值：

  # 先给每个样本分配对应的组分
  data['component'] = gmm.predict(scaled_data)
  # 标记异常值
  for idx in range(2):
      mask = data['component'] == idx
      lower = means[idx] - 3*stds[idx]
      upper = means[idx] + 3*stds[idx]
      data.loc[mask, 'is_outlier'] = (data['obs_value_nm'] < lower) | (data['obs_value_nm'] > upper)
  

• 概率阈值法：计算样本属于GMM模型的概率，低于阈值（如0.01）的视为异常：

  # 计算样本的对数似然，转换为概率
  log_likelihood = gmm.score_samples(scaled_data)
  probabilities = np.exp(log_likelihood)
  data['is_outlier'] = probabilities < 0.01
  

内容的提问来源于stack exchange，提问作者jai melcar