单变量高斯分布
N ( x ∣ μ , σ 2 ) = 1 ( 2 π σ 2 ) 1 / 2 exp { − 1 2 σ 2 ( x − μ ) 2 } mathcal{N}left(x mid mu, sigma^{2}right)=frac{1}{left(2 pi sigma^{2}right)^{1 / 2}} exp left{-frac{1}{2 sigma^{2}}(x-mu)^{2}right} N(x∣μ,σ2)=(2πσ2)1/21exp{−2σ21(x−μ)2}
多变量的高斯分布，D是维数
N ( x ∣ μ , Σ ) = 1 ( 2 π ) D / 2 1 ∣ Σ ∣ 1 / 2 exp { − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) } mathcal{N}(mathbf{x} mid boldsymbol{mu}, boldsymbol{Sigma})=frac{1}{(2 pi)^{D / 2}} frac{1}{|boldsymbol{Sigma}|^{1 / 2}} exp left{-frac{1}{2}(mathbf{x}-boldsymbol{mu})^{mathrm{T}} boldsymbol{Sigma}^{-1}(mathbf{x}-boldsymbol{mu})right} N(x∣μ,Σ)=(2π)D/21∣Σ∣1/21exp{−21(x−μ)TΣ−1(x−μ)}
对于一个混合高斯模型，可以写成多个单个高斯模型通过不同的权重进行叠加,其中 π k pi_k πk表示权重
p ( x ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) p(x)=sum_{k=1}^{K} pi_{k} mathcal{N}left(x mid mu_{k}, Sigma_{k}right) p(x)=k=1∑KπkN(x∣μk,Σk)
现在引入一个K维二进制变量z,表达形式如下
z k ∈ { 0 , 1 } , Σ k z k = 1 z_kin {0,1},Sigma_{k} z_{k}=1 zk∈{0,1},Σkzk=1
这 p ( z k = 1 ) p(z_k=1) p(zk=1)是高斯分布 N ( x ∣ μ k , Σ k ) mathcal{N}left(x mid mu_{k}, Sigma_{k}right) N(x∣μk,Σk)的先验概率,表达如下
p ( z k = 1 ) = π k p(z_k=1)=pi_k p(zk=1)=πk
先验的直观理解就是,不告诉我一个数据点在什么位置,只知道高斯模型的概率。以上式子有 0 ≤ π k ≤ 1 , ∑ k = 1 K π k = 1 0 leq pi_{k} leq 1, sum_{k=1}^{K} pi_{k}=1 0≤πk≤1,∑k=1Kπk=1的约束，联合起来， p ( z ) p(z) p(z)的表达如下
p ( z ) = ∏ k = 1 K π k z k p(z)=prod_{k=1}^{K} pi_{k}^{z_{k}} p(z)=k=1∏Kπkzk
其中 z = [ z 1 , . . . , z k , . . . , z K ] z=[z_1,...,z_k,...,z_K] z=[z1,...,zk,...,zK]。那么 p ( x ∣ z k = 1 ) p(x|z_k=1) p(x∣zk=1)就坍塌成一个高斯分布了。
p ( x ∣ z k = 1 ) = N ( x ∣ μ k , Σ k ) pleft(x mid z_{k}=1right)=mathcal{N}left(x mid mu_{k}, Sigma_{k}right) p(x∣zk=1)=N(x∣μk,Σk)
所以 p ( x ∣ z ) p(x|z) p(x∣z)表示如下
p ( x ∣ z ) = ∏ k = 1 K N ( x ∣ μ k , Σ k ) z k p(x mid z)=prod_{k=1}^{K} mathcal{N}left(x mid mu_{k}, Sigma_{k}right)^{z_{k}} p(x∣z)=k=1∏KN(x∣μk,Σk)zk
根据概率中的联合密度和边缘概率的关系 p ( x ) = ∑ z p ( x , z ) p(x)=sum_{z} p(x, z) p(x)=∑zp(x,z),可以得到
p ( x ) = ∑ z p ( x , z ) = ∑ z p ( z ) p ( x ∣ z ) = ∑ z ∏ k = 1 K π k z k ∏ k = 1 K N ( x ∣ μ k , Σ k ) z k = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) p(x)=sum_{z} p(x, z)=sum_{z}p(z)p(x|z)=sum_{z}prod_{k=1}^{K} pi_{k}^{z_{k}}prod_{k=1}^{K} mathcal{N}left(x mid mu_{k}, Sigma_{k}right)^{z_{k}}=sum_{k=1}^{K} pi_{k} mathcal{N}left(x mid mu_{k}, Sigma_{k}right) p(x)=z∑p(x,z)=z∑p(z)p(x∣z)=z∑k=1∏Kπkzkk=1∏KN(x∣μk,Σk)zk=k=1∑KπkN(x∣μk,Σk)
根据贝叶斯概率分布公式可以得到
p ( z ∣ x ) = p ( z , x ) p ( x ) = p ( z ) p ( x ∣ z ) ∑ j = 1 K p ( x ∣ z j ) p ( z j ) p(mathbf{z} mid mathbf{x})=frac{p(mathbf{z},mathbf{x})}{p(mathbf{x})}=frac{p(mathbf{z}) p(mathbf{x} mid mathbf{z})}{sum_{j=1}^{K} pleft(mathbf{x} mid z_{j}right) pleft(z_{j}right)} p(z∣x)=p(x)p(z,x)=∑j=1Kp(x∣zj)p(zj)p(z)p(x∣z)
令 γ ( z k ) = p ( z k = 1 ) gamma(z_k)=p(z_k=1) γ(zk)=p(zk=1),那么
γ ( z k ) ≡ p ( z k = 1 ∣ x ) = p ( z k = 1 ) p ( x ∣ z k = 1 ) ∑ j = 1 K p ( z j = 1 ) p ( x ∣ z j = 1 ) = π k N ( x n ∣ μ k , Σ k ) ∑ j = 1 K π j N ( x n ∣ μ j , Σ j ) begin{aligned} gammaleft(z_{k}right) equiv pleft(z_{k}=1 mid mathbf{x}right) &=frac{pleft(z_{k}=1right) pleft(mathbf{x} mid z_{k}=1right)}{sum_{j=1}^{K} pleft(z_{j}=1right) pleft(mathbf{x} mid z_{j}=1right)} \ &=frac{pi_{k} mathcal{N}left(mathbf{x}_{n} mid mu_{k}, Sigma_{k}right)}{sum_{j=1}^{K} pi_{j} mathcal{N}left(mathbf{x}_{n} mid mu_{j}, Sigma_{j}right)} end{aligned} γ(zk)≡p(zk=1∣x)=∑j=1Kp(zj=1)p(x∣zj=1)p(zk=1)p(x∣zk=1)=∑j=1KπjN(xn∣μj,Σj)πkN(xn∣μk,Σk)
物理意义
γ n k = p ( z n k = 1 ∣ x n ) gamma_{nk}=p(z_{nk}=1 | x_n) γnk=p(znk=1∣xn):一个点属于哪个类的可能性
上式当我们知道 { π k , μ k , Σ k } {pi_k,mu_k,Sigma_k} {πk,μk,Σk}和 x x x，我们可以得到 γ ( z k ) gamma(z_k) γ(zk)。但是实际是我只知道数据点 x 1 , . . . , x N {x_1,...,x_N} x1,...,xN和聚类的个数，让我们来估计 { π k , μ k , Σ k } , γ ( z k ) {pi_k,mu_k,Sigma_k},gamma(z_k) {πk,μk,Σk},γ(zk),此时我们需要最大似然估计
p ( x ) = ∑ k = 1 K π k N ( x ∣ μ k , Σ k ) ln p ( X ∣ π , μ , Σ ) = ∑ n = 1 N ln { ∑ k = 1 K π k N ( x n ∣ μ k , Σ k ) } begin{aligned} p(mathbf{x}) &=sum_{k=1}^{K} pi_{k} mathcal{N}left(mathbf{x} mid mu_{k}, Sigma_{k}right) \ ln p(mathbf{X} mid pi, mu, Sigma) &=sum_{n=1}^{N} ln left{sum_{k=1}^{K} pi_{k} mathcal{N}left(mathbf{x}_{mathbf{n}} mid mu_{k}, Sigma_{k}right)right} end{aligned} p(x)lnp(X∣π,μ,Σ)=k=1∑KπkN(x∣μk,Σk)=n=1∑Nln{k=1∑KπkN(xn∣μk,Σk)}
其中 X = [ x 1 , . . . , x N ] X=[x_1,...,x_N] X=[x1,...,xN],因此最大似然表达如下
π , μ , Σ = arg max π , μ , Σ ln p ( X ∣ π , μ , Σ ) = arg max π , μ , Σ ∑ n = 1 N ln { ∑ k = 1 K π k N ( x n ∣ μ k , Σ k ) } pi, mu, boldsymbol{Sigma}=underset{pi, mu, boldsymbol{Sigma}}{arg max } ln p(mathbf{X} mid pi, mu, boldsymbol{Sigma})=underset{pi, mu, boldsymbol{Sigma}}{arg max } sum_{n=1}^{N} ln left{sum_{k=1}^{K} pi_{k} mathcal{N}left(mathbf{x}_{n} mid mu_{k}, boldsymbol{Sigma}_{k}right)right} π,μ,Σ=π,μ,Σargmaxlnp(X∣π,μ,Σ)=π,μ,Σargmaxn=1∑Nln{k=1∑KπkN(xn∣μk,Σk)}
对于求解多个变量的最优问题，大多数是固定其中多个变量，只求解一个，进行迭代求解即可。
对上式进行求导可知，省略求导过程
0 = − ∑ n = 1 N π k N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) Σ K ( x n − μ k ) ⇔ 0 = − ∑ n = 1 N γ ( z n k ) Σ K ( x n − μ k ) 0 = - sumlimits_{n = 1}^N {frac{{{pi _k}{mathcal N}left( {{x_n}mid {mu _k},{Sigma _k}} right)}}{{sumnolimits_j {{pi _j}{mathcal N}left( {{x_n}mid {mu _j},{Sigma _j}} right)} }}} {Sigma _K}({x_n} - {mu _k})\ Leftrightarrow 0 = - sumlimits_{n = 1}^N {gamma ({z_{nk}})} {Sigma _K}({x_n} - {mu _k}) 0=−n=1∑N∑jπjN(xn∣μj,Σj)πkN(xn∣μk,Σk)ΣK(xn−μk)⇔0=−n=1∑Nγ(znk)ΣK(xn−μk)
化简得
μ k = ∑ n = 1 N γ ( z n k ) x n N k , N k = ∑ n = 1 N γ ( z n k ) {mu _k} = frac{{sumlimits_{n = 1}^N {gamma ({z_{nk}})} {x_n}}}{{{N_k}}},{N_k} = sumlimits_{n = 1}^N {gamma ({z_{nk}})} μk=Nkn=1∑Nγ(znk)xn,Nk=n=1∑Nγ(znk)
物理意义
N k N_k Nk:分配给聚类K的有效点数，即有多个点属于k这个类
μ k mu_k μk:所有点在k这个类上的加权平均，这个不同于KMeans简单的进行求平均，在这里引入了权重的概念
求导过程省略
Σ k = 1 N k ∑ n = 1 N γ ( z n k ) ( x n − μ k ) ( x n − μ k ) T ， N k = ∑ n = 1 N γ ( z n k ) {Sigma _k} = frac{1}{{{N_k}}}sumlimits_{n = 1}^N {gamma ({z_{nk}})({x_n} - {mu _k})} {({x_n} - {mu _k})^T}，{N_k} = sumlimits_{n = 1}^N {gamma ({z_{nk}})} Σk=Nk1n=1∑Nγ(znk)(xn−μk)(xn−μk)T，Nk=n=1∑Nγ(znk)
物理意义
Σ k Sigma_k Σk:是以 μ k mu_k μk为中心的所有点方差的加权平均值,权重为后验概率 γ ( z n k ) gamma(z_{nk}) γ(znk)
解这个问题时，存在一个限制条件$sumnolimits_{k = 1}^K {{pi _k} = 1} $。关于有限制条件的求导，使用拉格朗日求导法
ln p ( X ∣ π , μ , Σ ) + λ ( ∑ k = 1 K π k − 1 ) ln p(mathbf{X} mid pi, mu, Sigma) +lambda(sum limits_{k=1}^K pi_k-1) lnp(X∣π,μ,Σ)+λ(k=1∑Kπk−1)
对上式进行关于 μ k mu_k μk的求导
∑ n = 1 N N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) + λ = 0 sumlimits_{n = 1}^N {frac{{{mathcal N}left( {{x_n}mid {mu _k},{Sigma _k}} right)}}{{sumnolimits_j {{pi _j}{mathcal N}left( {{x_n}mid {mu _j},{Sigma _j}} right)} }}}+lambda=0 n=1∑N∑jπjN(xn∣μj,Σj)N(xn∣μk,Σk)+λ=0
首先求解 λ lambda λ,等式两边同时乘以 ∑ k = 1 K π k sum limits_{k=1}^K pi_k k=1∑Kπk,化简得
∑ k = 1 K π k ∑ n = 1 N N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) + ∑ k = 1 K π k λ = 0 sum limits_{k=1}^K pi_k {sumlimits_{n = 1}^N {frac{{{mathcal N}left( {{x_n}mid {mu _k},{Sigma _k}} right)}}{{sumnolimits_j {{pi _j}{mathcal N}left( {{x_n}mid {mu _j},{Sigma _j}} right)} }}} }+sum limits_{k=1}^K pi_klambda=0 k=1∑Kπkn=1∑N∑jπjN(xn∣μj,Σj)N(xn∣μk,Σk)+k=1∑Kπkλ=0
∑ k = 1 K ∑ n = 1 N π k N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) + ∑ k = 1 K π k λ = 0 sum limits_{k=1}^K {frac{{sumlimits_{n = 1}^N {pi_k}{mathcal N}left( {{x_n}mid {mu _k},{Sigma _k}} right)}}{{sumnolimits_j {{pi _j}{mathcal N}left( {{x_n}mid {mu _j},{Sigma _j}} right)} }}}+sum limits_{k=1}^K pi_klambda =0 k=1∑K∑jπjN(xn∣μj,Σj)n=1∑NπkN(xn∣μk,Σk)+k=1∑Kπkλ=0
( ∑ n = 1 N 1 ) + λ = 0 λ = − N (sum limits_{n=1}^N 1)+lambda=0\ lambda=-N (n=1∑N1)+λ=0λ=−N
再次将 λ lambda λ代入，进而求解 π k pi_k πk
∑ n = 1 N N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) + λ = 0 1 π k ∑ n = 1 N π k N ( x n ∣ μ k , Σ k ) ∑ j π j N ( x n ∣ μ j , Σ j ) − N = 0 1 π k ∑ n = 1 N γ n k − N = 0 π k = N k N sumlimits_{n = 1}^N {frac{{{mathcal N}left( {{x_n}mid {mu _k},{Sigma _k}} right)}}{{sumnolimits_j {{pi _j}{mathcal N}left( {{x_n}mid {mu _j},{Sigma _j}} right)} }}}+lambda=0 \ frac{1}{pi_k} sumlimits_{n = 1}^N {frac{{{pi_k}{mathcal N}left( {{x_n}mid {mu _k},{Sigma _k}} right)}}{{sumnolimits_j {{pi _j}{mathcal N}left( {{x_n}mid {mu _j},{Sigma _j}} right)} }}}-N=0\ frac {1}{pi_k} sum limits_{n=1}^Ngamma_{nk}-N=0\ pi_k=frac {N_k}{N} n=1∑N∑jπjN(xn∣μj,Σj)N(xn∣μk,Σk)+λ=0πk1n=1∑N∑jπjN(xn∣μj,Σj)πkN(xn∣μk,Σk)−N=0πk1n=1∑Nγnk−N=0πk=NNk
物理意义
π k pi_k πk:表示的是实际上有多少个点属于我这个类$div $所有点的数量
初始化每一个高斯模型的权重 π k pi_k πk,均值 μ k mu_k μk和方差 Σ k Sigma_k Σk
# π(alpha),初始化权重值，即初始时，平分权重
alpha = np.ones(self.n_clusters) / self.n_clusters# μ(mu),随机初始化期望值
mu = np.array([[.403, .237], [.714, .346], [.532, .472]])# ∑(sigma)，这是一个k*data.shape[1]*data.shape[1]三维矩阵，目的是一次性的将方差储存起来,其中初值为对角矩阵，
# 这里以3*2*2，∑为[[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]]],第一块(sigma[0])为k=0对应的方差,同理第二块，第三块也是这样
sigma = np.full((self.n_clusters, n_features, n_features), np.diag(np.full(n_features, 0.1)))
其中 self.n_clusters表示为聚类的个数，n_features表示数据的列数，这里是2，数据是以N*2表示的
E步:通过权重，均值，方差估计后验概率 γ n k = p ( z n k = 1 ∣ x n ) gamma_{nk}=p(z_{nk}=1 | x_n) γnk=p(znk=1∣xn),表示一个点属于哪个类的可能性
gamma = selfputeGamma(data, mu, sigma, alpha, self.multiGaussian) #维数N*k
计算 γ ( z n k ) gamma(z_{nk}) γ(znk)，调用computeGamma函数和multiGaussian函数，函数如下
def computeGamma(self, X, mu, sigma, alpha, multiGaussian):n_samples = X.shape[0]n_clusters = len(alpha)gamma = np.zeros((n_samples, n_clusters))p = np.zeros(n_clusters)g = np.zeros(n_clusters)for i in range(n_samples):for j in range(n_clusters):p[j] = multiGaussian(X[i], mu[j], sigma[j]) #对应公式中的 N(xn|μj，∑j)g[j] = alpha[j] * p[j] #对应公式中的πj*N(xn|μj，∑j)for k in range(n_clusters):if np.sum(g)!=0:gamma[i, k] = g[k] / np.sum(g) #对应公式的γnk的求解return gamma
其中 p [ j ] p[j] p[j]对应 N ( x ∣ μ k , Σ k ) mathcal{N}left(x mid mu_{k}, Sigma_{k}right) N(x∣μk,Σk), g [ j ] g[j] g[j]对应 π k N ( x ∣ μ k , Σ k ) pi_{k} mathcal{N}left(x mid mu_{k}, Sigma_{k}right) πkN(x∣μk,Σk)
def multiGaussian(self, x, mu, sigma):return 1 / ((2 * np.pi) * pow(np.linalg.det(sigma), 0.5)) * np.exp(-0.5 * (x - mu).dot(np.linalg.pinv(sigma)).dot((x - mu).T))
对应于前言的多元高斯分布函数
M步:计算最新的权重，均值，和方差
{ μ k = ∑ n = 1 N γ ( z n k ) x n N k Σ k = 1 N k ∑ n = 1 N γ ( z n k ) ( x n − μ k ) ( x n − μ k ) T π k = N k N N k = ∑ n = 1 N γ ( z n k ) N = ∑ n = 1 N 1 left{ begin{array}{l} {mu _k} = frac{{sumlimits_{n = 1}^N {gamma ({z_{nk}})} {x_n}}}{{{N_k}}}\ {Sigma _k} = frac{1}{{{N_k}}}sumlimits_{n = 1}^N {gamma ({z_{nk}})({x_n} - {mu _k})} {({x_n} - {mu _k})^T}\ {pi _k} = frac{{{N_k}}}{N}\ {N_k} = sumlimits_{n = 1}^N {gamma ({z_{nk}})} \ N = sumlimits_{n = 1}^N 1 end{array} right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧μk=Nkn=1∑Nγ(znk)xnΣk=Nk1n=1∑Nγ(znk)(xn−μk)(xn−μk)Tπk=NNkNk=n=1∑Nγ(znk)N=n=1∑N1
计算权重
alpha = np.sum(gamma, axis=0) / n_samples
对应公式 π k = N k N {pi _k} = frac{{{N_k}}}{N} πk=NNk
计算均值
#更新均值 μ(mu)
gamma_xn=data * gamma[:, i].reshape((n_samples, 1)) #维数N*k
mu_molecule=np.sum(gamma_xn, axis=0) #维数（k,）
mu_denominator=np.sum(gamma, axis=0)[i]
mu[i]=mu_molecule/mu_denominator
其中gamma_xn表示 γ ( z n k ) x n gamma ({z_{nk}}){x_n} γ(znk)xn,mu_molecule表示 ∑ n = 1 N γ ( z n k ) x n sum limits_{n=1}^N gamma(z_{nk})x_n n=1∑Nγ(znk)xn,mu_denominator表示 N k N_k Nk
计算协方差矩阵
sigma[i] = 0
for j in range(n_samples):#更新方差∑(sigma)xn_muk=data[j].reshape((1, n_features))-mu[i] #维数 1*ksigma[i]+=(xn_muk).T.dot(xn_muk)*gamma[j,i] #维数 k*k
sigma[i] = sigma[i] / (np.sum(gamma, axis=0)[i])
其中xn_muk表示为 x n − μ k x_n-mu_k xn−μk。注意这里和公式不一样，是 ( x n − μ k ) T ( x n − μ k ) (x_n-mu_k)^T(x_n-mu_k) (xn−μk)T(xn−μk)这样才可以表示为协方差矩阵，维数为2*2。
评估log函数的，如果收敛，则停止，否则返回到E步
效果如下，不同于KMeans的是，数据点有个渐变的过程，并不是非0即1
import numpy as np
from numpy import *
import pylab
import random, mathimport matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import multivariate_normalplt.style.use('seaborn')class GMM(object):def __init__(self, n_clusters, max_iter=50):self.n_clusters = n_clusters #聚类个数self.ITER = max_iter #最大迭代次数self.mu = 0 #初始化均值self.sigma = 0 #初始化协方差矩阵self.alpha = 0 #初始化权重#计算多元高斯分布函数，注意这里是(x-mu)*sigma^(-1)*(x-mu)^T公式没有错，因为这里算出来的是一个数，并且不能直接运用在三维中，缺了个系数def multiGaussian(self, x, mu, sigma):return 1 / ((2 * np.pi) * pow(np.linalg.det(sigma), 0.5)) * np.exp(-0.5 * (x - mu).dot(np.linalg.pinv(sigma)).dot((x - mu).T))#计算后验概率γnk,维数为data.shape[0]*k,这里以k=2为例，第一列放的是所有属于k=0类的概率，第二列放的是所有属于k=1类的概率def computeGamma(self, X, mu, sigma, alpha, multiGaussian):n_samples = X.shape[0]n_clusters = len(alpha)gamma = np.zeros((n_samples, n_clusters))p = np.zeros(n_clusters)g = np.zeros(n_clusters)for i in range(n_samples):for j in range(n_clusters):p[j] = multiGaussian(X[i], mu[j], sigma[j]) #对应公式中的 N(xn|μj，∑j)g[j] = alpha[j] * p[j] #对应公式中的πj*N(xn|μj，∑j)for k in range(n_clusters):if np.sum(g)!=0:gamma[i, k] = g[k] / np.sum(g) #对应公式的γnk的求解return gamma#GMM的核心代码def fit(self, data):n_samples = data.shape[0]n_features = data.shape[1]'''mu=data[np.random.choice(range(n_samples),self.n_clusters)]'''#初始化初值# π(alpha),初始化权重值，即初始时，平分权重alpha = np.ones(self.n_clusters) / self.n_clusters# μ(mu),随机初始化期望值mu = np.array([[.403, .237], [.714, .346], [.532, .472]])# ∑(sigma)，这是一个k*data.shape[1]*data.shape[1]三维矩阵，目的是一次性的将方差储存起来,其中初值为对角矩阵，# 这里以3*2*2，∑为[[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]],[[0.1,0],[0,0.1]]],第一块(sigma[0])为k=0对应的方差,同理第二块，第三块也是这样sigma = np.full((self.n_clusters, n_features, n_features), np.diag(np.full(n_features, 0.1)))for i in range(self.ITER):#E步，更新后验概率γnk，维数N*kgamma = selfputeGamma(data, mu, sigma, alpha, self.multiGaussian) #M步，更新三大参数：权重，均值，方差#计算权重π(alpha)，维数1*kalpha = np.sum(gamma, axis=0) / n_samplesfor i in range(self.n_clusters):#计算均值 μ(mu)gamma_xn=data * gamma[:, i].reshape((n_samples, 1)) #维数N*kmu_molecule=np.sum(gamma_xn, axis=0) #维数（k,）mu_denominator=np.sum(gamma, axis=0)[i]mu[i]=mu_molecule/mu_denominatorsigma[i] = 0for j in range(n_samples):#计算协方差矩阵∑(sigma)xn_muk=data[j].reshape((1, n_features)) #维数 1*k#维数 k*ksigma[i]+=(xn_muk-mu[i]).T.dot(xn_muk-mu[i])*gamma[j,i] sigma[i] = sigma[i] / (np.sum(gamma, axis=0)[i])self.mu = mu #更新权重π(alpha)self.sigma = sigma #更新均值μ(mu)self.alpha = alpha #更新协方差矩阵∑(sigma)#GMM算法中的最终预测点的分类函数def predict(self, data):pred = selfputeGamma(data, self.mu, self.sigma, self.alpha, self.multiGaussian)cluster_results = np.argmax(pred, axis=1)return cluster_results##初始化画布为4行2列
fig, ax = plt.subplots(4, 2)
## 可视化函数
## 输入：x 原始数的路径
# i 画布的行
# j 画布的列
# n_clusters 聚类的个数，认为指定
def show_result(x,i,j,n_clusters):ax[i][j].scatter(x[:, 0], x[:, 1], s=20, c="b", marker='o')ax[i][j].set_xlabel('x')ax[i][j].set_ylabel('y')gmm = GMM(n_clusters)gmm.fit(x)cat = gmm.predict(x)list_max = max(cat)if list_max == 1:for idx, value in enumerate(cat):if value == 0:ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="r", marker='*')elif value == 1:ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="g", marker='+')elif list_max == 2:for idx, value in enumerate(cat):if value == 0:ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="r", marker='*')elif value == 1:ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="g", marker='+')elif value == 2:ax[i][j+1].scatter(x[idx, 0], x[idx, 1], s=20, c="y", marker='^')ax[i][j+1].set_xlabel('x')ax[i][j+1].set_ylabel('y')## 读取数据的函数
def read_txt(path):filename = path # txt文件和当前脚本在同一目录下，所以不用写具体路径pos = []Efield = []with open(filename, 'r') as file_to_read:while True:lines = file_adline() # 整行读取数据if not lines:breakpassp_tmp, E_tmp = [float(i) for i in lines.split(",")] # 将整行数据分割处理，如果分割符是空格，括号里就不用传入参数，如果是逗号， 则传入‘，'字符。pos.append(p_tmp) # 添加新读取的数据Efield.append(E_tmp)passpos = np.array(pos) # 将数据从list类型转换为array类型。Efield = np.array(Efield)x=np.vstack(pos)y=np.vstack(Efield)X&#atenate((x,y),axis=1)return Xif __name__ == '__main__':x_circle = read_txt(&#")show_result(x_circle, 0, 0, 2)x_moons = read_txt(&#")show_result(x_moons, 1, 0, 2)x_blobs = read_txt(&#")show_result(x_blobs, 2, 0, 3)x_aniso = read_txt(&#")show_result(x_aniso, 3, 0, 3)plt.show()
仿真结果如下
数据集下载参见之前的KMeans文章
最大似然公式存在奇异性的问题：假设GMM的协方差矩阵为 Σ k = σ k 2 I Sigma_k=sigma^2_{k}I Σk=σk2I, I I I为单位向量，有一个数据点 x n = μ j x_n=mu_j xn=μj,那么高斯分布为
N ( x ∣ x n , σ j 2 I ) = 1 ( 2 π ) 1 / 2 1 σ j mathcal{N}left(x mid x_n, sigma_{j}^2Iright)=frac{1}{(2pi)^{1/2}} frac{1}{sigma_j} N(x∣xn,σj2I)=(2π)1/21σj1
那么如果 σ j sigma_j σj非常小，高斯分布函数就非常大，系统可能会崩溃，例如如果一个点就是一个高斯分布，那么 σ = 0 sigma=0 σ=0,N->0, ∏ p ( x i ) prod p(x_i) ∏p(xi)->0。因此，在GMM最大似然公式中，我们将高斯点放在单个点上，这会导致很大的问题。
解决方法：
1.采用MAP的方法，初始考虑 μ k , Σ k mu_k,Sigma_k μk,Σk，假如高斯分布的方差不能太小，就不会掉到奇点问题，但比较复杂
2.如果出现某个高斯点的方差很小，随机初始化另外的一个值。
其次GMM存在和KMeans同样的通病，需要人工给定聚类数。
本文发布于:2024-01-30 13:25:54，感谢您对本站的认可！
本文链接：https://www.4u4v.net/it/170659235620325.html
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系，我们将在24小时内删除。
| 留言与评论（共有 0 条评论） |
