一、高斯混合模型介绍高斯混合模型GMM通过高斯分布描述每个类别的数据分布。与K-means仅提供中心点不同GMM采用概率模型实现soft assignment可计算每个数据点属于各类别的概率。高斯混合模型的最大似然估计流程分为两步迭代E-step计算后验概率γ(zₙₖ)πₖN(xₙ|μₖ,Σₖ)/∑ⱼπⱼN(xₙ|μⱼ,Σⱼ)实现软分配M-step更新模型参数μₖ更新为当前簇的加权均值Σₖ更新为加权协方差πₖ更新为簇占比初始化需指定各高斯分布的权重、均值和协方差。迭代过程中协方差矩阵可能从各向同性圆形发展为各向异性椭圆最终收敛至数据分布的最优表示。这张图片是对高斯混合模型Gaussian Mixture Model, GMM的简明总结内容涵盖其时间复杂度、优势与劣势具体如下二、时间复杂度O(t⋅k⋅n⋅d) O(t \cdot k \cdot n \cdot d)O(t⋅k⋅n⋅d)其中$ t $EM 算法的迭代次数$ k $预设的高斯成分簇数量$ n $数据点总数$ d $每个数据点的维度特征数→ 表明 GMM 的计算开销随维度、样本量和簇数线性增长对大规模高维数据可能较慢。三、主要优势提供不确定性估计soft assignment不同于 K-means 的硬聚类GMM 给出每个点属于各簇的概率分布如点 x 以 70% 概率属簇 A30% 属簇 B更适合建模模糊边界或重叠结构。对离群点outliers更鲁棒因基于概率密度建模单个异常点对整体参数影响相对较小相比基于距离的算法如 K-means 更不易被拉偏。四、主要劣势需预先指定簇数kkk无自动确定最优kkk的内置机制常依赖 BIC/AIC 或交叉验证。仅适用于近似凸形/椭球形簇因每个成分是多元高斯分布等概率面为椭球难以拟合弯曲、带状或非凸结构如环形、月牙形。存在奇异性singularity问题当某个高斯成分坍缩到单个点协方差矩阵行列式 → 0对数似然趋于无穷导致 EM 算法崩溃——需通过正则化、协方差约束或先验如 Bayesian GMM缓解。你这段代码其实是一个「从零实现 演示」版的 GMM高斯混合模型用 EM 算法做聚类并用 KMeans 做初始化。我们按模块拆开讲一步步把它“拆开看懂”。代码相关一、代码介绍代码做的事情可以概括成生成三簇二维高斯分布的数据generate_X用 GMM 对这些数据做聚类GMM.fitGMM.predict把聚类结果画出来不同颜色表示不同簇灰色点表示均值二、GMM 类的核心属性classGMM(object):def__init__(self,n_clusters,max_iter50):self.__Kn_clusters self.__max_itermax_iter self.__posterioriNone# 后验概率 γ_{k,n}self.__muNone# 每个高斯分量的均值 μ_kself.__covNone# 每个高斯分量的协方差 Σ_kself.__prioriNone# 每个高斯分量的先验概率 π_k__K簇的个数高斯分量个数__posteriori大小为(K, N)第k行第n列表示样本n属于簇k的后验概率 (\gamma_{k,n})__mu形状(K, D)每一行是一个簇的均值向量__cov形状(K, D, D)每个簇一个协方差矩阵__priori形状(K, 1)每个簇的混合系数 (\pi_k)且 (\sum_k \pi_k 1)三、fitEM 算法的主循环deffit(self,data):N,_data.shape# 1. 用 KMeans 初始化 GMM 参数self.__init_kmeans(data)# 2. EM 迭代foriinrange(self.__max_iter):# E 步更新后验概率 γ_{k,n}self.__get_expectation(data)# 有效样本数 N_k ∑_n γ_{k,n}effective_countnp.sum(self.__posteriori,axis1)# M 步更新 μ_kself.__munp.asarray([np.dot(self.__posteriori[k],data)/effective_count[k]forkinrange(self.__K)])# M 步更新 Σ_kself.__covnp.asarray([np.dot((data-self.__mu[k]).T,np.dot(np.diag(self.__posteriori[k].ravel()),data-self.__mu[k]))/effective_count[k]forkinrange(self.__K)])# M 步更新 π_kself.__priori(effective_count/N).reshape((self.__K,1))1初始化调用self.__init_kmeans(data)用 KMeans 的结果来初始化均值μ_k KMeans 的簇中心协方差Σ_k 每个簇内样本的协方差先验π_k 每个簇样本数 / 总样本数后验γ_{k,n}先设为 02E 步计算后验概率在循环中先调用__get_expectation(data)它会根据当前的μ_k, Σ_k, π_k计算每个样本属于每个簇的概率后面单独讲。3M 步更新参数有效样本数对应代码effective_countnp.sum(self.__posteriori,axis1)# 形状 (K,)更新均值对应代码self.__munp.asarray([np.dot(self.__posteriori[k],data)/effective_count[k]forkinrange(self.__K)])这里self.__posteriori[k]是形状(N,)data是(N, D)点乘后得到(D,)再除以N_k。更新协方差对应代码self.__covnp.asarray([np.dot((data-self.__mu[k]).T,np.dot(np.diag(self.__posteriori[k].ravel()),data-self.__mu[k]))/effective_count[k]forkinrange(self.__K)])解释一下data - self.__mu[k]形状(N, D)np.diag(self.__posteriori[k].ravel())形状(N, N)对角线上是 γ_{k,n}中间那一乘相当于对每个样本的偏差(x_n - μ_k)乘上权重 γ_{k,n}最后左乘转置(D, N)得到(D, D)的协方差矩阵更新先验概率对应代码self.__priori(effective_count/N).reshape((self.__K,1))四、predict用后验概率做分类defpredict(self,data):self.__get_expectation(data)resultnp.argmax(self.__posteriori,axis0)returnresult先重新计算一次后验概率γ_{k,n}对每个样本n找出后验概率最大的簇k作为该样本的类别标签self.__posteriori形状是(K, N)axis0表示对每一列每个样本取最大值所在的行索引五、初始化方法随机 KMeans1随机初始化代码里没用到但实现了def__init_random(self,data):N,_data.shape self.__posteriorinp.zeros((self.__K,N))self.__mudata[np.random.choice(np.arange(N),sizeself.__K,replaceFalse)]self.__covnp.asarray([np.cov(data,rowvarFalse)]*self.__K)self.__priorinp.ones((self.__K,1))/self.__K随机选K个样本作为初始均值协方差全部用整体数据的协方差先验概率均匀分布2KMeans 初始化实际使用的def__init_kmeans(self,data):N,_data.shape k_meansKMeans(initk-means,n_clustersself.__K)k_means.fit(data)categoryk_means.labels_ self.__posteriorinp.zeros((self.__K,N))self.__muk_means.cluster_centers_ self.__covnp.asarray([np.cov(data[categoryk],rowvarFalse)forkinrange(self.__K)])value_countspd.Series(category).value_counts()self.__priorinp.asarray([value_counts[k]/Nforkinrange(self.__K)]).reshape((self.__K,1))用 KMeans 的结果来均值直接用簇中心协方差对每个簇内的数据单独算协方差先验每个簇的样本数 / 总样本数这样初始化通常比完全随机更稳定EM 更容易收敛到“好一点”的解。六、E 步计算后验概率__get_expectationdef__get_expectation(self,data):forkinrange(self.__K):self.__posteriori[k]multivariate_normal.pdf(data,meanself.__mu[k],covself.__cov[k])self.__posteriorinp.dot(np.diag(self.__priori.ravel()),self.__posteriori)self.__posteriori/np.sum(self.__posteriori,axis0)这一步对应 EM 算法中的 E-step计算似然( p(x_n | z_n k) )multivariate_normal.pdf(data,meanself.__mu[k],covself.__cov[k])对所有样本data计算在第k个高斯分布下的概率密度得到形状(N,)的数组存到self.__posteriori[k]中。乘上先验( \pi_k )self.__posteriorinp.dot(np.diag(self.__priori.ravel()),self.__posteriori)这里相当于对每一行每个簇乘上对应的π_k。归一化得到后验概率对应代码self.__posteriori/np.sum(self.__posteriori,axis0)对每一列每个样本做归一化使得对所有簇的概率和为 1。七、数据生成函数generate_Xdefgenerate_X(true_Mu,true_Var):num1,mu1,var1400,true_Mu[0],true_Var[0]X1np.random.multivariate_normal(mu1,np.diag(var1),num1)num2,mu2,var2600,true_Mu[1],true_Var[1]X2np.random.multivariate_normal(mu2,np.diag(var2),num2)num3,mu3,var31000,true_Mu[2],true_Var[2]X3np.random.multivariate_normal(mu3,np.diag(var3),num3)Xnp.vstack((X1,X2,X3))returnX一共三簇数据每簇是一个二维高斯分布均值true_Mu[i]协方差用np.diag(true_Var[i])构造对角矩阵也就是两个维度之间不相关每簇样本数分别是 400、600、1000最后用vstack拼成一个(2000, 2)的数据集八、主程序训练 可视化if__name____main__:K3true_Mu[[0.5,0.5],[5.5,2.5],[1,7]]true_Var[[1,3],[2,2],[6,2]]Xgenerate_X(true_Mu,true_Var)gmmGMM(n_clustersK)gmm.fit(X)categorygmm.predict(X)color[red,blue,green,cyan,magenta]labels[fCluster{k:02d}forkinrange(K)]forkinrange(K):plt.scatter(X[categoryk][:,0],X[categoryk][:,1],ccolor[k],labellabels[k])mugmm.get_mu()plt.scatter(mu[:,0],mu[:,1],s300,cgrey,markerP,labelCentroids)plt.xlabel(X)plt.ylabel(Y)plt.legend()plt.title(GMM Testcase)plt.show()流程生成数据三簇二维高斯创建 GMM 模型GMM(n_clusters3)训练gmm.fit(X)执行 EM 算法预测类别category gmm.predict(X)画图不同簇用不同颜色散点用gmm.get_mu()得到每个簇的均值用灰色大点标出来