EM algorithm for parameter estimation of Gaussian mixture model

1 #Coding:utf-82 ImportNumPy as NP3 4 defQQ (Y,alpha,mu,sigma,k,gama):#Calculate Q function5gsum=[]6n=Len (y)7      forKinchRange (K):8Gsum.append (Np.sum ([gama[j,k] forJinchrange (n)]))9     returnNp.sum ([G*np.log (AK) forG,akinchZip (Gsum,alpha)]) +TenNp.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)) One                      forJinchRange (n)] forKinchRange (K)])#the formula in "statistical learning method" is 9.29 wrong A  - defPhi (Mu,sigma,y):#Calculate Phi -     return1/(Np.sqrt (2*np.pi*sigma) *np.exp (-(Y-MU) **2/2/Sigma)) the  - defGama (ALPHA,MU,SIGMA,I,K):#Calculate Gama -Sumak=np.sum ([[A*phi (M,s,i)] forA,m,sinchzip (alpha,mu,sigma)]) -     returnAlpha[k]*phi (mu[k],sigma[k],i)/Sumak +  - defDatan (length,k):#Generate Data +Y=[np.random.normal (5*j,j+5,length/k) forJinchrange (k)] A     returny at  - defEM (y,k,iter=1000):#Kmeans Algorithm -n =Len (y) -sigma=[10]*K -mu=Range (K) -Alpha=np.ones (K) inqqold,qqnew=0,0 -      forItinchRange (ITER): toGama2=Np.ones ((n,k)) +          forKinchRange (K): -              forIinchrange (n): thegama2[i,k]=Gama (alpha,mu,sigma,y[i],k) *          forKinchRange (K): $Sum_gama=np.sum ([Gama2[j,k] forJinchrange (n)])Panax NotoginsengMu[k]=np.sum ([Gama2[j,k]*y[j] forJinchRange (n)])/Sum_gama -Sigma[k]=np.sum ([gama2[j,k]* (Y[j]-mu[k]) **2 forJinchRange (n)])/Sum_gama thealpha[k]=sum_gama/N +qqnew=QQ (Y,ALPHA,MU,SIGMA,K,GAMA2) A         ifABS (qqold-qqnew) <0.000001: the              Break +Qqold=qqnew -     returnAlpha,mu,sigma $  $N = 500 -k=2 -Data=Datan (n,k) theY=np.reshape (data, (1, N)) -A,b,c =EM (y[0], k)Wuyi PrintA,b,c the #iter=180 - #[ 0.57217609 0.42782391] [4.1472879054766887, 0.72534713118155769] [44.114682884921415, 24.676116557533351] Wu  -Sigma = 6#data on the Web AboutMIU1 = 40 $MIU2 = 20 -X = Np.zeros ((1, N)) -  forIinchxrange (N): -     ifNp.random.random () > 0.5: AX[0, I] = NP.RANDOM.RANDN () * Sigma +miu1 +     Else: theX[0, I] = NP.RANDOM.RANDN () * Sigma +MIU2 -A,b,c =EM (x[0], k) $ PrintA,b,c the #iter=114 the #[ 0.44935959 0.55064041] [40.561782615819361, 21.444533254494189] [33.374144230703514, 51.459622219329155]

EM algorithm for parameter estimation of Gaussian mixture model