EM算法與混合高斯模型
EM算法與混合高斯模型
close all;
clear;
clc;
%% Sample Generate
N=5000;
a_real =[3/10,5/10,2/10];
mu_real = [7,12;12,7;14,15];
cov_real(:,:,1) = [1,0;0,1];
cov_real(:,:,2) = [3,1;1,3];
cov_real(:,:,3) = [3,1;1,3];
X_1 = mvnrnd(mu_real(1,:),cov_real(:,:,1),N*a_real(1));
X_2 = mvnrnd(mu_real(2,:),cov_real(:,:,2),N*a_real(2));
X_3 = mvnrnd(mu_real(3,:),cov_real(:,:,2),N*a_real(3));
X=[X_1;X_2;X_3];
X = X(randperm(size(X,1)),:);
%% Sample Ploting
x = 0:0.5:20;
y = 0:0.5:20;
[x y]=meshgrid(x,y);
mesh = [x(:),y(:)];
z_real = a_real(1)*mvnpdf(mesh,mu_real(1,:),cov_real(:,:,1))+...
a_real(2)* mvnpdf(mesh,mu_real(2,:),cov_real(:,:,2))+...
a_real(3)* mvnpdf(mesh,mu_real(3,:),cov_real(:,:,3));
subplot(2,3,1);
plot(X_1(:,1),X_1(:,2),‘x‘,X_2(:,1),X_2(:,2),‘o‘,X_3(:,1),X_3(:,2),‘<‘)
subplot(2,3,2);
contour(x,y,reshape(z_real,size(x,2),size(y,2)));
subplot(2,3,3);
surf(x,y,reshape(z_real,size(x,2),size(y,2)));
subplot(2,3,4);
plot(X(:,1),X(:,2),‘x‘);
%%subplot(2,2,3);
%%surf(x,y,reshape(z_real,size(x,2),size(y,2)));
%%shading interp;
%% Parameter Initialization
a = [1/3, 1/3,1/3];
cov(:,:,1) = [1,0;0,1];
cov(:,:,2) = [1,0;0,1];
cov(:,:,3) = [1,0;0,1];
mu_y_init = (max(X(:,1))+min(X(:,1)))/2;
mu_x1_init = max(X(:,2))/4+3*min(X(:,2))/4;
mu_x2_init = 2*max(X(:,2))/4+2*min(X(:,2))/4;
mu_x3_init = 3*max(X(:,2))/4+1*min(X(:,2))/4;
w = zeros(size(X,1),length(a)); %%w(i,j),i是樣本數量,j是聚類數量,w(i,j)的值表示樣本i屬於聚類j的概率,最大值表示i的聚類
mu = [mu_x1_init,mu_y_init;mu_x2_init,mu_y_init;mu_x3_init,mu_y_init];
%% EM Implementation
iter = 40;
for i = 1:iter
%% Expectaion: 根據現有的(最新得到的)高斯分布和先驗概率,計算每一個樣本屬於每一個聚類的概率
for j = 1 : length(a)
w(:,j)=a(j)*mvnpdf(X,mu(j,:),cov(:,:,j));
end
w=w./repmat(sum(w,2),1,size(w,2)); %%為了方便./計算,建立一個和w一樣規模的矩陣,矩陣的列是相同的,每一行是w(:,1)+w(:,2)的值
%% Maximum: 根據現有的(最新得到的)的每一個樣本屬於每一個聚類的概率,計算先驗概率和高斯參數
a = sum(w,1)./size(w,1); %%把w的每一行加起來就能得到每一個聚類的先驗概率
mu = w‘*X; %%分別得到X每一維對於每一個聚類的期望,mu(i,j),i是維數,j是聚類數
mu= mu./repmat((sum(w,1))‘,1,size(mu,2));
for j = 1 : length(a)
vari = repmat(w(:,j),1,size(X,2)).*(X- repmat(mu(j,:),size(X,1),1)); %%得到一個特定聚類X每一維的方差矩陣,乘以w,(相當於選擇出屬於該聚類的X)
cov(:,:,j) = (vari‘*vari)/sum(w(:,j),1);
end
end
%% Estimation
[c estimate] = max(w,[],2);
%% Estimation Plotting
z = a(1)*mvnpdf(mesh,mu(1,:),cov(:,:,1))+...
a(2)* mvnpdf(mesh,mu(2,:),cov(:,:,2))+...
a(3)* mvnpdf(mesh,mu(3,:),cov(:,:,3));
subplot(2,3,5);
contour(x,y,reshape(z,size(x,2),size(y,2)));
one = find(estimate==1);
two = find(estimate == 2);
three = find(estimate == 3);
% Plot Examples
subplot(2,3,6);
plot(X(one, 1), X(one, 2), ‘x‘,X(two, 1), X(two, 2), ‘o‘, X(three, 1), X(three, 2), ‘<‘);
%plot(X(neg, 1), X(neg, 2), ‘o‘);
%% Build-in EM Implementation
%options = statset(‘Display‘,‘final‘);
%gm = gmdistribution.fit(X,2,‘Options‘,options);
%subplot(2,2,4);
%ezcontour(@(x,y)pdf(gm,[x y]),[0 20],[0 20]);
matlab中各種高斯相關函數
matlab, 高斯函數, 高斯分布
最常見的是產生服從一維標準正態分布的隨機數
- n=100;
- x=randn(1,n)
實現服從任意一維高斯分布的隨機數
- u=10;
- sigma=4;
- x=sigma*randn(1,n)+u
產生服從多元高斯分布的隨機變量函數mvnrnd,[multivarite normal random]
- n=100; %產生隨機數的個數
- mu=[1 -1];
- Sigma=[.9,.4;.4,.3];
- r=mvnrnd(mu,Sigma,n);
將產生的隨機數繪制在二維平面
- scatter(r(:,1),r(:,2));
1474457688429.jpg
當然mvnrnd函數還可以產生更高維數的高斯隨機數,具體參見matlab help。
產生多元高斯分布概率密度函數
Y=mvnpdf(X,[MU,Sigma])
其中可省參數MU,Sigma默認值分別為零向量和單位陣,X是的矩陣,N是樣本個數,D是樣本維數。
- mu = [1 -1]; Sigma = [.9 .4; .4 .3];
- [X1,X2] = meshgrid(linspace(-1,3,25)‘, linspace(-3,1,25)‘);
- X = [X1(:) X2(:)];
- p = mvnpdf(X, mu, Sigma);
- surf(X1,X2,reshape(p,25,25));
和下面代碼產生的趨勢相同
- mu = [1 -1];
- Sigma = [.9 .4; .4 .3];
- [X,Y] = meshgrid(linspace(-1,3,25)‘, linspace(-3,1,25)‘);
- for
i=1:25
-
for
j=1:25
- XY=[X(i,j),Y(i,j)];
- Z(i,j)=exp(-0.5*(XY-mu)/Sigma*(XY-mu)‘);
-
end
- end
- surf(X,Y,Z);
1474460161935.jpg
高斯分布函數
Y=mvncdf(X,[Mu],[Sigma]) , cumulative probability of the multivariate norm distribution with mean Mu and covariance Sigma.
具體使用看代碼
- mu = [1 -1]; Sigma = [.9 .4; .4 .3];
- [X1,X2] = meshgrid(linspace(-1,3,25)‘, linspace(-3,1,25)‘);
- X = [X1(:) X2(:)];
- p = mvncdf(X, mu, Sigma);
- surf(X1,X2,reshape(p,25,25));
1474460555428.jpg
高斯隸屬函數
gaussmf(X,[Sigma,Mu])
- x = (0:0.1:10)‘;
- y1 = gaussmf(x, [0.5 5]);
- y2 = gaussmf(x, [1 5]);
- y3 = gaussmf(x, [2 5]);
- y4 = gaussmf(x, [3 5]);
- plot(x, [y1 y2 y3 y4]);
EM算法與混合高斯模型