Блог

Добрый день!

Я задавал этот вопрос много раз без какого-либо ответа или комментария (даже при перекрестной проверке)!

В 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,],mu=mu ,sigma_list=sigma_list,W_k=w_k,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,],mu,sigma_list,W_k,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(j) Reduce(" ",lapply(1:nrow(data), function(i) l[[i]][j]*(data[i,]-mu[j,])%*%t(data[i,]-mu[j,])/sum_P_Cj_Xi[j]))) #updating clusters cov matrices         
print(list(mu,sigma,W_j)) #printing for debuggging
return(list(mu,sigma,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,mu=mu1 ,sigma_list=sigma,W_k=w_j,n_clusters=k) 
   
if(t==0) break                                   
t = t-1    
}    

                              
}
  

Эта реализация основана на уравнениях, описанных в части 2 этой работы: https://arxiv.org/ftp/arxiv/papers/1603/1603.07879.pdf

Мой вопрос прост, в соответствии с алгоритмом EM :

Являются ли полученные результаты / выходные данные обновленных параметров (mu, sigma, w_k) правильными / (соответствуют подходу em)?

Надеюсь, мой вопрос прост и понятен!

Заранее благодарим вас за помощь!