參考資料：機器學習課程的ppt……

Mixture Models

我們將研究混合模型，包括高斯混合模型和伯努利混合模型。

關鍵思想是引入潛變數，它允許從更簡單的分佈形成複雜的分佈。·

我們將看到，混合模型可以用具有離散的潛在變數（在有向的圖形模型中）來解釋。

在後面的課堂上，我們還會看到連續的潛在變數。

K-Means Clustering

k-群集分析

首先，我們來看看下面的問題：在多維空間中識別資料點的簇或組。

我們希望把資料劃分成K簇，其中給出k。

我們觀察到由N維觀測組成的資料集。

其次，我們介紹了D維向量，原型我們可以認為K代表聚類中心。

我們的目標是：找到資料點到叢集的分配。-每個資料點到其最接近的原型的平方距離的總和是最小值。

·對於每個資料點xn，我們引入長度為K的二進位制向量rn（K的1/K編碼），它指示資料點xn被分配給哪個K簇。

定義目標（失真測度）：

它表示每個資料點到其指定的原型k的距離的平方和。

我們的目標是找到rnk和聚類中心uk的值，以便最小化目標J。

Iterative Algorithm

定義迭代過程以最小化：

關於給定的k，將j相對於RNK（E步驟）最小化：

簡單地說，將第n個數據點Xn分配到它最接近的叢集中心。

給出給定的RNK，相對於k（m步驟）最小化J：

其中n是分配給群集K的點的數目。

集合k等於分配給群集K的所有資料點的平均值。

保證了收斂到區域性最小值（不是全域性最小值）。

舉例

 在舊的資料集上使用k-均值（k＝2）的例子，收斂步驟如下：

參考資料：

https://www.cnblogs.com/cfantaisie/archive/2011/08/20/2147075.html

matlab程式碼：

如果理解了上面的內容，寫起來一小時內就可以完成，為何不自己試一試呢。

函式：

function [data, mu, var, weight] = CreateSample(M, dim, N)
% 生成實驗樣本集，由M組正態分佈的資料構成
% % GMM模型的原理就是僅根據資料估計引數：每組正態分佈的均值、方差，
% 以及每個正態分佈函式在GMM的權重alpha。
% 在本函式中，這些引數均為隨機生成，
% 
% 輸入
%   M    : 高斯函式個數
%   dim  : 資料維數
%   N    : 資料總個數
% 返回值
%   data : dim-by-N, 每列為一個數據
%   miu  : dim-by-M, 每組樣本的均值，由本函式隨機生成
%   var  : 1-by-M, 均方差，由本函式隨機生成
%   weight: 1-by-M, 每組的權值，由本函式隨機生成
% ---------------------------------------------------- 
%
% 隨機生成不同組的方差、均值及權值
weight = rand(1,M);
weight = weight / norm(weight, 1); % 歸一化，保證總合為1
var = double(mod(int16(rand(1,M)*100),10) + 1);  % 均方差，取1~10之間，採用對角矩陣
mu = double(round(randn(dim,M)*100));            % 均值，可以有負數
 
for i = 1: M
  if i ~= M
    n(i) = floor(N*weight(i));
  else
    n(i) = N - sum(n);
  end
end
 
% 以標準高斯分佈生成樣本值，並平移到各組相應均值和方差
start = 0;
for i=1:M
  X = randn(dim, n(i));
  X = X.* var(i) + repmat(mu(:,i),1,n(i));
  data(:,(start+1):start+n(i)) = X;
  start = start + n(i);
end
save('d:\data.mat', 'data');

function [MU_pre,SIGMA_pre,Alpha_Pre,Center_Pre]=CreatePre(Gao_siNum,dimention);
% 生成隨機的MU,SIGMA和權重
% 輸入
%   Gao_siNum    : 高斯函式個數
%   dimention    : 資料維數
% 返回值
%   MU_pre : dim-Num, 每組樣本的均值，由本函式隨機生成
%   SIGMA_pre  : dim-M, 均方差，由本函式隨機生成
%   Alpha_Pre  : 1-M, 權重
%   Center_Pre : 2-M，每個點的中心 
% ---------------------------------------------------- 
% 

MU_pre=normrnd(10,5,dimention,Gao_siNum);
SIGMA_pre=normrnd(10,5,1,Gao_siNum);
Alpha_Pre=normrnd(10,5,1,Gao_siNum);
Center_Pre=normrnd(30,100,2,Gao_siNum);
% MU_pre=normrnd(rand(1),rand(1),dimention,Gao_siNum);
% SIGMA_pre=normrnd(rand(1),rand(1,1),dimention,Gao_siNum);
% Alpha_Pre=normrnd(rand(1,1),rand(1,1),1,Gao_siNum);

主程式：

close all
% %% 畫圖
% num=60;%每個集合的樣本數
% x=1:1:num;
% MU1=4;
% MU2=6;
% MU3=2;
% SIGMA=2;
% y1=normrnd(MU1,SIGMA,1,num);
% y2=normrnd(MU2,SIGMA,1,num);
% y3=normrnd(MU3,SIGMA,1,num);
% %% 畫出原影象
% figure();
% hold on
% scatter(x,y1);
% scatter(x,y2);
% scatter(x,y3);
% hold off
%% 建立生成資料並且繪圖
Gao_siNum=4;
dimention=2;
sampleNum=180;
[data, MU, SIGMA, weight] = CreateSample(Gao_siNum, dimention, sampleNum);  % 生成測試資料
draw_x=data(1,:);%x軸
draw_y=data(2,:);%y軸
figure();
scatter(draw_x,draw_y);
hold on
scatter(MU(1,:),MU(2,:));
hold off
%% 進行區分GMM_EM演算法
[MU_pre,SIGMA_pre,Alpha_Pre,Center_Pre]=CreatePre(Gao_siNum,dimention);
hold on
scatter(Center_Pre(1,:),Center_Pre(2,:));
legend('data','real center',' pre_trained center');
hold off
%% EM 迭代停止條件
maxStep=2000;
%% 初始化引數
[dim, N] = size(data);
nbStep = 0;
Epsilon = 0.0001;
distance=zeros(Gao_siNum,sampleNum);
distance_min=zeros(1,sampleNum);
distance_min_Index=zeros(1,sampleNum);
while (nbStep < 1200)
   nbStep=nbStep+1;
   %計算每個點到各自中心的衡量，需要一個dimention*sampleNum大小的矩陣來儲存
   for i=1:sampleNum
         for j=1:Gao_siNum
             %(x1-x2)^2+(y1-y2)^2
             distance(j,i)=sqrt((data(1,i)-Center_Pre(1,j))^2+(data(2,i)-Center_Pre(2,j))^2);
         end
   end
     %% E-步驟 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     for i=1:sampleNum
             distance_min(1,i)=min(distance(:,i));
             for j=1:Gao_siNum
                 if distance(j,i)==distance_min(1,i);
                     distance_min_Index(1,i)=j;%將第n個數據點Xn分配到它最接近的叢集中心。
                 end
             end
     end
     %% M-步驟 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %給出給定的RNK，相對於k（m步驟）最小化J：重新貼標籤
        %先把每個類的對應標籤找出來，然後再計算均值。
        find_dimention1= find(distance_min_Index==1); %查詢對應的類
        
        find_dimention1(1)=1;
        n=length(find_dimention1);
        Center_Pre(1,1)=sum(data(1,find_dimention1))/n;
        Center_Pre(2,1)=sum(data(2,find_dimention1))/n;
        find_dimention2= find(distance_min_Index==2); %查詢對應的類   
        
        find_dimention2(1)=1;
        n=length(find_dimention2);
        Center_Pre(1,2)=sum(data(1,find_dimention2))/n;
        Center_Pre(2,2)=sum(data(2,find_dimention2))/n;
        find_dimention3= find(distance_min_Index==3); %查詢對應的類
        
        find_dimention3(1)=1;
        n=length(find_dimention3);
        Center_Pre(1,3)=sum(data(1,find_dimention3))/n;
        Center_Pre(2,3)=sum(data(2,find_dimention3))/n;
        find_dimention4= find(distance_min_Index==4); %查詢對應的類
        n=length(find_dimention4);
        find_dimention4(1)=1;
        Center_Pre(1,4)=sum(data(1,find_dimention4))/n;
        Center_Pre(2,4)=sum(data(2,find_dimention4))/n;
%         for j=1:Gao_siNum
%             n=length(find_dimention(:,j));
%             Center_Pre(1,j)=sum(data(1,find_dimention(:,j)))/n;
%             Center_Pre(2,j)=sum(data(2,find_dimention(:,j)))/n;
%         end
%%
cost=0;
for j=1:Gao_siNum
cost=cost+sum(distance(:,j));
end
end
%%
figure();
hold on
scatter(draw_x,draw_y,'y');
scatter(MU(1,:),MU(2,:),'b');
scatter(Center_Pre(1,:),Center_Pre(2,:),'g');
legend('data','real center',' pre_trained center');
hold off

成果：