购物网站建站系统,课堂资源管理网站开发需求分析,让互联网之光点亮生活,代理app推广高斯混合模型GMM及EM迭代求解算法#xff08;含代码实现#xff09;
高斯分布与高斯混合模型
高斯分布
高斯分布大家都很熟悉了#xff0c;下面是一元高斯分布的概率密度函数#xff08;Probability Density Function#xff0c;PDF#xff09;#xff1a; P(x)N(μ,…高斯混合模型GMM及EM迭代求解算法含代码实现
高斯分布与高斯混合模型
高斯分布
高斯分布大家都很熟悉了下面是一元高斯分布的概率密度函数Probability Density FunctionPDF P(x)N(μ,σ2)12πσexp(−(x−μ)22σ2)P(x)N(\mu,\sigma^2)\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{(x-\mu)^2}{2\sigma^2}) P(x)N(μ,σ2)2πσ1exp(−2σ2(x−μ)2) 其中 μ\muμ 和 σ2\sigma^2σ2 分别是该高斯分布的均值和方差而如果是多元高斯分布则为 P(x)1(2π)D2∣Σ∣12exp(−(x−μ)TΣ−1(x−μ)2)P(x)\frac{1}{(2\pi)^{\frac{D}{2}}|\Sigma|^{\frac{1}{2}}}\exp(-\frac{(x-\mu)^T\Sigma^{-1}(x-\mu)}{2}) P(x)(2π)2D∣Σ∣211exp(−2(x−μ)TΣ−1(x−μ)) 其中 μ\muμ 、Σ\SigmaΣ 和 DDD 分别是均值、协方差矩阵和数据维度。
高斯混合模型
高斯混合模型顾名思义就是由多个单高斯模型组合而成的模型。具体来说它是具有如下形式的概率分布模型 P(x∣θ)∑k1Kαkϕ(x∣θk)P(x|\theta)\sum_{k1}^K\alpha_k\phi(x|\theta_k) P(x∣θ)k1∑Kαkϕ(x∣θk) 其中
KKK 是该混合模型中单高斯模型的个数αk\alpha_kαk 是各个单高斯模型的权重系数它满足αk≥0\alpha_k\ge0αk≥0 且 ∑k1Kαk1\sum_{k1}^K\alpha_k1∑k1Kαk1ϕ(x∣θk)\phi(x|\theta_k)ϕ(x∣θk) 是第 kkk 个单高斯模型的概率密度也称为 kkk 个分模型θk\theta_kθk 表示第 kkk 个单高斯模型的参数均值和方差即 θk(μk,σk2)\theta_k(\mu_k,\sigma_k^2)θk(μk,σk2)
为什么是高斯
对于一般地混合模型我们可以选用任何自己认为合适的单概率分布模型。这
里使用高斯混合模型是因为由中心极限定理可以得知将现实世界中的数据分布假设为高斯分布是比较合理的并且高斯分布具备很好的数学性质以及良好的计算性能。
为什么要混合
现实世界中我们的样本通常都是会有多个特征来描述此时如果使用单个高斯模型来对问题进行建模显然表达能力是不足的。而如果我们使用多个高斯模型按照一定的权重参数将它们组合将大大提高整个模型的表达能力。
具体来说我们知道单个高斯模型是单峰的如图所示即只有某一个区间的概率特别高两端都是逐渐降低的。而如果我们使用两个高斯模型混合得到的模型表达能力更强可以处理现实世界中的复杂分布。
理论上混合高斯模型的概率密度函数曲线可以是任意形状的非线性函数。 图1单高斯模型与混合高斯模型的概率密度函数 高斯混合模型的三种理解角度几何角度 从几何的角度来看高斯混合模型在图像上如图1所示它由多个高斯分布叠加而成是多个高斯分布的加权平均值。 P(x∣θ)∑k1Kαkϕ(x∣θk)P(x|\theta)\sum_{k1}^K\alpha_k\phi(x|\theta_k) P(x∣θ)k1∑Kαkϕ(x∣θk) 在这个角度下这里的 αk\alpha_kαk 就是每个单高斯模型的权重。 混合模型的角度 高斯混合模型中有两个变量 X(x1,x2,…,xN)X(x_1,x_2,\dots,x_N)X(x1,x2,…,xN)观测数据Z(z1,z2,…,zN)Z(z_1,z_2,\dots,z_N)Z(z1,z2,…,zN)隐变量 这两个变量合称为完全数据。 关于什么是隐变量这是笔者从网络上找到的一个例子感觉这个例子和 GMM 中的引入的隐变量意思非常接近可结合下面样本生成的角度来理解并且比较直观 举个例子吧 一个人拿着n个袋子里面有m种颜色不同的球。现在这个人随机地抓球规则如下 先随机挑一个袋子 从这个袋子中随机挑一个球 如果你站在这个人旁边你目睹了整个过程这个人选了哪个袋子、抓出来的球是什么颜色的。然后你把每次选择的袋子和抓出来的球的颜色都记录下来样本观察值那个人不停地抓你不停地记。最终你就可以通过你的记录推测出每个袋子里每种球颜色的大致比例。并且你记录的越多推测的就越准中心极限定理。 然而抓球的人觉得这样很不爽于是决定不告诉你他从哪个袋子里抓的球只告诉你抓出来的球的颜色是什么。这时候“选袋子”的过程由于你看不见其实就相当于是一个隐变量。 隐变量在很多地方都是能够出现的。现在我们经常说的隐变量主要强调它的“latent”。所以广义上的隐变量主要就是指“不能被直接观察到但是对系统的状态和能观察到的输出存在影响的一种东西”。 在 GMM 中 ZZZ 就是我们引入一个隐变量 zi,i[1,2,…,N]z_i,i[1,2,\dots,N]zi,i[1,2,…,N] 用它来表示对应的观测样本 xix_ixi 是属于哪一个高斯分布具体来说是样本 xix_ixi 属于每一个高斯分布的概率类似软分类的思想。 这样很自然地ZZZ 应当是一个离散随机变量并服从多项分布。其分布律为 ZZZC1C_1C1C2C_2C2…\dots…CKC_KCKP(z)P(z)P(z)p1p_1p1p2p_2p2…\dots…pKp_KpK其中 CkC_kCk 表示第 kkk 个单高斯模型 ∑k1Kpk1\sum_{k1}^Kp_k1∑k1Kpk1 。 现在我们推导在混合模型的角度下高斯混合模型的概率分布 P(X)∫zP(X,Z)dz∑ZP(X,ZCk)∑k1KP(ZCk)P(X∣ZCk)∑k1KpkN(X∣μk,σk2)∑k1Kpkϕ(x∣θk)\begin{align} P(X)\int_zP(X,Z)dz\\ \sum_ZP(X,ZC_k)\\ \sum_{k1}^KP(ZC_k)P(X|ZC_k)\\ \sum_{k1}^Kp_kN(X|\mu_k,\sigma_k^2)\\ \sum_{k1}^Kp_k\phi(x|\theta_k) \end{align} P(X)∫zP(X,Z)dzZ∑P(X,ZCk)k1∑KP(ZCk)P(X∣ZCk)k1∑KpkN(X∣μk,σk2)k1∑Kpkϕ(x∣θk) 推导过程解释由联合概率密度函数 P(X,Z)P(X,Z)P(X,Z) 求边缘概率密度函数 P(X)P(X)P(X)最直接的想法就是将另一个随机变量 ZZZ 直接积掉又由于这里的 ZZZ 是一个离散的随机变量因此应该是积分应该写为求和然后由公式 P(X,Y)P(Y)P(X∣Y)P(X,Y)P(Y)P(X|Y)P(X,Y)P(Y)P(X∣Y) 进行分解此时前一项 P(ZCk)P(ZC_k)P(ZCk) 由分布律知就是 pkp_kpk 而后一项条件概率 P(X∣ZCk)P(X|ZC_k)P(X∣ZCk) 表示的是在选定第 kkk 个高斯分布之后的概率密度函数那自然就是其本身的概率密度函数 N(X∣μk,σk2)N(X|\mu_k,\sigma_k^2)N(X∣μk,σk2)。 这里的 pkp_kpk 表示的是多项分布 ZZZ 的概率取值这里的 pkp_kpk 其实对应的就是上一种角度中的权重系数 αk\alpha_kαk 但是在两种不同的理解角度中有不同的物理含义。 从混合模型的角度来看GMM是由离散的多项分布连续的高斯分布组成。 样本生成的角度 GMM 是一个生成模型从样本生成的角度来说 对于每一个样本点 xix_ixi 我们可以认为它是通过这样的过程得到的 先通过一个多项KKK 项分布选择一个单高斯模型相当于掷一个有 KKK 个面的骰子然后从对应的高斯分布中进行采样得到 xix_ixi。 为什么MLE无法求出高斯混合模型的解析解 高斯混合模型的参数 高斯混合模型中的参数就是我们公式中的 θ\thetaθ 是哪些呢很明显的参数包括每个单高斯模型对应的权重系数 αk\alpha_kαk 和每个单高斯模型自身的均值方差 θk(μk,σk2)\theta_k(\mu_k,\sigma_k^2)θk(μk,σk2) 。 即混合高斯模型中的参数 θ(α1,α2,…,αK;θ1,θ2,…,θK)(α1,α2,…,αK;μ1,μ2,…,μK;σ1,σ2,…,σK)\theta(\alpha_1,\alpha_2,\dots,\alpha_K;\theta_1,\theta_2,\dots,\theta_K)(\alpha_1,\alpha_2,\dots,\alpha_K;\mu_1,\mu_2,\dots,\mu_K;\sigma_1,\sigma_2,\dots,\sigma_K) θ(α1,α2,…,αK;θ1,θ2,…,θK)(α1,α2,…,αK;μ1,μ2,…,μK;σ1,σ2,…,σK) 观测数据 x1,x2,…,xnx_1,x_2,\dots,x_nx1,x2,…,xn 由参数为 θ\thetaθ 的高斯混合模型生成 高斯混合模型的MLE尝试 对于单高斯模型我们可以直接用极大似然估计 MLE 来求解。下面我们对混合高斯模型也进行类似的尝试看能不能行得通 θ^MLEargmaxθP(X)argmaxθlogP(X)argmaxθlogP(x1)P(x2)…P(XN)argmaxθlog∏i1NP(xi)argmaxθ∑i1NlogP(xi)argmaxθ∑i1Nlog∑k1Kpkϕ(x∣θk)\begin{align} \hat{\theta}_{MLE}\arg\max_\theta P(X)\\ \arg\max_\theta \log P(X)\\ \arg\max_\theta\log P(x_1)P(x_2)\dots P(X_N)\arg\max_\theta\log \prod_{i1}^NP(x_i)\\ \arg\max_\theta\sum_{i1}^N\log P(x_i)\\ \arg\max_\theta\sum_{i1}^N\log \sum_{k1}^Kp_k\phi(x|\theta_k) \end{align} θ^MLEargθmaxP(X)argθmaxlogP(X)argθmaxlogP(x1)P(x2)…P(XN)argθmaxlogi1∏NP(xi)argθmaxi1∑NlogP(xi)argθmaxi1∑Nlogk1∑Kpkϕ(x∣θk) 至此下一步应该是对 θk\theta_kθk 各个值进行求偏导然后令导数为 0。 但是我们发现在对数 log\loglog 里面还有一个连加符号。对于对数里的连乘我们可以利用对数的性质直接拿出来变为连加但是对于对数里本身的连加符号是非常难处理的。因此到这一步之后无法直接求出混合高斯模型的解析解。 接下来我们将介绍如何利用 EM 算法来迭代求解高斯混合模型的参数 θ\thetaθ 。 EM算法迭代求解高斯混合模型 EM算法的一般步骤 EM 算法是一种非监督期望最大化算法。其结合了极大似然和迭代求解的方法去预估数据的分布。EM 算法主要用来解决含有隐变量的混合模型的参数估计非常适合求解高斯混合模型。 这里直接给出一般的算法步骤详情见[TODO] EM 算法通过迭代求 L(θ)logP(X∣θ)L(\theta)\log P(X|\theta)L(θ)logP(X∣θ) 的极大似然估计每次迭代包含两部E步求期望M步求极大化。 算法流程 输入观测变量数据 XXX隐变量数据 ZZZ联合分布 P(X,Z∣θ)P(X,Z|\theta)P(X,Z∣θ) 条件分布P(Z∣X,θ)P(Z|X,\theta)P(Z∣X,θ) 输出模型参数 θ\thetaθ 步骤 选择参数的初值 θ0\theta^{0}θ0 开始迭代 E 步记 θt\theta^{t}θt 为第 ttt 次迭代参数 θ\thetaθ 的估计值在第 i1i1i1 次迭代的 E 步计算 Q(θ,θt)EZ[logP(X,Z∣θ)∣X,θ(t)]∑ZlogP(X,Z∣θ)P(Z∣X,θ(t))\begin{align} Q(\theta,\theta^{t})\mathbb{E}_Z[\log P(X,Z|\theta)|X,\theta^{(t)}]\\ \sum_{Z}\log P(X,Z|\theta)P(Z|X,\theta^{(t)}) \end{align} Q(θ,θt)EZ[logP(X,Z∣θ)∣X,θ(t)]Z∑logP(X,Z∣θ)P(Z∣X,θ(t)) 这里 P(Z∣X,θ(t))P(Z|X,\theta^{(t)})P(Z∣X,θ(t)) 为给定观测数据 XXX 和当前参数估计 θ(t)\theta^{(t)}θ(t) 下隐变量数据 ZZZ 的条件概率分布 M 步求使 Q(θ,θ(t))Q(\theta,\theta^{(t)})Q(θ,θ(t)) 极大化的 θ\thetaθ 确定第 t1t1t1 次迭代的参数估计值 θt1\theta^{t1}θt1 : θ(t1)argmaxθQ(θ,θt)\theta^{(t1)}\arg\max_{\theta}Q(\theta,\theta^{t}) θ(t1)argθmaxQ(θ,θt) 重复 2、3 两步直到收敛。 函数 Q(θ,θ(t))Q(\theta,\theta^{(t)})Q(θ,θ(t)) 是 EM 算法的核心称为 QQQ 函数。 使用EM算法求解GMM推导 1 写出完全数据的对数似然函数 这里隐变量的定义与之前不同是因为 回忆我们之前介绍的生成模型的角度观测数据 xi,i1,2,…,Nx_i,i1,2,\dots,Nxi,i1,2,…,N 是这样产生的首先根据概率 αk\alpha_kαk 选择第 kkk 个单高斯分布模型 ϕ(x∣θk)\phi(x|\theta_k)ϕ(x∣θk) 然后根据 ϕ(x∣θk)\phi(x|\theta_k)ϕ(x∣θk) 生成观测数据 xix_ixi 。注意观测数据 xix_ixi 是已知的反映观测数据 xix_ixi 来自第 kkk 个单高斯模型的的数据是未知的用隐变量 zikz_{ik}zik 表示其定义如下 zik{1,第i个观测来自第k个分模型0,否则i1,2,…,N;k1,2,…,Kz_{ik}\begin{cases}1,\ \ \ 第i个观测来自第k个分模型 \\ 0,\ \ \ 否则\end{cases}\\ i1,2,\dots,N;\ \ \ k1,2,\dots,K zik{1, 第i个观测来自第k个分模型0, 否则i1,2,…,N; k1,2,…,K 隐变量 zikz_{ik}zik 与观测数据 XXX 共同组成了完全数据记为 (xi,zi1,zi2,…,ziK),i1,2,…,N(x_i,z_{i1},z_{i2},\dots,z_{iK}),\ \ \ i1,2,\dots,N(xi,zi1,zi2,…,ziK), i1,2,…,N 。 记 GMM 模型 P(x∣θ)∑k1Kαkϕ(x∣θk)P(x|\theta)\sum_{k1}^K\alpha_k\phi(x|\theta_k) P(x∣θ)k1∑Kαkϕ(x∣θk) 则有完全数据的似然函数 P(x,z∣θ)∏i1NP(xi,zi1,zi2,…,ziK∣θ)∏k1K∏i1N[αkϕ(xi∣θk)]zik∏k1Kαknk∏i1N[ϕ(xi∣θk)]zik∏k1Kαknk∏i1N[12πσkexp(−(xi−μk)22σk2)]zik\begin{align} P(x,z|\theta)\prod_{i1}^NP(x_i,z_{i1},z_{i2},\dots,z_{iK}|\theta)\\ \prod_{k1}^K\prod_{i1}^N[\alpha_k\phi(x_i|\theta_k)]^{z_{ik}}\\ \prod_{k1}^K\alpha_k^{n_k}\prod_{i1}^N[\phi(x_i|\theta_k)]^{z_{ik}}\\ \prod_{k1}^K\alpha_k^{n_k}\prod_{i1}^N[\frac{1}{\sqrt{2\pi}\sigma_k}\exp(-\frac{(x_i-\mu_k)^2}{2\sigma_k^2})]^{z_{ik}} \end{align} P(x,z∣θ)i1∏NP(xi,zi1,zi2,…,ziK∣θ)k1∏Ki1∏N[αkϕ(xi∣θk)]zikk1∏Kαknki1∏N[ϕ(xi∣θk)]zikk1∏Kαknki1∏N[2πσk1exp(−2σk2(xi−μk)2)]zik 其中nl∑i1Nzikn_l\sum_{i1}^Nz_{ik}nl∑i1Nzik ∑k1KnkN\sum_{k1}^Kn_kN∑k1KnkN。 那么完全数据的对输入似然函数为 logP(x,z∣θ)∑k1K{(nklogαk∑j1Nzik[log(12π)−logσk−12σl2(xi−μk)2]}\log P(x,z|\theta)\sum_{k1}^K\{(n_k\log\alpha_k\sum_{j1}^Nz_{ik}[\log(\frac{1}{\sqrt{2\pi}})-\log \sigma_k-\frac{1}{2\sigma_l^2}(x_i-\mu_k)^2]\} logP(x,z∣θ)k1∑K{(nklogαkj1∑Nzik[log(2π1)−logσk−2σl21(xi−μk)2]} 2 E步求期望写出 QQQ 函数 Q(θ,θ(t)E(logP(x,z∣θ)∣x,θ(t))E{∑k1K{(nklogαk∑j1Nzik[log(12π)−logσk−12σl2(xi−μk)2]}}∑k1K{∑i1N(Ezik)logαk∑i1N(Ezik[log12π−logσk−12σk2)(xi−μk)2]}\begin{align} Q(\theta,\theta^{(t)}\mathbb{E}(\log P(x,z|\theta)|x,\theta^{(t)})\\ \mathbb{E}\{\sum_{k1}^K\{(n_k\log\alpha_k\sum_{j1}^Nz_{ik}[\log(\frac{1}{\sqrt{2\pi}})-\log \sigma_k-\frac{1}{2\sigma_l^2}(x_i-\mu_k)^2]\}\}\\ \sum_{k1}^K\{\sum_{i1}^N(\mathbb{E}z_{ik})\log \alpha_k\sum_{i1}^N(\mathbb{E}z_{ik}[\log \frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{1}{2\sigma_k^2})(x_i-\mu_k)^2]\} \end{align} Q(θ,θ(t)E(logP(x,z∣θ)∣x,θ(t))E{k1∑K{(nklogαkj1∑Nzik[log(2π1)−logσk−2σl21(xi−μk)2]}}k1∑K{i1∑N(Ezik)logαki1∑N(Ezik[log2π1−logσk−2σk21)(xi−μk)2]} 这里需要计算 E(zik∣x,θ)\mathbb{E}(z_{ik}|x,\theta)E(zik∣x,θ)记为 z^ik\hat{z}_{ik}z^ik z^ikE(zik∣x,θ)P(zik1∣x,θ)P(zik1,xi∣θ)∑k1KP(zik1,xi∣θ)P(xi∣zik1,θ)P(zik1∣θ)∑k1KP(xi∣zik1,θ)P(zik1,θ)αkϕ(xi∣θk)∑k1Kαkϕ(xi∣θk),i1,2,…,N;k1,2,…,K\begin{align} \hat{z}_{ik}\mathbb{E}(z_{ik}|x,\theta)P(z_{ik}1|x,\theta)\\ \frac{P(z_{ik}1,x_i|\theta)}{\sum_{k1}^KP(z_{ik}1,x_i|\theta)}\\ \frac{P(x_i|z_{ik}1,\theta)P(z_{ik}1|\theta)}{\sum_{k1}^KP(x_i|z_{ik}1,\theta)P(z_{ik}1,\theta)}\\ \frac{\alpha_k\phi(x_i|\theta_k)}{\sum_{k1}^K\alpha_k\phi(x_i|\theta_k)},\ \ i1,2,\dots,N;\ k1,2,\dots,K \end{align} z^ikE(zik∣x,θ)P(zik1∣x,θ)∑k1KP(zik1,xi∣θ)P(zik1,xi∣θ)∑k1KP(xi∣zik1,θ)P(zik1,θ)P(xi∣zik1,θ)P(zik1∣θ)∑k1Kαkϕ(xi∣θk)αkϕ(xi∣θk), i1,2,…,N; k1,2,…,K z^ik\hat{z}_{ik}z^ik 是当前模型参数下第 iii 个观测数据来自第 kkk 个模型的概率称为分模型 kkk 对观测数据 xix_ixi 的响应程度。 将 z^ikEzik\hat{z}_{ik}\mathbb{E}z_{ik}z^ikEzik 和 nk∑i1NEzikn_k\sum_{i1}^N\mathbb{E}z_{ik}nk∑i1NEzik 代入得 Q(θ,θ(t))∑k1K{∑i1N(nklogαk∑i1N(z^ik[log12π−logσk−12σk2)(xi−μk)2]}Q(\theta,\theta^{(t)})\sum_{k1}^K\{\sum_{i1}^N(n_k\log \alpha_k\sum_{i1}^N(\hat{z}_{ik}[\log \frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{1}{2\sigma_k^2})(x_i-\mu_k)^2]\} Q(θ,θ(t))k1∑K{i1∑N(nklogαki1∑N(z^ik[log2π1−logσk−2σk21)(xi−μk)2]} 3 M步求极大 M 步是求函数 Q(θ,θ(t))Q(\theta,\theta^{(t)})Q(θ,θ(t)) 取得极大值时的 θ\thetaθ 作为 θ(t1)\theta^{(t1)}θ(t1) θ(t1)argmaxθQ(θ,θ(t))\theta^{(t1)}\arg\max_{\theta}Q(\theta,\theta^{(t)}) θ(t1)argθmaxQ(θ,θ(t)) 我们之前介绍过GMM 模型要估计的参数为 θ(α1,α2,…,αK;θ1,θ2,…,θK)(α1,α2,…,αK;μ1,μ2,…,μK;σ1,σ2,…,σK)\theta(\alpha_1,\alpha_2,\dots,\alpha_K;\theta_1,\theta_2,\dots,\theta_K)(\alpha_1,\alpha_2,\dots,\alpha_K;\mu_1,\mu_2,\dots,\mu_K;\sigma_1,\sigma_2,\dots,\sigma_K) θ(α1,α2,…,αK;θ1,θ2,…,θK)(α1,α2,…,αK;μ1,μ2,…,μK;σ1,σ2,…,σK) 对于 μk\mu_kμk 和 σk2\sigma_k^2σk2 我们只需根据上式分别对他们求偏导并令其为 0 即可 而对于 α^k\hat{\alpha}_kα^k 还需要注意约束条件 ∑k1Kαk1\sum_{k1}^K\alpha_k1∑k1Kαk1 。 结果如下 μ^k∑i1Nz^ikxi∑i1Nz^ikσ^k2∑i1Nz^ik(xi−μk)2∑i1Nz^ikα^knkN∑i1Nz^ikN\hat{\mu}_k\frac{\sum_{i1}^N\hat{z}_{ik}x_i}{\sum_{i1}^N\hat{z}_{ik}}\\ \hat{\sigma}_{k}^2\frac{\sum_{i1}^N\hat{z}_{ik}(x_i-\mu_k)^2}{\sum_{i1}^N\hat{z}_{ik}}\\ \hat{\alpha}_k\frac{n_k}{N}\frac{\sum_{i1}^N\hat{z}_{ik}}{N} μ^k∑i1Nz^ik∑i1Nz^ikxiσ^k2∑i1Nz^ik∑i1Nz^ik(xi−μk)2α^kNnkN∑i1Nz^ik 重复上述 2、3 步直到收敛。 GMM模型的EM算法步骤 至此可以给出GMM模型的EM算法步骤 选择参数的初值 θ0\theta^{0}θ0 开始迭代 E 步记 θt\theta^{t}θt 为第 ttt 次迭代参数 θ\thetaθ 的估计值在第 i1i1i1 次迭代的 E 步计算 Q(θ,θ(t))∑k1K{∑i1N(nklogαk∑i1N(z^ik[log12π−logσk−12σk2)(xi−μk)2]}Q(\theta,\theta^{(t)})\sum_{k1}^K\{\sum_{i1}^N(n_k\log \alpha_k\sum_{i1}^N(\hat{z}_{ik}[\log \frac{1}{\sqrt{2\pi}}-\log\sigma_k-\frac{1}{2\sigma_k^2})(x_i-\mu_k)^2]\} Q(θ,θ(t))k1∑K{i1∑N(nklogαki1∑N(z^ik[log2π1−logσk−2σk21)(xi−μk)2]} 这里 P(Z∣X,θ(t))P(Z|X,\theta^{(t)})P(Z∣X,θ(t)) 为给定观测数据 XXX 和当前参数估计 θ(t)\theta^{(t)}θ(t) 下隐变量数据 ZZZ 的条件概率分布 M 步求使 Q(θ,θ(t))Q(\theta,\theta^{(t)})Q(θ,θ(t)) 极大化的 θ\thetaθ 确定第 t1t1t1 次迭代的参数估计值 θt1\theta^{t1}θt1 : μ^k∑i1Nz^ikxi∑i1Nz^ikσ^k2∑i1Nz^ik(xi−μk)2∑i1Nz^ikα^knkN∑i1Nz^ikN\hat{\mu}_k\frac{\sum_{i1}^N\hat{z}_{ik}x_i}{\sum_{i1}^N\hat{z}_{ik}}\\ \hat{\sigma}_{k}^2\frac{\sum_{i1}^N\hat{z}_{ik}(x_i-\mu_k)^2}{\sum_{i1}^N\hat{z}_{ik}}\\ \hat{\alpha}_k\frac{n_k}{N}\frac{\sum_{i1}^N\hat{z}_{ik}}{N} μ^k∑i1Nz^ik∑i1Nz^ikxiσ^k2∑i1Nz^ik∑i1Nz^ik(xi−μk)2α^kNnkN∑i1Nz^ik 重复 2、3 两步直到收敛。 代码 给出参考实现[代码][[https://github.com/wl-lei/upload/blob/master/homework/GMM.py] import numpy as np
import randomdef calc_prob(X, K, pMu, pSigma):N X.shape[0]D X.shape[1]Px np.zeros((N, K))for i in range(K):Xshift X-np.tile(pMu[i], (N, 1))lambda_flag np.e**(-5)conv pSigma[i]lambda_flag*np.eye(D)inv_pSigma np.linalg.inv(conv)tmp np.sum(np.dot(Xshift, inv_pSigma)*Xshift, axis1)coef (2*np.pi)**(-D/2)*np.sqrt(np.linalg.det(inv_pSigma))Px[:, i] coef*np.e**(-1/2*tmp)return Pxdef gmm(X, K): #用array来处理threshold np.e**(-15)N X.shape[0]D X.shape[1]rndp random.sample(np.arange(N).tolist(),K)centroids X[rndp,:]pMu centroidspPi np.zeros((1, K))pSigma np.zeros((K, D, D))dist np.tile(np.sum(X*X, axis1).reshape(N,1), (1, K))np.tile(np.sum(pMu*pMu, axis1), (N, 1))-2*np.dot(X, pMu.T)labels np.argmin(dist,axis1)for i in range(K):index labels iXk X[index,:]pPi[:,i] (Xk.shape[0])/NpSigma[i] np.cov(Xk.T)Loss -float(inf)while True:Px calc_prob(X, K, pMu, pSigma)pGamma Px*np.tile(pPi, (N, 1))pGamma pGamma/np.tile(np.sum(pGamma, axis1).reshape(N,1), (1, K))Nk np.sum(pGamma, axis0)pMu np.dot(np.dot(np.diag(1/Nk), pGamma.T), X)pPi Nk/Nfor i in range(K):Xshift X-np.tile(pMu[i], (N, 1))pSigma[i] np.dot(Xshift.T, np.dot(np.diag(pGamma[:, i]), Xshift))/Nk[i]L np.sum(np.log(np.dot(Px, pPi.T)), axis0)if L-Loss threshold:breakLoss Lreturn Px,pMu,pSigma,pPiif __name__ __main__: Data_list []with open(data.txt, r) as file:for line in file.readlines(): point [] point.append(float(line.split()[0])) point.append(float(line.split()[1])) Data_list.append(point) Data np.array(Data_list)Px,pMu,pSigma,pPi gmm(Data, 2)print(Px,pMu,pSigma,pPi)Ref 高斯混合模型GMM高斯混合模型白板推导隐变量是什么详解EM算法与混合高斯模型(Gaussian mixture model, GMM) 高斯混合模型(GMM)的两种详解及简化
