在本例中，将探讨高斯混合模型（GMM）在非高斯随机变量混合的数据集上的表现。数据集由100个点组成，这些点大致按照带有噪声的正弦曲线分布。由于数据并非来自高斯随机变量的混合，因此没有关于高斯分量数量的真实值。
首先使用期望最大化（EM）算法拟合了一个包含10个分量的经典高斯混合模型。接着，利用变分推断方法拟合了一个具有狄利克雷过程先验的贝叶斯高斯混合模型，该模型的浓度先验值较低，使得模型倾向于使用较少的活跃分量。这种模型将重点放在数据集的整体结构上，即通过非对角协方差矩阵建模的具有交替方向的点群。这些交替方向大致捕捉了原始正弦信号的交替特性。
第三个模型也是一个具有狄利克雷过程先验的贝叶斯高斯混合模型，但这次浓度先验的值较高，给予模型更多的自由度来模拟数据的细粒度结构。结果是，这个模型具有更多的活跃分量，类似于随意决定固定分量数量为10的第一个模型。
哪个模型是最好的，这是一个主观判断的问题：希望模型只捕捉大部分数据结构的大图景，同时忽略细节，还是更喜欢那些紧密跟随信号高密度区域的模型？
import numpy as np
from sklearn import mixture
import matplotlib.pyplot as plt
from scipy import linalg
import itertools
# 定义颜色迭代器
color_iter = itertools.cycle(["navy", "c", "cornflowerblue", "gold", "darkorange"])
# 绘制结果的函数
def plot_results(X, Y, means, covariances, index, title):
splot = plt.subplot(5, 1, 1 + index)
for i, (mean, covar, color) in enumerate(zip(means, covariances, color_iter)):
v, w = linalg.eigh(covar)
v = 2.0 * np.sqrt(2.0) * np.sqrt(v)
u = w[0] / linalg.norm(w[0])
if not np.any(Y == i):
continue
plt.scatter(X[Y == i, 0], X[Y == i, 1], 0.8, color=color)
# 绘制高斯分量的椭圆
angle = np.arctan(u[1] / u[0])
angle = 180.0 * angle / np.pi # 转换为度
ell = plt.patches.Ellipse(mean, v[0], v[1], angle=180.0 + angle, color=color)
ell.set_clip_box(splot.bbox)
ell.set_alpha(0.5)
splot.add_artist(ell)
plt.xlim(-6.0, 4.0 * np.pi - 6.0)
plt.ylim(-5.0, 5.0)
plt.title(title)
plt.xticks(())
plt.yticks(())
# 参数设置
n_samples = 100
# 按照正弦曲线生成随机样本
np.random.seed(0)
X = np.zeros((n_samples, 2))
step = 4.0 * np.pi / n_samples
for i in range(X.shape[0]):
x = i * step - 6.0
X[i, 0] = x + np.random.normal(0, 0.1)
X[i, 1] = 3.0 * (np.sin(x) + np.random.normal(0, 0.2))
plt.figure(figsize=(10, 10))
plt.subplots_adjust(bottom=0.04, top=0.95, hspace=0.2, wspace=0.05, left=0.03, right=0.97)
# 使用EM算法拟合高斯混合模型，包含10个分量
gmm = mixture.GaussianMixture(n_components=10, covariance_type="full", max_iter=100)
gmm.fit(X)
plot_results(X, gmm.predict(X), gmm.means_, gmm.covariances_, 0, "期望最大化")
# 使用变分推断拟合具有狄利克雷过程先验的贝叶斯高斯混合模型
dpgmm = mixture.BayesianGaussianMixture(n_components=10, covariance_type="full", weight_concentration_prior=1e-2, weight_concentration_prior_type="dirichlet_process", mean_precision_prior=1e-2, covariance_prior=1e0 * np.eye(2), init_params="random", max_iter=100, random_state=2)
dpgmm.fit(X)
plot_results(X, dpgmm.predict(X), dpgmm.means_, dpgmm.covariances_, 1, "具有狄利克雷过程先验的贝叶斯高斯混合模型，浓度先验为0.01")
# 绘制采样结果
X_s, y_s = dpgmm.sample(n_samples=2000)
plt.subplot(5, 1, 4)
for i, color in zip(range(dpgmm.n_components), color_iter):
if not np.any(y_s == i):
continue
plt.scatter(X_s[y_s == i, 0], X_s[y_s == i, 1], 0.8, color=color)
plt.xlim(-6.0, 4.0 * np.pi - 6.0)
plt.ylim(-5.0, 5.0)
plt.title("具有狄利克雷过程先验的贝叶斯高斯混合模型，浓度先验为0.01，采样2000个样本")
plt.show()