经典机器学习模型(九)EM算法在高斯混合模型中的应用

经典机器学习模型(九)EM算法在高斯混合模型中的应用

EM算法的推导可以参考:

经典机器学习模型(九)EM算法的推导

若随机变量X服从一个数学期望为 μ μ μ、方差为 σ 2 σ^2 σ2的正态分布,可以记为 N ( μ , σ 2 ) N(μ,σ2) N(μσ2)

在这里插入图片描述

如果我们知道一组样本服从高斯分布,那么根据极大似然估计,就很容得出参数 u , σ 2 u,σ^2 u,σ2的值,求解过程如下:

在这里插入图片描述

注:

上图用均值而不是平均值,这两者的差别在于:期望还可以称为均值,是对应于总体分布的,但是平均值一般指用样本计算的算术平均数。

高斯混合模型,既然是高斯混合,那么一定是要有高斯分布的,接着还有另一个分布和它混合在一起。

假如我们有一个标准正态分布,即均值为0,方差1,另一个高斯分布的均值为5,方差为1,我们按照1:1叠加合并,也就是分别对应的占比是0.5,就得到了下图所示的混合高斯分布。

在这里插入图片描述

假如,我们现在有K个高斯分布:
N ( u 1 , σ 1 2 ) , N ( u 2 , σ 2 2 ) , . . . , N ( u K , σ K 2 ) 对应的比例为: a 1 , a 2 , . . . , a K , 意味着我们给每个高斯分布设一个权重 比例之和为 1 ,即 a 1 + a 2 + . . . + a K = 1 N(u_1,\sigma_1^2),N(u_2,\sigma_2^2),...,N(u_K,\sigma_K^2) \\ 对应的比例为:a_1,a_2,...,a_K,意味着我们给每个高斯分布设一个权重 \\ 比例之和为1,即a_1+a_2+...+a_K =1 N(u1,σ12),N(u2,σ22),...,N(uK,σK2)对应的比例为:a1,a2,...,aK,意味着我们给每个高斯分布设一个权重比例之和为1,即a1+a2+...+aK=1
如果我们想要估计出混合高斯分布里面每一个分布的参数 u u u σ 2 \sigma^2 σ2,此时就麻烦了。因为我们只能观测到样本,并不知道哪一个样本属于哪一个高斯分布。

如果,我们知道每一个样本属于哪一个高斯分布,那么就很简单了,利用极大似然估计,就能直接求出每一个高斯分布所拥有的样本的均值和方差,因此关键点就在于找出每一个样本属于哪一个高斯分布

1 高斯混合模型公式推导

1.1 高斯混合模型E步推导

根据上文,我们现在并不知道对应的样本到底是哪个高斯分布,这就是对应的隐变量

这恰好就是EM算法能够实现的情况。

所以,我们就先准备高斯混合模型的E步,这一步要完成的是两个目的:

  • 一是补全信息,也就是下文中即将要出场的隐变量 γ j k \gamma_{jk} γjk
  • 二是确定 M 步的目标函数。

1.1.1 隐变量的表示

我们设 γ j k = { 1 , y j 来自第 k 个高斯分布 N ( u k , σ k 2 ) 0 , 其他 例如观测的样本来自于第 1 个高斯分布,可以表示为: γ j 1 = { 1 , 观测样本 y j 来自第 1 个高斯分布 N ( u 1 , σ 1 2 ) 0 , 其他 那么,此时我们有 y 1 , y 2 , . . . , y N 观测样本,那么隐变量就可以表示为: y 1 , γ 11 , γ 12 , γ 13 , . . . , γ 1 K y 2 , γ 21 , γ 22 , γ 23 , . . . , γ 2 K . . . . . . y N , γ N 1 , γ N 2 , γ N 3 , . . . , γ N K 对应的 y 1 , y 2 , . . . , y N 就是「观测量」,是已知的。而后面的「隐变量」 γ j i 自然是未知的。 此时对应的参数 θ = [ a 1 , u 1 , σ 1 2 a 2 , u 2 , σ 2 2 . . . . . . a K , u K , σ K 2 ] 我们设\gamma_{jk}=\begin{cases} 1, & y_j来自第k个高斯分布N(u_k,\sigma_k^2)\\ 0,& \text{其他} \end{cases} \\ 例如观测的样本来自于第1个高斯分布,可以表示为:\\ \gamma_{j1}=\begin{cases} 1, & 观测样本y_j来自第1个高斯分布N(u_1,\sigma_1^2)\\ 0,& \text{其他} \end{cases}\\ 那么,此时我们有y_1,y_2,...,y_N观测样本,那么隐变量就可以表示为:\\ y_1,\gamma_{11},\gamma_{12},\gamma_{13},...,\gamma_{1K} \\ y_2,\gamma_{21},\gamma_{22},\gamma_{23},...,\gamma_{2K} \\ ......\\ y_N,\gamma_{N1},\gamma_{N2},\gamma_{N3},...,\gamma_{NK} \\ 对应的y_1,y_2,...,y_N就是「观测量」,是已知的。而后面的「隐变量」\gamma_{ji}自然是未知的。\\ 此时对应的参数\theta=\\ [ a_{1},u_{1},\sigma_{1}^2 \\ a_{2},u_{2},\sigma_{2}^2 \\ ......\\ a_{K},u_{K},\sigma_{K}^2 ] 我们设γjk={1,0,yj来自第k个高斯分布N(uk,σk2)其他例如观测的样本来自于第1个高斯分布,可以表示为:γj1={1,0,观测样本yj来自第1个高斯分布N(u1,σ12)其他那么,此时我们有y1,y2,...,yN观测样本,那么隐变量就可以表示为:y1,γ11,γ12,γ13,...,γ1Ky2,γ21,γ22,γ23,...,γ2K......yN,γN1,γN2,γN3,...,γNK对应的y1,y2,...,yN就是「观测量」,是已知的。而后面的「隐变量」γji自然是未知的。此时对应的参数θ=[a1,u1,σ12a2,u2,σ22......aK,uK,σK2]

1.1.2 隐变量的求解

γ j k \gamma_{jk} γjk是在当前模型参数下第 j j j个观测数据来自第 k k k个分模型的概率,那么如何求出 γ j k \gamma_{jk} γjk呢?

我们可以写成期望的形式:
γ j k = E ( γ j k ∣ y j , θ i ) = 1 ∗ P ( γ j k = 1 ∣ y j , θ i ) + 0 ∗ P ( γ j k = 0 ∣ y j , θ i ) = P ( γ j k = 1 ∣ y j , θ i ) \gamma_{jk}=E(\gamma_{jk}|y_j,\theta^i)\\ =1*P(\gamma_{jk}=1|y_j,\theta^i) + 0*P(\gamma_{jk}=0|y_j,\theta^i)\\ =P(\gamma_{jk}=1|y_j,\theta^i) γjk=E(γjkyj,θi)=1P(γjk=1∣yj,θi)+0P(γjk=0∣yj,θi)=P(γjk=1∣yj,θi)
然后,我们借助贝叶斯公式:
γ j k = P ( γ j k = 1 ∣ y j , θ i ) = P ( y j ∣ γ j k = 1 , θ i ) P ( γ j k = 1 ∣ θ i ) P ( y j , θ i ) 再利用全概率公式 p ( A ) = ∑ p ( A ∣ B i ) p ( B i ) ,将分母展开 = P ( y j ∣ γ j k = 1 , θ i ) P ( γ j k = 1 ∣ θ i ) ∑ k = 1 K P ( y j ∣ γ j k = 1 , θ i ) P ( γ j k = 1 ∣ θ i ) 上式中 θ i 是上一轮的参数,是已知的。 观察上式, P ( y j ∣ γ j k = 1 , θ i ) 是指高斯分布参数 θ i 已知, 样本 y j 属于第 k 个分布,那么可以用高斯密度函数进行表示 ∅ ( y j ∣ u k i , σ k 2 ( i ) ) P ( γ j k = 1 ∣ θ i ) 表示已知上一轮的参数 θ i ,那么属于第 k 个分布的概率很容易得出为 a k ( i ) 因此,我们可以求得隐变量的值,如下式: γ j k = ∅ ( y j ∣ u k i , σ k 2 ( i ) ) a k ( i ) ∑ k = 1 K ∅ ( y j ∣ u k i , σ k 2 ( i ) ) a k ( i ) 我们令 n k = ∑ j = 1 N γ j k , 表示对所有的样本进行统计,如果属于第 k 个分布,那么 γ j k = 1 n k 的意思就是属于第 k 个分布的样本数 \gamma_{jk}=P(\gamma_{jk}=1|y_j,\theta^i)\\ =\frac{P(y_j|\gamma_{jk}=1,\theta^i)P(\gamma_{jk}=1|\theta^i)}{P(y_j,\theta^i)}\\ 再利用全概率公式p(A)=\sum p(A|B_i)p(B_i),将分母展开\\ =\frac{P(y_j|\gamma_{jk}=1,\theta^i)P(\gamma_{jk}=1|\theta^i)}{\sum_{k=1}^{K} P(y_j|\gamma_{jk}=1,\theta^i)P(\gamma_{jk}=1|\theta^i)}\\ 上式中\theta^i是上一轮的参数,是已知的。\\ 观察上式,P(y_j|\gamma_{jk}=1,\theta^i)是指高斯分布参数\theta^i已知,\\样本y_j属于第k个分布,那么可以用高斯密度函数进行表示\emptyset(y_j|u_k^{i},\sigma_k^{2(i)}) \\ P(\gamma_{jk}=1|\theta^i)表示已知上一轮的参数\theta^i,那么属于第k个分布的概率很容易得出为a_k^{(i)}\\ 因此,我们可以求得隐变量的值,如下式:\\ \gamma_{jk}=\frac{\emptyset(y_j|u_k^{i},\sigma_k^{2(i)})a_k^{(i)}}{\sum_{k=1}^{K}\emptyset(y_j|u_k^{i},\sigma_k^{2(i)})a_k^{(i)}} \\ 我们令n_k=\sum_{j=1}^{N}\gamma_{jk},表示对所有的样本进行统计,如果属于第k个分布,那么\gamma_{jk}=1 \\n_k的意思就是属于第k个分布的样本数 γjk=P(γjk=1∣yj,θi)=P(yj,θi)P(yjγjk=1,θi)P(γjk=1∣θi)再利用全概率公式p(A)=p(ABi)p(Bi),将分母展开=k=1KP(yjγjk=1,θi)P(γjk=1∣θi)P(yjγjk=1,θi)P(γjk=1∣θi)上式中θi是上一轮的参数,是已知的。观察上式,P(yjγjk=1,θi)是指高斯分布参数θi已知,样本yj属于第k个分布,那么可以用高斯密度函数进行表示(yjuki,σk2(i))P(γjk=1∣θi)表示已知上一轮的参数θi,那么属于第k个分布的概率很容易得出为ak(i)因此,我们可以求得隐变量的值,如下式:γjk=k=1K(yjuki,σk2(i))ak(i)(yjuki,σk2(i))ak(i)我们令nk=j=1Nγjk,表示对所有的样本进行统计,如果属于第k个分布,那么γjk=1nk的意思就是属于第k个分布的样本数
根据以上推导的公式,我们就能在已知上一轮参数 θ ( i ) \theta^{(i)} θ(i)的情况下,求得隐变量 γ j k \gamma_{jk} γjk的值。

注:

我们现在在已知上一轮参数 θ ( i ) \theta^{(i)} θ(i)的情况下,求出了隐变量 γ j k \gamma_{jk} γjk的值,其实我们现在已经知道哪一个样本属于哪一个高斯分布了,因此只要求得每一个高斯分布所有样本的均值和方差,就求出了下一轮的的参数 θ ( i + 1 ) \theta^{(i+1)} θ(i+1)
我们现在求出了隐变量 γ j k 的值 , 属于第 k 个高斯分布的样本数我们表示为 n k = ∑ j = 1 N γ j k 下一轮的 a k 是第 k 个高斯分布样本所占的比例,很容易得出: a k = n k N 下一轮的 u k 是第 k 个高斯分布的均值 , 很容易得出: u k = ∑ j = 1 N r j k y j n k 下一轮的 σ 2 是第 k 个高斯分布的方差 , 很容易得出: σ k 2 = ∑ j = 1 N r j k ( y j − u k ) 2 n k 我们现在求出了隐变量\gamma_{jk}的值,属于第k个高斯分布的样本数我们表示为n_k=\sum_{j=1}^{N}\gamma_{jk}\\ 下一轮的a_k是第k个高斯分布样本所占的比例,很容易得出:a_k=\frac{n_k}{N}\\ 下一轮的u_k是第k个高斯分布的均值,很容易得出:u_k=\frac{\sum_{j=1}^{N}r_{jk}y_j}{n_k}\\ 下一轮的\sigma^{2}是第k个高斯分布的方差,很容易得出:\sigma_k^{2}=\frac{\sum_{j=1}^{N}r_{jk}(y_j-u_k)^2}{n_k} 我们现在求出了隐变量γjk的值,属于第k个高斯分布的样本数我们表示为nk=j=1Nγjk下一轮的ak是第k个高斯分布样本所占的比例,很容易得出:ak=Nnk下一轮的uk是第k个高斯分布的均值,很容易得出:uk=nkj=1Nrjkyj下一轮的σ2是第k个高斯分布的方差,很容易得出:σk2=nkj=1Nrjk(yjuk)2
我们下面通过另一个角度,即通过联合概率分布构造目标函数 Q ( θ , θ ( i ) ) Q(\theta, \theta^{(i)}) Q(θ,θ(i))进行求解。

1.1.3 联合概率分布

此时的联合概率分布,因为是连续数据,所以我们可以用联合概率密度来表示,对应的值应该是连乘的形式:

在这里插入图片描述

然接着我们生成完全数据的对数似然函数
l n L ( θ ) = ∑ k = 1 K ( n k l n a k + ∑ j = 1 N l n [ ∅ ( y j ∣ θ k ) ] r j k ) = ∑ k = 1 K ( n k l n a k + ∑ j = 1 N r j k l n [ ∅ ( y j ∣ θ k ) ] ) 其中 n k = ∑ j = 1 N γ j k lnL(\theta)=\sum_{k=1}^{K}(n_klna_k+\sum_{j=1}^{N}ln[\emptyset(y_j|\theta_k)]^{r_{jk}})\\ =\sum_{k=1}^{K}(n_klna_k+\sum_{j=1}^{N}{r_{jk}}ln[\emptyset(y_j|\theta_k)])\\ 其中n_k=\sum_{j=1}^{N}\gamma_{jk} lnL(θ)=k=1K(nklnak+j=1Nln[(yjθk)]rjk)=k=1K(nklnak+j=1Nrjkln[(yjθk)])其中nk=j=1Nγjk
根据以上推导,我们从E步就得到了目标函数为:
Q ( θ , θ ( i ) ) = l n L ( θ ) = ∑ k = 1 K ( n k l n a k + ∑ j = 1 N r j k l n [ ∅ ( y j ∣ θ k ) ] ) 约束条件: ∑ k = 1 K a k = 1 其中 n k = ∑ j = 1 N γ j k , 表示属于第 k 个分布的样本数 Q(\theta, \theta^{(i)})=lnL(\theta)=\sum_{k=1}^{K}(n_klna_k + \sum_{j=1}^{N}r_{jk}ln[\emptyset(y_j|\theta_k)])\\ 约束条件:\sum_{k=1}^{K} a_k=1\\ 其中n_k=\sum_{j=1}^{N}\gamma_{jk},表示属于第k个分布的样本数 Q(θ,θ(i))=lnL(θ)=k=1K(nklnak+j=1Nrjkln[(yjθk)])约束条件:k=1Kak=1其中nk=j=1Nγjk,表示属于第k个分布的样本数

1.2 高斯混合模型的M步

通过EM算法的E步,我们得到了目标函数 Q ( θ , θ ( i ) ) Q(\theta, \theta^{(i)}) Q(θ,θ(i)),我们要最大化目标函数
θ ( i + 1 ) = a r g m a x Q ( θ , θ ( i ) ) s . t . ∑ k = 1 K a k = 1 求解有约束问题的最优化问题,首先利用拉格朗日乘子法转换为无约束条件 Q 0 = a r g m a x [ Q ( θ , θ ( i ) ) + λ ( ∑ k = 1 K a k − 1 ) ] Q 0 = ∑ k = 1 K ( n k l n a k + ∑ j = 1 N r j k l n [ ∅ ( y j ∣ θ k ) ] ) + λ ( ∑ k = 1 K a k − 1 ) 将高斯密度函数代入,得到 = ∑ k = 1 K ( n k l n a k + ∑ j = 1 N r j k ( − 1 2 l n 2 π − 1 2 l n σ k 2 − 1 2 σ k 2 ( y j − u k 2 ) ) ) + λ ( ∑ k = 1 K a k − 1 ) \theta^{(i+1)}=argmax \quad Q(\theta, \theta^{(i)}) \\ s.t. \quad \sum_{k=1}^{K} a_k=1\\ 求解有约束问题的最优化问题,首先利用拉格朗日乘子法转换为无约束条件\\ Q_0=argmax [Q(\theta, \theta^{(i)})+\lambda(\sum_{k=1}^{K} a_k-1)]\\ Q_0=\sum_{k=1}^{K}(n_klna_k+\sum_{j=1}^{N}{r_{jk}}ln[\emptyset(y_j|\theta_k)])+\lambda(\sum_{k=1}^{K} a_k-1) \\ 将高斯密度函数代入,得到\\ =\sum_{k=1}^{K}(n_klna_k+\sum_{j=1}^{N}{r_{jk}}(-\frac{1}{2}ln2π-\frac{1}{2}ln \sigma_k^2-\frac{1}{2\sigma_k^2}(y_j-u_k^2)))+\lambda(\sum_{k=1}^{K} a_k-1) θ(i+1)=argmaxQ(θ,θ(i))s.t.k=1Kak=1求解有约束问题的最优化问题,首先利用拉格朗日乘子法转换为无约束条件Q0=argmax[Q(θ,θ(i))+λ(k=1Kak1)]Q0=k=1K(nklnak+j=1Nrjkln[(yjθk)])+λ(k=1Kak1)将高斯密度函数代入,得到=k=1K(nklnak+j=1Nrjk(21ln2π21lnσk22σk21(yjuk2)))+λ(k=1Kak1)
然后,我们分别对几个参数求偏导,并令其为0。
∂ Q 0 ∂ a k = n k a k + λ = 0 怎么求 λ 呢 ? 我们可以利用约束条件 ∑ k = 1 K a k = 1 。 易知 n k = − λ a k , 那么 ∑ k = 1 K n k = − λ ∑ k = 1 K a k 左边式子就是总样本数 N ,右边就变成了 − λ ,因此 − λ = N 因此 a k = n k − λ = n k N ,和 1.2 中注下我们计算的一样。 \frac{\partial Q_0}{\partial a_k}=\frac{n_k}{a_k}+\lambda=0\\ 怎么求\lambda呢?我们可以利用约束条件\sum_{k=1}^{K} a_k=1。\\ 易知n_k=-\lambda a_k,\\ 那么\sum_{k=1}^{K} n_k=-\lambda \sum_{k=1}^{K} a_k\\ 左边式子就是总样本数N,右边就变成了-\lambda,因此-\lambda=N\\ 因此a_k=\frac{n_k}{-\lambda}=\frac{n_k}{N},和1.2中注下我们计算的一样。 akQ0=aknk+λ=0怎么求λ?我们可以利用约束条件k=1Kak=1易知nk=λak那么k=1Knk=λk=1Kak左边式子就是总样本数N,右边就变成了λ,因此λ=N因此ak=λnk=Nnk,和1.2中注下我们计算的一样。
同理,我们对 u k u_k uk求偏导:
∂ Q 0 ∂ u k = ∑ j = 1 N r j k ( + 2 2 σ k 2 ( y j − u k ) ) = 0 方差 σ 2 > 0 , 因此 ∑ j = 1 N r j k ( y j − u k ) = 0 因此 u k = ∑ j = 1 N r j k y j ∑ j = 1 N r j k = ∑ j = 1 N r j k y j n k ,和 1.2 中注下我们计算的一样。 \frac{\partial Q_0}{\partial u_k}=\sum_{j=1}^{N}{r_{jk}}(+\frac{2}{2\sigma_k^2}(y_j-u_k))=0\\ 方差\sigma^2>0,因此\sum_{j=1}^{N}{r_{jk}}(y_j-u_k)=0\\ 因此u_k=\frac{\sum_{j=1}^{N}r_{jk}y_j}{\sum_{j=1}^{N}r_{jk}}=\frac{\sum_{j=1}^{N}r_{jk}y_j}{n_k},和1.2中注下我们计算的一样。 ukQ0=j=1Nrjk(+2σk22(yjuk))=0方差σ2>0,因此j=1Nrjk(yjuk)=0因此uk=j=1Nrjkj=1Nrjkyj=nkj=1Nrjkyj,和1.2中注下我们计算的一样。
同理,我们对 σ k 2 \sigma_k^2 σk2求偏导:
∂ Q 0 ∂ σ k 2 = ∑ j = 1 N r j k ( − 1 2 σ K 2 + 1 2 σ k 4 ( y j − u k ) 2 ) = 0 同理,方差 σ k 2 > 0 , ∑ j = 1 N r j k ( − σ k 2 + ( y j − u k ) 2 ) = 0 因此 σ k 2 = ∑ j = 1 N r j k ( y j − u k ) 2 ∑ j = 1 N r j k = ∑ j = 1 N r j k ( y j − u k ) 2 n k ,和 1.2 中注下我们计算的一样。 \frac{\partial Q_0}{\partial \sigma_k^2}=\sum_{j=1}^{N}{r_{jk}}(-\frac{1}{2\sigma_K^2}+\frac{1}{2\sigma_k^4}(y_j-u_k)^2)=0\\ 同理,方差\sigma_k^2>0,\sum_{j=1}^{N}{r_{jk}}(-\sigma_k^2+(y_j-u_k)^2)=0 \\ 因此\sigma_k^2=\frac{\sum_{j=1}^{N}r_{jk}(y_j-u_k)^2}{\sum_{j=1}^{N}r_{jk}}=\frac{\sum_{j=1}^{N}r_{jk}(y_j-u_k)^2}{n_k},和1.2中注下我们计算的一样。 σk2Q0=j=1Nrjk(2σK21+2σk41(yjuk)2)=0同理,方差σk2>0,j=1Nrjk(σk2+(yjuk)2)=0因此σk2=j=1Nrjkj=1Nrjk(yjuk)2=nkj=1Nrjk(yjuk)2,和1.2中注下我们计算的一样。

1.3 高斯混合模型举例

我们根据上面公式的推导,现在应该很容易理解下面EM算法求解高斯混合模型参数的过程了。

在这里插入图片描述

我们通过一个例子,再重复以下上述求解过程:

假设有6个未知分布的数据 x 1 , x 2 , . . . , x 6 {x1,x2,...,x6} x1,x2,...,x6 ,采用2个高斯分布线性组合对其进行拟合。

在这里插入图片描述
在这里插入图片描述

2 高斯混合模型的应用

1、数据分类: 由于观测数据的分布可以由多个单高斯分布进行近似。通过判断观测数据属于哪个单高斯分布,从而对数据进行分类。

2、异常检测: 数据的异常检测在金融风控、设备状态监测等场景应用广泛。GMM首先对观测数据进行多高斯分布建模,并通过EM算法估计单高斯分布参数,如果测试数据(在线数据)在隶属的高斯分布中的概率密度值小于阈值,则判定为异常点。

3、图像分割:在医学图像分割中应用广泛,通常用于前景以及不同类别的背景之间的分割,本质上与数据分类相同。

4、图像生成: 近期在各大平台上火热的AI绘画,为生成式模型作品。GMM能够用于图像生成关键在于其能够通过不同的高斯分布的线性组合对任意图像分布进行拟合。因此,我们通过GMM求出了图像分布,就可以从图像分布中任意采样生成与原数据相似但不同的图像。

2.1 异常检测

  • 我们从正态分布中生成了800个点,从均匀分布中生成了20个点作为异常点。然后为了对比模型效果我们记录了数据是不是异常值的标签(真实数据集一般没有)。

  • 准确率99.87%,只有1个点预测错误。而且这个点,和正常值太接近了,预测错了完全可以理解。

  • 高斯混合模型是基于分布的,而我们则模拟数据集就是从不同分布里面生成的,所以正好对上了它的专长,所以效果会比较好。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.mixture import GaussianMixture
from sklearn.metrics import accuracy_score


def get_data():
    # 生成数据:正常数据和异常数据
    np.random.seed(77)
    X_normal = 0.3 * np.random.randn(800, 2)                     # 正常数据
    X_abnormal = np.random.uniform(low=-4, high=4, size=(20, 2)) # 异常数据
    X = np.vstack([X_normal, X_abnormal])

    # 标准化数据
    scaler = StandardScaler()
    scaled = scaler.fit(X_normal)
    X_scaled = scaled.transform(X)
    # 正常点label为1,异常点label为-1
    true_labels = np.hstack([np.ones(len(X_normal)), -1 * np.ones(len(X_abnormal))])
    return X_normal, X_abnormal, X, X_scaled, true_labels



def plot_result(y_pred, model_name=''):
    # Splitting the data into predicted normal and abnormal by gmm
    X_k_normal = X[y_pred == 1]
    X_k_abnormal = X[y_pred == -1]

    # Correctly and incorrectly predicted points by gmm
    correctly_predicted_k = y_pred == true_labels
    incorrectly_predicted_k = ~correctly_predicted_k
    accuracy_k_simple = accuracy_score(true_labels, y_pred)
    print(accuracy_k_simple)

    plt.figure(figsize=(15, 5))
    # Original data plot
    plt.subplot(1, 3, 1)
    plt.scatter(X_normal[:, 0], X_normal[:, 1], color='blue', label='Normal', s=20, marker='x')
    plt.scatter(X_abnormal[:, 0], X_abnormal[:, 1], color='red', label='Abnormal', s=20, marker='x')
    plt.title("Original Data")
    plt.legend()

    # Predicted data plot
    plt.subplot(1, 3, 2)
    plt.scatter(X_k_normal[:, 0], X_k_normal[:, 1], color='b', label='Predicted Normal', s=20, marker='x')
    plt.scatter(X_k_abnormal[:, 0], X_k_abnormal[:, 1], color='r', label='Predicted Abnormal', s=20, marker='x')
    plt.title(f"Predicted Data ({model_name})")
    plt.legend()

    # Prediction accuracy plot
    plt.subplot(1, 3, 3)
    plt.scatter(X[correctly_predicted_k, 0], X[correctly_predicted_k, 1], color='gold', label='Correctly Predicted',
                s=20, marker='x')
    plt.scatter(X[incorrectly_predicted_k, 0], X[incorrectly_predicted_k, 1], color='purple',
                label='Incorrectly Predicted', s=20, marker='x')
    plt.title(f"Prediction Accuracy ({model_name})")
    plt.legend()

    plt.show()


if __name__ == '__main__':
    X_normal, X_abnormal, X, X_scaled, true_labels = get_data()
    # 构建gmm模型进行异常点检测
    gmm = GaussianMixture(n_components=2, random_state=77)
    gmm_labels = gmm.fit_predict(X_scaled)
    plot_result(np.where(gmm_labels == 0, 1, -1), model_name='GaussianMixture')

在这里插入图片描述

2.2 数据聚类

import numpy as np
import matplotlib.pyplot as plt
from sklearn.mixture import GaussianMixture

# 生成随机数据
np.random.seed(0)
n_samples = 500
X1 = np.random.randn(n_samples, 2)                                # 均值0,方差为1
X2 = np.random.normal(5, 1, 2 * n_samples).reshape(n_samples, 2)  # 均值5,方差为1

X = np.vstack([X1, X2])

# 创建GMM模型对象并拟合数据
n_clusters = 2
gmm = GaussianMixture(n_components=n_clusters).fit(X)

# 预测每个样本所属的聚类
labels = gmm.predict(X)

# 可视化聚类结果
plt.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis')
plt.show()

在这里插入图片描述

2.3 手写数字生成

  • 我们以手写体数字生成为例,利用GMM拟合图像数据分布,并从拟合后的分布中对数据进行采样,从而生成新的手写体数字。
  • 我们先下载数据集,并保存本地,方便加载数据。
from sklearn.datasets import fetch_openml
import pickle

# 获取数据集
# False表示以原始格式返回,每个特征是一个单独的数组。True表示返回Pandas
# DataFrame对象
# 自0.24.0(2020 年 12 月)以来,as_frame参数为auto(而不是之前的False默认选项)
mnist = fetch_openml('mnist_784', version=1, as_frame=False)

# 保存数据集到本地
with open('mnist_data.pkl', 'wb') as f:
    pickle.dump(mnist, f)

首先,加载手写体数字数据集,并对图像进行可视化。从70000张数据中选择2000张作为观测数据。

import pickle
from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA


def get_data():
    # 从本地导入数据集
    with open('mnist_data.pkl', 'rb') as f:
        mnist = pickle.load(f)
        # 一共7w张图,每个有784个特征。每个特征是28 x 28 像素中的一个点的数值,在0(白)~ 255(黑)之间。
        X, y = mnist["data"], mnist["target"]
    return X, y

def show_images(train_images_part, train_labels_part, picture_type):
    # 显示前10张图片
    fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(8, 4))
    for ax, image, label in zip(axes.ravel(), train_images_part[:10], train_labels_part[:10]):
        ax.imshow(image, cmap=plt.cm.gray_r)
        ax.set_xticks([])
        ax.set_yticks([])
        if label:
            ax.set_title('Label: {}'.format(label))
    plt.suptitle(picture_type, fontsize=16)
    plt.tight_layout()
    plt.show()

if __name__ == '__main__':
    train_images, train_labels = get_data()
    # 1、选取2000张图片
    train_images_part = train_images[:2000]
    train_labels_part = train_labels[:2000]

    # 展示原始图片
    show_images(train_images_part.reshape(-1, 28, 28), train_labels_part, 'original picture')

在这里插入图片描述

其次,由于图像为28*28=784维,GMM对高维数据拟合难以收敛,因此,对数据通过PCA进行降维。我们保留85%的主成分,将图像维度降为56维。降维后的图像如下图所示,图像虽部分边缘模糊,但保留了原始图像的主要特征。

    # 2、对图片进行pca降维,将图像维度降为56维
    pca = PCA(0.85, whiten=True)
    train_images_pca = pca.fit_transform(train_images_part)
    # 将经过降维处理的图像转换回原始的高维空间(784维)
    train_images_pca_rec = pca.inverse_transform(train_images_pca)

    show_images(train_images_pca_rec.reshape(-1, 28, 28), train_labels_part, 'pca picture')

在这里插入图片描述

最后,我们通过GMM对(2000,56)数据进行拟合,并从中采样10组数据,重建出生成新图像。如下图所示,我们从GMM拟合数据分布中随机采样了10组点,并通过升维重建出了原始图像。可以发现,新生成的图像与原始图像具有较高相似度。

# 3、构建高斯混合模型,并采样10组数据生成新的图片
    gmm = GaussianMixture(110, covariance_type='full', random_state=42)
    gmm.fit(train_images_pca)
    # 采样10组数据生成新的图片
    # data_new[0]为生成的图像
    data_new = gmm.sample(10)
    # 将56维图像转换回原始的高维空间(784维)
    digits_new = pca.inverse_transform(data_new[0])

    show_images(digits_new.reshape(-1, 28, 28), [None] * 10, 'generate picture')

在这里插入图片描述

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/521164.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

cmake学习笔记1

基础概念 CMake是什么? CMake是一个元构建系统(meta build-system),用于生产其他构建系统文件(如Makefile或Ninja)。 基础操作方式 CMake使用一个CMakeLists.txt文件描述配置,然后使用cmake驱动这个文件生成对应构建系统文件。…

【数据结构】ArrayList详解

目录 前言 1. 线性表 2. 顺序表 3. ArrayList的介绍和使用 3.1 语法格式 3.2 添加元素 3.3 删除元素 3.4 截取部分arrayList 3.5 其他方法 4. ArrayList的遍历 5.ArrayList的扩容机制 6. ArrayList的优缺点 结语 前言 在集合框架中,ArrayList就是一个…

【Linux】环境基础开发工具使用——vim使用

Linux 软件包管理器 yum 什么是软件包 1.在 Linux 下安装软件 , 一个通常的办法是下载到程序的源代码 , 并进行编译 , 得到可执行程序 . 2.但是这样太麻烦了 , 于是有些人把一些常用的软件提前编译好 , 做成软件包 ( 可以理解成 windows 上的安装程序) 放在一个服务器…

C#,简单,精巧,实用的文件夹时间整理工具FolderTime

点击下载本文软件(5积分): https://download.csdn.net/download/beijinghorn/89071073https://download.csdn.net/download/beijinghorn/89071073 百度网盘(不需积分): https://pan.baidu.com/s/1FwCsSz…

ThreadLocal上传下载文件

文章目录 ThreadLocal1.基本介绍1.什么是ThreadLocal?2.示意图 2.快速入门1.创建普通java项目2.编写代码1.T1.java2.T1Service.java3.T2Dao.java4.Dog.java 3.结果 3.ThreadLocal源码解读1.set方法2.set方法总结3.get方法 上传下载文件1.基本介绍1.基本说明2.文件上…

Spring Cloud介绍

一、SpringCloud总体概述 Cloud Foundry Service Broker:通用service集成进入Cloud Foundry Cluster:服务集群 Consul:注册中心 Security:安全认证 Stream:消息队列 Stream App Starters:Spring Cloud Stre…

Redis 客户端

Redis 客户端 客户端-服务器结构 Redis 同 Mysql 一样,也是一个客户端-服务器结构的程序,结构如下图: 注:Redis 客户端和服务器可以在同一个主机上,也可以在不同主机上 Redis 客户端的多种形态 自带的命令行客户端&…

【Qt 学习笔记】详解Qt中的信号和槽

博客主页:Duck Bro 博客主页系列专栏:Qt 专栏关注博主,后期持续更新系列文章如果有错误感谢请大家批评指出,及时修改感谢大家点赞👍收藏⭐评论✍ 详解Qt中的信号与槽 文章编号:Qt 学习笔记 / 12 文章目录…

【Node.js】短链接

原文链接:Nodejs 第六十二章(短链接) - 掘金 (juejin.cn) 短链接是一种缩短长网址的方法,将原始的长网址转换为更短的形式。短链接的主要用途之一是在社交媒体平台进行链接分享。由于这些平台对字符数量有限制,长网址可…

旋转花键有哪些优缺点?

旋转花键是在花键外筒的外径上装上专用的轴承外套,使之运转动作,适用于水平多关节机械手臂(SCARA)、产业用机器人、自动装载机、镭射加工机、搬送装置、机械加工中心的ATC装置等各项设备。 目前,旋转花键的应用越来越普…

redis 哨兵

文章目录 前言主从复制的问题怎么人工恢复故障主节点 Redis Setinel 架构使用 docker 来配置哨兵结构安装 docker编排 redis 主从节点编排 redis 哨兵节点 观察哨兵模式的作用主从切换的具体流程小结 前言 redis 主从复制模式下, 一旦主节点出现故障, 不能提供服务的时候, 就需…

刷题之Leetcode283题(超级详细)

283.移动零 283. 移动零https://leetcode.cn/problems/move-zeroes/ 给定一个数组 nums,编写一个函数将所有 0 移动到数组的末尾,同时保持非零元素的相对顺序。 请注意 ,必须在不复制数组的情况下原地对数组进行操作。 示例 1: 输入: nu…

第十讲 Query Execution Part 1

1 处理模型【Processing Model】 DBMS 的处理模型【Processing Model】定义了系统如何执行【execute】查询计划【Query Plan】。 针对不同的工作负载进行不同的权衡。 方法1:迭代器模型【Iterator Model】 方法2:物化模型【Materialization Model】 方…

创建和启动线程

概述 Java语言的JVM允许程序运行多个线程,使用java.lang.Thread类代表线程,所有的线程对象都必须是Thread类或其子类的实例。 Thread类的特性 每个线程都是通过某个特定Thread对象的run()方法来完成操作的,因此把run()方法体称为线程执行体。…

数据结构之堆底层实现的循序渐进

题外话 把没写的都补回来! 正题 堆 概念 堆是一棵完全二叉树,因此可以层序的规则采用顺序的方式来高效存储, 大根堆:指根结点比左右孩子都大的堆 小根堆:指根结点比左右孩子都小的堆 性质 1.堆中某个节点的值总是不大于或不小于其父节点的值 2…

CCIE-14-MPLS_and_BGP

目录 实验条件网络拓朴 环境配置开始配置配置MPLSR1访问R6检测结果R6访问R1检测结果 实验条件 网络拓朴 环境配置 在我的资源里可以下载&#xff08;就在这篇文章的开头也可以下载&#xff09; 开始配置 R1<->R2&#xff1a;EBGP R2<->R5&#xff1a;IBGP&…

蓝桥杯备考3

P8196 [传智杯 #4 决赛] 三元组 题目描述 给定一个长度为 n 的数列 a&#xff0c;对于一个有序整数三元组 (i,j,k)&#xff0c;若其满足 1≤i≤j≤k≤n 并且&#xff0c;则我们称这个三元组是「传智的」。 现在请你计算&#xff0c;有多少有序整数三元组是传智的。 输入格式…

小米手机澎湃OS,不Root查看电池健康

首先&#xff0c;在键盘拨号界面&#xff0c;输入*#*#284#*#*&#xff0c;会调用问题反馈APP来生成当前系统的故障日志&#xff0c;如果提示你需要授权什么就点确认 稍等几分钟&#xff0c;会得到一个压缩包&#xff0c;保存在目录MIUI/debug_log下 这里为了方便&#xff0c;我…

肖恩带你学C语言·文件操作(上)

1. 为什么使用文件 如果没有文件&#xff0c;我们写的程序的数据是存储在电脑的内存中&#xff0c;如果程序退出&#xff0c;内存回收&#xff0c;数据就丢失了&#xff0c;等再次运行程序&#xff0c;是看不到上次程序的数据的&#xff0c;如果要将数据进行持久化的保存&…

打造自然资源“一张图”管理平台,推动生态文明建设新篇章

在信息化时代的浪潮下&#xff0c;自然资源管理正面临着前所未有的挑战与机遇。传统的资源管理模式已经难以满足当前生态环境保护与经济发展的双重需求&#xff0c;我们需要一个全新的平台&#xff0c;一个集信息集成、智能分析、决策支持于一体的自然资源“一张图”管理平台。…