%%聚类算法测试
%%
%%生成数据集
clear
close all
mu = [2 3];
SIGMA = [1 0; 0 2];
n1 = 100;
r = mvnrnd(mu,SIGMA,n1);
plot(r(:,1),r(:,2),'r+');
hold on;
mu = [5 6];
SIGMA = [ 1 0; 0 2];
n2 = 100;
r2 = mvnrnd(mu,SIGMA,n2);
plot(r2(:,1),r2(:,2),'*')
title('original data')
N = n1 + n2;
data = [r;r2];
y = [ones(n1,1);ones(n2,1)+1];
k = randperm(N);
p1 = data(randi(n1,1),:);
p2 = data(randi(n2,1) + n1,:);
data = data(k,:);
%%
%%mixture of gaussian cluster算法
C = 2;
y = y(k);
R = 100;
alpha1 = 0.5;
alpha2 = 0.5;
nu1 = p1;
nu2 = p2;
sigma1 = [0.1,0;0,0.1];
sigma2 = [0.1,0;0,0.1];
gamma = zeros(N,2);
P = zeros(N,2);
while(R)
for i = 1:N
gamma(i,1) = 1./((2*pi)^1*sqrt(det(sigma1))).*exp(-0.5*(data(i,:)-nu1)*sigma1^-1*(data(i,:)-nu1).');
gamma(i,2) = 1./((2*pi)^1*sqrt(det(sigma2))).*exp(-0.5*(data(i,:)-nu2)*sigma2^-1*(data(i,:)-nu2).');
end
P(:,1) = alpha1*gamma(:,1)./(alpha1*gamma(:,1)+alpha2*gamma(:,2));
P(:,2) = alpha1*gamma(:,2)./(alpha1*gamma(:,1)+alpha2*gamma(:,2));
nu1 = P(:,1).'*data./sum(P(:,1));
nu2 = P(:,2).'*data./sum(P(:,2));
sigma11 = zeros(2,2);
sigma22 = zeros(2,2);
for i = 1:N
sigma11 = sigma11 + P(i,1)*(data(i,:)-nu1).'*(data(i,:)-nu1);
sigma22 = sigma22 + P(i,2)*(data(i,:)-nu2).'*(data(i,:)-nu2);
end
sigma1 = sigma11./sum(P(:,1));
sigma2 = sigma22./sum(P(:,2));
alpha1 = sum(P(:,1))/N;
alpha2 = sum(P(:,2))/N;
R = R - 1;
end
[~,index] = max(P.');
C1 = data(index == 1,:);
C2 = data(index == 2,:);
nu = [nu1;nu2];
figure(2)
plot(C1(:,1),C1(:,2),'r+');
hold on;
plot(C2(:,1),C2(:,2),'*')
plot(nu(:,1),nu(:,2),'g^')
title('mixture of gaussian cluster')
