demomix.m

这里包含了聚类的工具箱还有很详细的文档说明

clear all

N = 200;
Ki = 4;
d = 2; % The dimension shouln't ne changed in this demo

x = [];
for k = 1:Ki
    theta = rand(1)*pi/2;
    Mr = [cos(theta) -sin(theta); sin(theta) cos(theta)]*diag([1 (rand(1)+0.1)]);
    x = [x (rand(d,1)*ones(1,N) + 0.12*Mr*randn(d,N))];
end
N = N*Ki;

%load dataem; [d,N] = size(x);
%load dataem2

K = 5;

[xi,yi] = meshgrid(linspace(-0.5,1.5,100),linspace(-0.5,1.5,100));
xax = [xi(:) yi(:)]';

% Initialization
for k=1:K
    mui{k} = rand(d,1);
    vari{k} = 0.3*diag(rand(d,1));
end
lambi = rand(1,K); lambi = lambi/sum(lambi);


z = zeros(size(xi));
p = cell(1,K);
for k = 1:K
    p{k} = z;
end

figure(1)
hold off
plot(x(1,:),x(2,:),'.r')


for iter = 1:500
    
    z(:) = 0;
    for k = 1:K
        p{k}(:) = lambi(k)*gaussian(xax,mui{k},vari{k}+eps)';
        z(:) = z(:) + p{k}(:);
    end
    %z(:) = log(z(:));
    figure(1)
    hold off
    plot(x(1,:),x(2,:),'.r')
    hold on
    if 0
        contz = linspace(min(z(:)),max(z(:)),10);
        contour(xi,yi,z,contz,'b')
        axis([-0.4 1.4 -0.4 1.4])
    else
        for k = 1:K
            lim = max(p{k}(:))/(lambi(k)*K)/4;
            contour(xi,yi,p{k},[lim lim],'b');
            axis([-0.4 1.4 -0.4 1.4])
        end
    end
    
    pause(0.01);
    %pause
    
    
    [lambo,muo,varo] = emgaussmix(x,lambi,mui,vari);
    for k=1:K
        lambi(k) = lambo(k);
        mui{k} = muo{k};
        vari{k} = varo{k};
    end
    
end