高斯混合模型R代码

最新推荐文章于 2024-08-23 15:30:37 发布

水令青争傻傻

最新推荐文章于 2024-08-23 15:30:37 发布

阅读量3.4k

点赞数 3

分类专栏：聚类文章标签：机器学习 r语言

原文链接：http://blog.pluskid.org/?p=39

版权

聚类专栏收录该内容

0 篇文章 1 订阅

订阅专栏

高斯混合模型R代码

关于本文的说明

关于本文的说明

作为一个基础很差的小白，主要参考了几位大神的博客，将matlab代码改成了R代码，算是记录一下我的学习。本文没有获得原作者授权，仅作为学习分享的心得。

原文链接如下

http://blog.pluskid.org/?p=39
https://blog.csdn.net/haoxiaodao/article/details/39480883
https://blog.csdn.net/zhouxinxin0202/article/details/79417371

一些关键的公式

1

$\gamma(i,k) = \frac{\pi_kN(x_i | \mu_k,\Sigma_k)}{\sum_{j=1}^K \pi_kN(x_i | \mu_j,\Sigma_j)}$

上式即每个数据由第k个高斯分布生成概率。注意此时 $x_i$ 理解为每一行的样本。
式子中 $\mu_k$ , $\Sigma_k$ 即每次需要估计和更新的参数。比较容易想到这两个参数的估计分别可以用上次划分后该类样本的平均向量和协方差矩阵得到。

2

$N_k = \sum_{i=1}^N\gamma(i,k) ..............2.1$
上式中 $N_k$ 反映的是什么？
$\pi_k = N_k/N ...................2.2$
上式中 $\pi_k$ 反映的是第k个高斯分布被选中的概率
$\mu_k = \frac{1}{N_k}{\sum_{i=1}^N\gamma(i,k)x_i} ..........2.3$
上式即对均值向量的估计。
$\Sigma_k = \frac{1}{N_k}{\sum_{i=1}^N\gamma(i,k)(x_i-\mu_k)(x_i-\mu_k)^T} .............2.4$
上式即对协方差矩阵的估计。

$N_j(x|\mu_j,\Sigma_j) = \frac{1}{(2\pi)^K|\Sigma_j|}{exp{[-\frac{1}{2}(x-\mu_j)^T\Sigma_j^{-1}(x-\mu_j)]}} .........2.5$
即第j个高斯分布的概率产生

R代码(分段注解)

(基础太差，所以比较啰嗦，完整代码见最后)
输入相关参数:

# 首先接受原始的矩阵X和类别数目K_or_centroids（或类别的中心向量）
N = dim(X)[1]# 统计样本数量
D = dim(X)[2]# 统计样本的维度

if(length(K_or_centroids)==1){
    K = K_or_centroids
    # 随机选择中心点
    rndp = sample(1:N)# 感觉直接sample出k个应该也是可以的--------
    centroids = X[rndp[1:K],]# 这个中心的选取就是打乱后的前K个样本
  }else{
    K = nrow(K_or_centroids)
    centroids = K_or_centroids
  }

初始化高斯分布的参数

# 初始化高斯分布的参数大小：
  init_pars = function(centroids,K,D){
    pMiu = centroids# 先将之前随机得到的中心即作为均值向量
    pPi = matrix(0,nrow = 1, ncol = K)# 构造pi_k
    pSigma = array(0,dim=c(D,D,K))# 构造sigma
    
    # 按照欧式距离分配原来的样本分别到K个分类
    tst1 = apply(X^2,1,sum)
    tst1 = rep(tst1,K)
    tst1 = matrix(tst1,ncol=3)
    
    tst2 = t(apply(pMiu^2,1,sum))
    tst2 = rep(tst2,N)
    tst2 = matrix(tst2,ncol=3,byrow=TRUE)
    
    tst3 = 2*as.matrix(X)%*%t(pMiu)
    # 关于这里为什么产生的是欧式距离矩阵，推导可以见
    # https://blog.csdn.net/IT_forlearn/article/details/100022244
    distmat = tst1+tst2-tst3
    
    # 寻找每行的最小值---也就是到哪一类最近即将样本分为该类
    labels = apply(distmat,1,which.min)
    for(k in 1:K){
      Xk = X[labels==k,]# 某个中心对应的样本
      pPi[1,k] = nrow(Xk)/N# 利用每类的样本量来估计该分布的贡献大小
      pSigma[,,k] = cov(Xk)# 其实就是通过那些样本的协方差矩阵估计的
    }
    return(list(pMiu,pPi,pSigma))
    # 至此我们通过已有的中心向量能够估计出上述的三个参数,当然pMiu并没有改变
  }

计算概率

calc_prob = function(){
    Px = matrix(0,N,K)# 即N个样本在当前K类条件下的概率
    for(k in 1:K){
      # 首先计算每个样本（行向量）与第k个中心的差值
      Xshift = as.matrix(X) - matrix(rep(as.matrix(pMiu[k,]),N),nrow=N,byrow = TRUE)
      # 上一行感觉写的比较麻烦，其实可以直接减
      
      inv_pSigma = solve(pSigma[,,k])# 求第k类的逆协方差矩阵
      # Xshift即X-μ---对应公式
      tmp = apply((Xshift%*%inv_pSigma)*Xshift,1,sum)
      coef = (2*pi)^(-D/2)*sqrt(det(inv_pSigma))

      # 第k个高斯混合模型得出的Px似然吗？
      Px[,k] = coef*exp(-0.5*tmp)# 公式2.5
    }
    # 最终得到了在当前高斯参数情况下的各样本的概率似然
   return(Px) 
  }

主要部分，实现E-M算法

Lprev = -Inf # 原文为-lnf

while(TRUE){
    # 根据初始化的参数首先计算出Px
    Px = calc_prob()
    
    # 当前已知k个高斯分布的贡献和各样本概率
    pGamma = Px*matrix(rep(pPi,times=N),nrow=N,byrow = TRUE)
    tst = rep(apply(pGamma,1,sum),times = K)# 这个应该就是分母上所有的分布的概率
    pGamma = pGamma/matrix(tst,nrow = N,byrow = FALSE)# 公式1
    
	#“或者说，我们可以看作 x_i 这个值其中有gamma 这部分是由 Component k 所生成的”
	#“集中考虑所有的数据点，现在实际上可以看作 Component 生成了gamma(N, k)x_N 这些点”
	
    # 每个组分的新参数
    Nk = apply(pGamma,2,sum)# 公式2.1---所有样本分别由K个分布按照pPi和Px生成
    pMiu = diag(1/Nk)%*%t(pGamma)%*%X# 公式2.3
    pPi = Nk/N# 公式2.2
    for (kk in 1:K){
      tst = rep(as.numeric(pMiu[kk,]),times=N)
      # 这里也是只取了前三个，所以其实将变动的范围都固定到了前三？
      Xshift = as.matrix(X)-matrix(tst,nrow=N,byrow = TRUE)
      pSigma[,,k] = (t(Xshift)%*%(diag(pGamma[,kk])%*%Xshift))/Nk[kk]
    }
    
    # 检查收敛情况---由log-likelyhood所产生,gamma的分子部分
    L = sum(log(Px%*%matrix(pPi,nrow=3)))
    if(L-Lprev < threshold){
      break
    }
    Lprev = L

完整代码如下

GMM_copy = function(X,K_or_centroids,nargout=1){
	init_pars = function(centroids,K,D){
   	 	pMiu = centroids
    	pPi = matrix(0,nrow = 1, ncol = K)
    	pSigma = array(0,dim=c(D,D,K))
    
    	tst1 = apply(X^2,1,sum)
    	tst1 = rep(tst1,K)
    	tst1 = matrix(tst1,ncol=3)
    	tst2 = t(apply(pMiu^2,1,sum))
    	tst2 = rep(tst2,N)
    	tst2 = matrix(tst2,ncol=3,byrow=TRUE)
    	tst3 = 2*as.matrix(X)%*%t(pMiu)
    	distmat = tst1+tst2-tst3
    	labels = apply(distmat,1,which.min)

	    for(k in 1:K){
      		Xk = X[labels==k,]
      		pPi[1,k] = nrow(Xk)/N
      		pSigma[,,k] = cov(Xk)
    	}
    	return(list(pMiu,pPi,pSigma))
  	}
  
  	calc_prob = function(){
    	Px = matrix(0,N,K)
    	for(k in 1:K){
      		Xshift = as.matrix(X) - matrix(rep(as.matrix(pMiu[k,]),N),nrow=N,byrow = TRUE)
      		inv_pSigma = solve(pSigma[,,k])# 求逆矩阵
      		tmp = apply((Xshift%*%inv_pSigma)*Xshift,1,sum)
      		coef = (2*pi)^(-D/2)*sqrt(det(inv_pSigma))
      		Px[,k] = coef*exp(-0.5*tmp)
    	}
		return(Px)}




  # nargout表示是否输出所有的变量
  threshold = 10^(-15)
  N = dim(X)[1]
  D = dim(X)[2]
  
  if(length(K_or_centroids)==1){
    K = K_or_centroids
    rndp = sample(1:N)
    centroids = X[rndp[1:K],]
  }else{
    K = nrow(K_or_centroids)
    centroids = K_or_centroids
  }
  
  paras = init_pars(centroids,K,D)
  pMiu = paras[[1]]
  pPi = paras[[2]]
  pSigma = paras[[3]]
  
  Lprev = -Inf

  
  while(TRUE){
    Px = calc_prob()
    
    pGamma = Px*matrix(rep(pPi,times=N),nrow=N,byrow = TRUE)
    tst = rep(apply(pGamma,1,sum),times = K)
    pGamma = pGamma/matrix(tst,nrow = N,byrow = FALSE)
    
    Nk = apply(pGamma,2,sum)
    pMiu = diag(1/Nk)%*%t(pGamma)%*%X
    pPi = Nk/N
    for (kk in 1:K){
      tst = rep(as.numeric(pMiu[kk,]),times=N)
      Xshift = as.matrix(X)-matrix(tst,nrow=N,byrow = TRUE)
      pSigma[,,k] = (t(Xshift)%*%(diag(pGamma[,kk])%*%Xshift))/Nk[kk]
    }
    
    L = sum(log(Px%*%matrix(pPi,nrow=3)))
    if(L-Lprev < threshold){
      break
    }
    Lprev = L
    
    if(nargout==1){
      return(list(Px))
    }else{
      Model = list()
      Model$Miu = pMiu
      Model$Sigma = pSigma
      Model$Pi = pPi
      return(list(Px,model))
    }
  }
  
  
}
# 测试
data(iris)
test_result = apply(GMM_copy(iris[,-ncol(iris)],3),1,which.max)
table(test_result,iris[,4])

后记

如果我们要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这 K 个 Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数πk，选中了 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。

通过这种方式计算欧式距离是第一次知道，相比以前一直两两跑循环更简单；
截至完成，尚未了解Nk可以直接加和得到的原理
关于协方差矩阵奇异，请移步原链接文章