0

I'm trying to implement EM algorithm. I'm using a while loop to iterate while the change in my output matrix is larger than some threshold but I receive an error saying "Error in while (delta >= tau) { : missing value where TRUE/FALSE needed" Below is my code.

Multinomial<-function(H,K,tau) {
H=H+.01
count = 0
a_mat = 0
a_data = list()     #list to store old and new A matrices
delta = Inf     #initialize delta

#initialize mix weights and centroid matrices

c_mat = matrix(1/K,nrow(H),K)
rand_row = H[sample(1:40000,1),]        
t_mat = matrix(rand_row/sum(rand_row),ncol(H),1)
for(i in 1:(K-1)) {
    rand_row = H[sample(1:40000,1),]
    t_mat = cbind(t_mat,rand_row/sum(rand_row))
}   

while(delta>=tau) {
    count=count + 1
    a_data[[1]] = 0         #initialize for new A matrix
    a_data[[2]] = a_mat     #old A matrix

    #E-Step: computes assignment probability matrix
    phi_mat = exp(H%*%log(t_mat))       
    num_mat = (c_mat*phi_mat)         #numerator of a
    num_mat[num_mat==NaN]<-.00000001
    nrsums=rowSums(num_mat)
    a_mat = num_mat/nrsums  #assignment prob matrix
    a_data[[1]] = a_mat

    #M-Step: compute new mix weights and centroids
    acol_sum=.colSums(a_mat,nrow(a_mat),ncol(a_mat),na.rm=T)
    c_k = acol_sum/nrow(H)
    c_mat = matrix(c_k,nrow(H),K)
    b_mat = t(H)%*%a_mat
    t_mat = b_mat/colSums(b_mat)

    #compute measure of change (delta)
    if(count>1){
        delta = norm(a_data[[1]]-a_data[[2]],'O')
    }   
}
David Heredia
  • 1
  • 1
    Have you searched for the error message? The answer is always the same. There is an `NA` in either `delta` or `tau`. – Roland Nov 13 '13 at 08:14
  • Have you tried `backtracek()` or `browser()`? – Roman Luštrik Nov 13 '13 at 10:12
  • Looks like your code always has `a_data[[1]]` and `a_data[[2]]` both equal to `a_mat` when they hit the `delta=norm...` statement. That won't solve your question, which is more likely due to entering an illegal value for `tau` when calling the function. – Carl Witthoft Nov 13 '13 at 12:32
  • It occurs to me that if your matrix contains an NA, then `norm` will return NA. – Carl Witthoft Nov 13 '13 at 12:57

0 Answers0