高斯混合模型参数估计的EM算法

 1 # coding:utf-8
 2 import numpy as np
 3 
 4 def qq(y,alpha,mu,sigma,K,gama):#计算Q函数
 5     gsum=[]
 6     n=len(y)
 7     for k in range(K):
 8             gsum.append(np.sum([gama[j,k] for j in range(n)]))
 9     return np.sum([g*np.log(ak) for g,ak in zip(gsum,alpha)])+10            np.sum([[np.sum(gama[j,k]*(np.log(1/np.sqrt(2*np.pi))-np.log(np.sqrt(sigma[k]))-1/2/sigma[k]*(y[j]-mu[k])**2))
11                     for j in range(n)] for k in range(K)])  #《统计学习方法》中公式9.29有误
12 
13 def phi(mu,sigma,y): #计算phi
14     return 1/(np.sqrt(2*np.pi*sigma)*np.exp(-(y-mu)**2/2/sigma))
15 
16 def gama(alpha,mu,sigma,i,k): #计算gama
17     sumak=np.sum([[a*phi(m,s,i)] for a,m,s in zip(alpha,mu,sigma)])
18     return alpha[k]*phi(mu[k],sigma[k],i)/sumak
19 
20 def dataN(length,k):#生成数据
21     y=[np.random.normal(5*j,j+5,length/k) for j in range(k)]
22     return y
23 
24 def EM(y,K,iter=1000): #kmeans算法
25     n = len(y)
26     sigma=[10]*K
27     mu=range(K)
28     alpha=np.ones(K)
29     qqold,qqnew=0,0
30     for it in range(iter):
31         gama2=np.ones((n,K))
32         for k in range(K):
33             for i in range(n):
34                 gama2[i,k]=gama(alpha,mu,sigma,y[i],k)
35         for k in range(K):
36             sum_gama=np.sum([gama2[j,k] for j in range(n)])
37             mu[k]=np.sum([gama2[j,k]*y[j] for j in range(n)])/sum_gama
38             sigma[k]=np.sum([gama2[j,k]*(y[j]-mu[k])**2 for j in range(n)])/sum_gama
39             alpha[k]=sum_gama/n
40         qqnew=qq(y,alpha,mu,sigma,K,gama2)
41         if abs(qqold-qqnew)<0.000001:
42             break
43         qqold=qqnew
44     return alpha,mu,sigma
45 
46 N = 500
47 k=2
48 data=dataN(N,k)
49 y=np.reshape(data,(1,N))
50 a,b,c = EM(y[0], k)
51 print a,b,c
52 # iter=180
53 #[ 0.57217609  0.42782391] [4.1472879054766887, 0.72534713118155769] [44.114682884921415, 24.676116557533351]
54 
55 sigma = 6  #网上的数据
56 miu1 = 40
57 miu2 = 20
58 X = np.zeros((1, N))
59 for i in xrange(N):
60     if np.random.random() > 0.5:
61         X[0, i] = np.random.randn() * sigma + miu1
62     else:
63         X[0, i] = np.random.randn() * sigma + miu2
64 a,b,c = EM(X[0], k)
65 print a,b,c
66 # iter=114
67 #[ 0.44935959  0.55064041] [40.561782615819361, 21.444533254494189] [33.374144230703514, 51.459622219329155]