function [eta] = Estep (Xs, p, mix_p) N = size(Xs,1); D = size(Xs,2); K = size(p,1); log_eta = zeros(N,K); for i=1:N log_eta(i,:)=(log(mix_p)+sum(log(p(:,logical(Xs(i,:)))),axis=2)+sum(log(1-p(:,logical(1-Xs(i,:)))),axis=2))'; end eta = exp(log_eta - max(log_eta,[],axis=2)); %model.eta = exp(log_eta); eta = eta./sum(eta,axis=2); endfunction function [p, mix_p] = Mstep (Xs, eta, alpha1,alpha2) K = size(eta,2); N = size(Xs,1); D = size(Xs,2); p = zeros(K,D); Nk = sum(eta)'; mix_p = (Nk+alpha2)/(sum(Nk)+alpha2*K); for k=1:K p(k,:) = (sum(spdiags(eta(:,k),0,N,N)*Xs)+alpha1)/(alpha1*D+Nk(k)); end p = max(p,eps); endfunction function [clusters] = MoBlabels_T (Xs, p, mix_p) eta = Estep_T(Xs,p,mix_p); [mval clusters] =max(eta,[],axis=2); endfunction