如何检查R中已实现的期望最大化算法的正确性？

问题描述

下午好！

我问了这个问题很多次，没有任何回复或评论（即使经过交叉验证）！

在R下，我为高斯混合模型实现了期望最大化算法：

k=3   # number of known clusters 
w_k=rep(1,k)/k  # we initialize clusters weights by 1/k for each one 
n_j=rep(0,k)    # will be used later
print(w_k)      # printing
data=as.matrix(iris[1:150,-5]) # numerical datasets with 4-dimensions/axis,fifth-axis contains clusters labels.

means=sample(1:dim(data)[1],k,replace=FALSE) # We choose "k" vectors to be intial clusters means. means contains rows indices
mu=as.matrix(iris[means,-5])  # We retreive those means shuffled at random in a matrix of k rows and 4 columns.
sigma=cov(data)      # We compute covariance matrix for the whole dataset. 
print(sigma)
print(solve(sigma))
sigma_list=rep(list(sigma),k) # For each of the k clusters,we itinialize it's covariance matrix by sigma. 
 
# w_k,mu and sigma_list are the 3 parameters to update during the M-step of Esperance maximization ( EM-algorithm).

# Function to compute P(Xi/Cj).
# example : P_Xi_Cj(x=data[1,],mu = mu[3,sigma= sigma) // P(X1/C3)  liklihood of obser1 given cluster3.
#i.e : liklihood of observation given a cluster
P_Xi_Cj <- function(Xi,Cj,sigma= sigma) {
    x=Xi
    mu=Cj
    d=length(x)
    if (length(as.vector(sigma))==1 ){
        if(length(as.vector(x))==1){
            if(length(as.vector(mu))==1){
             return(as.numeric(1/sqrt(abs(1/sigma)*(2*pi)^d)*exp(-1/2*(x)%*%(1/sigma)%*%(x))))
             break  
            }
        }  }
    
    tr1=matrix(x-mu,ncol=1)
    tr2=matrix(x-mu,nrow=1)
    inverse=solve(sigma)
    total=tr2 %*% inverse %*% tr1
    return(as.numeric(1/sqrt(abs(det(sigma))*(2*pi)^d)*exp(-1/2*total)))
    
  }



# Computing P(Cj/Xi) : what is the more likely cluster given the observation ? 

P_Cj_Xi<-function(Xi,mu=mu,sigma_list=sigma_list,W_k=w_k,n_clusters=k){
k=n_clusters
n_j=rep(0,n_clusters)
P_C_Xi=rep(0,n_clusters)    
r=matrix(NA,length(Xi),length(Xi))    
r=rep(list(r),n_clusters )    
r=lapply(1:n_clusters,function(i) r[[i]]=solve(matrix(unlist(sigma_list[i]),ncol=length(Xi))))     
n_j=sapply(1:n_clusters,function(i) -1/2*(rbind(Xi-mu[i,]))%*%as.matrix(r[[i]])%*%(cbind(Xi-mu[i,])))          
total=sum(((sapply(1:n_clusters,function(i) abs(det(as.matrix(r[[i]])))))^(-1/2))*exp(n_j)*W_k)
P_C_Xi=((sapply(1:n_clusters,function(i) abs(det(as.matrix(r[[i]])))))^(-1/2))*exp(n_j)*W_k/total
names(P_C_Xi)=paste("P(Cj/Xi)",1:n_clusters)
               
return(P_C_Xi)        
}
# example :  P_Cj_Xi(Xi=data[1,n_clusters=k)  i.e P(Cj/X1) j=1,..,n_clusters
                
M_step<-function(data=data,mu,sigma_list,W_k,n_clusters=k){
l=lapply(1:nrow(data),function(i) P_Cj_Xi(Xi=data[i,n_clusters=k) )  # E-step 
#print(l)       
sum_P_Cj_Xi=as.vector(Reduce("+",l)) 
#print("sum_P_Cj_Xi")
#print(sum_P_Cj_Xi)          
W_j=sum_P_Cj_Xi/nrow(data) #updating clusters weights (7)
mu=t(sapply(1:k,function(j) Reduce("+",lapply(1:nrow(data),function(i) l[[i]][j]*data[i,]/sum_P_Cj_Xi[j])))) #updating clusters means 
sigma=lapply(1:k,function(i) l[[i]][j]*(data[i,]-mu[j,])%*%t(data[i,])/sum_P_Cj_Xi[j]))) #updating clusters cov matrices         
print(list(mu,sigma,W_j)) #printing for debuggging
return(list(mu,W_j)) 
                                                
}                

max=6                            
t <-max  # number of total iterations
while (t <= max) {
if(t==0) break 
if(t==max){
mu1=mu
sigma=sigma_list 
w_j=w_k  
tmp=M_step(data=data,mu=mu1,sigma_list=sigma,W_k=w_j,n_clusters=k)
t=t-1      
}else{
print(c("iteration : ",t))         
mu1=tmp[[1]]
sigma=tmp[[2]]
w_j=tmp[[3]]   
tmp=M_step(data=data,n_clusters=k) 
   
if(t==0) break                                   
t = t-1    
}    

                              
}

此实现基于本文第二部分中解释的方程式：https://arxiv.org/ftp/arxiv/papers/1603/1603.07879.pdf

我的问题很简单，根据EM算法：

获得的更新参数（mu，sigma，w_k）的结果/输出是否正确/（与em方法匹配）？

我希望我的问题简单明了！

谢谢您的帮助！

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）