?function?varargout?=?gmm（X?K_or_centroids）
%?============================================================
%?Expectation-Maximization?iteration?implementation?of
%?Gaussian?Mixture?Model.
%
%?PX?=?GMM（X?K_OR_CENTROIDS）
%?[PX?MODEL]?=?GMM（X?K_OR_CENTROIDS）
%
%??-?X:?N-by-D?data?matrix.
%??-?K_OR_CENTROIDS:?either?K?indicating?the?number?of
%???????components?or?a?K-by-D?matrix?indicating?the
%???????choosing?of?the?initial?K?centroids.
%
%??-?PX:?N-by-K?matrix?indicating?the?probability?of?each
%???????component?generating?each?point.
%??-?MODEL:?a?structure?containing?the?parameters?for?a?GMM:
%???????MODEL.Miu:?a?K-by-D?matrix.
%???????MODEL.Sigma:?a?D-by-D-by-K?matrix.
%???????MODEL.Pi:?a?1-by-K?vector.
%?============================================================
?
????threshold?=?1e-15;
????[N?D]?=?size（X）;
?
????if?isscalar（K_or_centroids）
????????K?=?K_or_centroids;
????????%?randomly?pick?centroids
????????rndp?=?randperm（N）;
????????centroids?=?X（rndp（1:K）?:）;
????else
????????K?=?size（K_or_centroids?1）;
????????centroids?=?K_or_centroids;
????end
?
????%?initial?values
????[pMiu?pPi?pSigma]?=?init_params（）;?%初始化
?
????Lprev?=?-inf;?%inf表示正無究大，-inf表示為負無究大
????while?true
????????Px?=?calc_prob（）;
?
????????%?new?value?for?pGamma
????????pGamma?=?Px?.*?repmat（pPi?N?1）;
????????pGamma?=?pGamma?./?repmat（sum（pGamma?2）?1?K）;?%求每個樣本由第K個聚類，也叫“component“生成的概率
?
????????%?new?value?for?parameters?of?each?Component
????????Nk?=?sum（pGamma?1）;
????????pMiu?=?diag（1./Nk）?*?pGamma‘?*?X;?%重新計算每個component的均值
????????pPi?=?Nk/N;?%更新混合高斯的加權系數
????????for?kk?=?1:K?%重新計算每個component的協方差
????????????Xshift?=?X-repmat（pMiu（kk?:）?N?1）;
????????????pSigma（:?:?kk）?=?（Xshift‘?*?...
????????????????（diag（pGamma（:?kk））?*?Xshift））?/?Nk（kk）;
????????end
?
????????%?check?for?convergence
????????L?=?sum（log（Px*pPi‘））;?%求混合高斯分布的似然函數
????????if?L-Lprev?????????????break;?%似然函數收斂則退出
????????end
????????Lprev?=?L;
????end
?
????if?nargout?==?1?%如果返回是一個參數的話，那么varargout=Px;
????????varargout?=?{Px};
????else?%否則，返回[Px?model]其中model是結構體
????????model?=?[];
????????model.Miu?=?pMiu;
????????model.Sigma?=?pSigma;
????????model.Pi?=?pPi;
????????varargout?=?{Px?model};
????end
?
????function?[pMiu?pPi?pSigma]?=?init_params（）
????????pMiu?=?centroids;
????????pPi?=?zeros（1?K）;
????????pSigma?=?zeros（D?D?K）;
?
????????%?hard?assign?x?to?each?centroids
????????distmat?=?repmat（sum（X.*X?2）?1?K）?+?...?%distmat第j行的第i個元素表示第j個數據與第i個聚類點的距離，如果數據有4個，聚類2個，那么distmat就是4*2矩陣
????????????repmat（sum（pMiu.*pMiu?2）‘?N?1）?-?2*X*pMiu‘;?%sum（A，2）結果為列向量，第i個元素是第i行的求和
????????[dummy?labels]?=?min（distmat?[]?2）;?%返回列向量dummy和labels，dummy向量記錄distmat的每行的最小值，labels向量記錄每行最小值的列號，即是第幾個聚類，labels是N×1列向量，N為樣本數
?
????????for?k=1:K
????????????Xk?=?X（labels?==?k?:）;?%把標志為同一個聚類的樣本組合起來
????????????pPi（k）?=?size（Xk?1）/N;?%求混合高斯模型的加權系數，pPi為1*K的向量
????????????pSigma（