mixtures4.m

无限混合模型的非监督学习算法的例程和MATLAB代码

                      semi_indic = semi_indic(2:k,:);
                else 
                     if comp==k
                        estmu = estmu(:,1:k-1);
                        estcov = estcov(:,:,1:k-1);
                        estpp = estpp(1:k-1);
                        semi_indic = semi_indic(1:k-1,:);
                     else
                        estmu = [estmu(:,1:comp-1) estmu(:,comp+1:k)];
                        newcov = zeros(dimens,dimens,k-1);
                        for kk=1:comp-1
                            newcov(:,:,kk) = estcov(:,:,kk);
                        end
                        for kk=comp+1:k
                            newcov(:,:,kk-1) = estcov(:,:,kk);
                        end
                        estcov = newcov;
                        estpp = [estpp(1:comp-1) estpp(comp+1:k)];
                        semi_indic = semi_indic([1:comp-1,comp+1:k],:);
                     end
                end
            
                % since we've just killed a component, k must decrease
                k=k-1; 
             end    
             
             if (covoption == 2)|(covoption == 3)
                for kk = 1:k
                   semi_indic(kk,:) = multinorm(y,estmu(:,kk),estcov(:,:,kk)); 
                end
             end
             
             if killed==0
                % if the component was not killed, we update the corresponding
                % indicator variables...
                semi_indic(comp,:) = multinorm(y,estmu(:,comp),estcov(:,:,comp)); 
                % ...and go on to the next component
                comp = comp + 1;  
             end
             % if killed==1, it means the in the position "comp", we now
             % have a component that was not yet visited in this sweep, and
             % so all we have to do is go back to the M setp without 
             % increasing "comp"
             
    end % this is the end of the innermost "while comp <= k" loop 
        % which cycles through the components
        
    % increment the iterations counter            
    countf = countf + 1;
    
    clear indic
    clear semi_indic
    for i=1:k
        semi_indic(i,:) = multinorm(y,estmu(:,i),estcov(:,:,i));
        indic(i,:) = semi_indic(i,:)*estpp(i);
    end

    if k~=1
       % if the number of surviving components is not just one, we compute
       % the loglikelihood from the unnormalized assignment variables
       loglike(countf) = sum(log(realmin+sum(indic)));
    else
       % if it is just one component, it is even simpler
       loglike(countf) = sum(log(realmin+indic));
    end
    
    % compute and store the description length and the current number of components
    dlength = -loglike(countf) + (nparsover2*sum(log(estpp))) + ...
              (nparsover2 + 0.5)*k*log(npoints); 
    dl(countf) = dlength;
    kappas(countf) = k;
    
    % compute the change in loglikelihood to check if we should stop
    deltlike = loglike(countf) - loglike(countf-1);
    if (verb~=0)
        disp(sprintf('deltaloglike/th = %0.7g', abs(deltlike/loglike(countf-1))/th));
    end      
    if (abs(deltlike/loglike(countf-1)) < th)
       % if the relative change in loglikelihood is below the threshold, we stop CEM2
       cont=0;
    end

  end % this end is of the inner loop: "while(cont)"
  
  figure 
  if dimens == 1
   [hh,xx] = hist(y,bins);
   barplot = bar(xx,hh/npoints/(xx(2)-xx(1)));
   set(barplot,'EdgeColor',[1 1 1]);
   set(barplot,'FaceColor',[0.75 0.75 0.75]);
   set(gca,'FontName','Times'); set(gca,'FontSize',14);
   drawnow
   hold on
   miny = min(y);
   maxy = max(y);
   plotgrid = miny:(maxy-miny)/200:maxy;
   mix = zeros(size(plotgrid));
   for comp=1:k
       mix = mix + estpp(comp)*uninorm(plotgrid,estmu(comp),estcov(comp));
       plot(plotgrid,estpp(comp)*uninorm(plotgrid,estmu(comp),estcov(comp)),'k')
   end
   plot(plotgrid,mix,'Color','red','LineWidth',3);
   text(plotgrid(5),0.9*max(mix),sprintf('k=%d',k),...
       'FontName','Times','FontSize',18,'FontWeight','bold');
   drawnow
   hold off
 else
  plot(y(axis1,:),y(axis2,:),'.','Color',[0.5 0.5 0.5]);
  placex = get(gca,'Xlim'); placey = get(gca,'Ylim');
  text(placex(1)+0.2,placey(2)-0.2,sprintf('k=%d',k),...
       'FontName','Times','FontSize',22);
  hold on
  axis equal
  set(gca,'FontName','Times','FontSize',22);
  for comp=1:k
      elipsnorm(estmu([axis1,axis2],comp),estcov([axis1,axis2],[axis1,axis2],comp),2)
  end
  drawnow
  hold off
 end
  

    
  % now check if the latest description length is the best; 
  % if it is, we store its value and the corresponding estimates 
  if dl(countf) < mindl
       bestpp = estpp;
       bestmu = estmu;
       bestcov = estcov;
       bestk = k;
       mindl = dl(countf);
  end
  
  % at this point, we may try smaller mixtures by killing the 
  % component with the smallest mixing probability and then restarting CEM2,
  % as long as k is not yet at kmin
  if k>kmin
     [minp indminp] = min(estpp);      
     % what follows is the book-keeping associated with removing one component
     if indminp==1
        estmu = estmu(:,2:k);
        estcov = estcov(:,:,2:k);
        estpp = estpp(2:k);
     else 
        if indminp==k
           estmu = estmu(:,1:k-1);
           estcov = estcov(:,:,1:k-1);
           estpp = estpp(1:k-1);
        else
           estmu = [estmu(:,1:indminp-1) estmu(:,indminp+1:k)];
           newcov = zeros(dimens,dimens,k-1);
           for kk=1:indminp-1
               newcov(:,:,kk) = estcov(:,:,kk);
           end
           for kk=indminp+1:k
               newcov(:,:,kk-1) = estcov(:,:,kk);
           end
           estcov = newcov;
           estpp = [estpp(1:indminp-1) estpp(indminp+1:k)];
        end
     end
     k=k-1;
     
     % we renormalize the mixing probabilities after killing the component
     estpp = estpp/sum(estpp);
     
     % and register the fact that we have forced one component to zero
     transitions2 = [transitions2 countf];
          
     %increment the iterations counter
     countf = countf+1
     
     % ...and compute the loglikelihhod function and the description length
     clear indic
     clear semi_indic
     for i=1:k
        semi_indic(i,:) = multinorm(y,estmu(:,i),estcov(:,:,i));
        indic(i,:) = semi_indic(i,:)*estpp(i);
     end
     if k~=1
        loglike(countf) = sum(log(realmin+sum(indic)));
     else
        loglike(countf) = sum(log(realmin+indic));
     end
     dl(countf) = -loglike(countf) + (nparsover2*sum(log(estpp))) + ...
              (nparsover2 + 0.5)*k*log(npoints); 

     
     kappas(countf) = k;
 
   else %this else corresponds to "if k > kmin"
        %of course, if k is not larger than kmin, we must stop      
        k_cont = 0;
   end
   
end % this is the end of the outer loop "while(k_cont)" 

lastpp = estpp;
lastmu = estmu;
lastcov = estcov;

% finally, we plot the results
figure 
  if dimens == 1
   [hh,xx] = hist(y,bins);
   barplot = bar(xx,hh/npoints/(xx(2)-xx(1)));
   set(barplot,'EdgeColor',[1 1 1]);
   set(barplot,'FaceColor',[0.75 0.75 0.75]);
   set(gca,'FontName','Times'); set(gca,'FontSize',14);
   drawnow
   hold on
   miny = min(y);
   maxy = max(y);
   plotgrid = miny:(maxy-miny)/200:maxy;
   mix = zeros(size(plotgrid));
   for comp=1:bestk
       mix = mix + bestpp(comp)*uninorm(plotgrid,bestmu(comp),bestcov(comp));
       plot(plotgrid,bestpp(comp)*uninorm(plotgrid,bestmu(comp),bestcov(comp)),'k')
   end
   plot(plotgrid,mix,'Color','red','LineWidth',3);
   text(plotgrid(5),0.9*max(mix),sprintf('k=%d',bestk),...
       'FontName','Times','FontSize',18,'FontWeight','bold');
   drawnow
   hold off
else
plot(y(axis1,:),y(axis2,:),'.','Color',[0.5 0.5 0.5]);
hold on
axis equal
set(gca,'FontName','Times','FontSize',22);
placex = get(gca,'Xlim'); placey = get(gca,'Ylim');
text(placex(1)+1,placey(2)-1,sprintf('k=%d',length(bestpp)),...
     'FontName','Times','FontSize',22);
for comp=1:length(bestpp)
    elipsnorm(bestmu([axis1,axis2],comp),bestcov([axis1,axis2],[axis1,axis2],comp),2)
end
drawnow
hold off
end 
thefig=gcf;

% now, we plot the evolution of the description length
% and signal where components where killed
figure
plot(dl,'LineWidth',2)
title('Description Length')
ax = axis;
for i=1:length(transitions1)
    line([transitions1(i) transitions1(i)],[ax(3) ax(4)])
end
for i=1:length(transitions2)
    line([transitions2(i) transitions2(i)],[ax(3) ax(4)],'LineStyle',':')
end
figure(thefig)