轉載請註明出處，該文章的官方來源：
高斯混合分佈 | Teaching ML
現有的高斯模型有單高斯模型（SGM）和高斯混合模型（GMM）兩種。從幾何上講，單高斯分佈模型在二維空間上近似於橢圓，在三維空間上近似於橢球。 在很多情況下，屬於同一類別的樣本點並不滿足“橢圓”分佈的特性，所以我們需要引入混合高斯模型來解決這種情況。
1 單高斯模型
多維變數X服從高斯分佈時，它的機率密度函式PDF定義如下：
在上述定義中，x是維數為D的樣本向量，mu是模型期望，sigma是模型協方差。對於單高斯模型，可以明確訓練樣本是否屬於該高斯模型，所以我們經常將mu用訓練樣本的均值代替，將sigma用訓練樣本的協方差代替。 假設訓練樣本屬於類別C，那麼上面的定義可以修改為下面的形式：
這個公式表示樣本屬於類別C的機率。我們可以根據定義的機率閾值來判斷樣本是否屬於某個類別。
2 高斯混合模型
高斯混合模型，顧名思義，就是資料可以看作是從多個高斯分佈中生成出來的。從中心極限定理可以看出，高斯分佈這個假設其實是比較合理的。 為什麼我們要假設資料是由若干個高斯分佈組合而成的，而不假設是其他分佈呢？實際上不管是什麼分佈，只K取得足夠大，這個XX Mixture Model就會變得足夠複雜，就可以用來逼近任意連續的機率密度分佈。只是因為高斯函式具有良好的計算效能，所GMM被廣泛地應用。
每個GMM由K個高斯分佈組成，每個高斯分佈稱為一個元件（Component），這些元件線性加成在一起就組成了GMM的機率密度函式**（1）**：
根據上面的式子，如果我們要從GMM分佈中隨機地取一個點，需要兩步：
隨機地在這K個元件之中選一個，每個元件被選中的機率實際上就是它的係數pi_k；
選中了元件之後，再單獨地考慮從這個元件的分佈中選取一個點。
怎樣用GMM來做聚類呢？其實很簡單，現在我們有了資料，假定它們是由GMM生成出來的，那麼我們只要根據資料推出GMM的機率分佈來就可以了，然後GMM的K個元件實際上就對應了K個聚類了。 在已知機率密度函式的情況下，要估計其中的引數的過程被稱作“引數估計”。
我們可以利用最大似然估計來確定這些引數，GMM的似然函式**（2）**如下：
可以用EM演算法來求解這些引數。EM演算法求解的過程如下：
1
E-步
。求資料點由各個元件生成的機率（並不是每個元件被選中的機率）。對於每個資料
來說，它由第k個元件生成的機率為公式**（3）**：
在上面的機率公式中，我們假定mu和sigma均是已知的，它們的值來自於初始化值或者上一次迭代。
2
M-步
。估計每個元件的引數。由於每個元件都是一個標準的高斯分佈，可以很容易分佈求出最大似然所對應的引數值，分別如下公式**（4）
,
（5）
,
（6）
,
（7）**：
3 原始碼分析
3。1 例項
在分析原始碼前，我們還是先看看高斯混合模型如何使用。
import org。apache。spark。mllib。clustering。GaussianMixture
import org。apache。spark。mllib。clustering。GaussianMixtureModel
import org。apache。spark。mllib。linalg。Vectors
// 載入資料
val data = sc。textFile（“data/mllib/gmm_data。txt”）
val parsedData = data。map（s => Vectors。dense（s。trim。split（‘ ’）。map（_。toDouble）））。cache（）
// 使用高斯混合模型聚類
val gmm = new GaussianMixture（）。setK（2）。run（parsedData）
// 儲存和載入模型
gmm。save（sc， “myGMMModel”）
val sameModel = GaussianMixtureModel。load（sc， “myGMMModel”）
// 列印引數
for （i <- 0 until gmm。k） {
println（“weight=%f\nmu=%s\nsigma=\n%s\n” format
（gmm。weights（i）， gmm。gaussians（i）。mu， gmm。gaussians（i）。sigma））
}
由上面的程式碼我們可以知道，使用高斯混合模型聚類使用到了GaussianMixture類中的run方法。下面我們直接進入run方法，分析它的實現。
3。2 高斯混合模型的實現
3。2。1 初始化
在run方法中，程式所做的第一步就是初始化權重（上文中介紹的pi）及其相對應的高斯分佈。
val （weights， gaussians） = initialModel match {
case Some（gmm） => （gmm。weights， gmm。gaussians）
case None => {
val samples = breezeData。takeSample（withReplacement = true， k * nSamples， seed）
（Array。fill（k）（1。0 / k）， Array。tabulate（k） { i =>
val slice = samples。view（i * nSamples， （i + 1） * nSamples）
new MultivariateGaussian（vectorMean（slice）， initCovariance（slice））
}）
}
}
在上面的程式碼中，當initialModel為空時，用所有值均為1。0/k的陣列初始化權重，用值為MultivariateGaussian物件的陣列初始化所有的高斯分佈（即上文中提到的元件）。 每一個MultivariateGaussian物件都由從資料集中抽樣的子集計算而來。這裡用樣本資料的均值和方差初始化MultivariateGaussian的mu和sigma。
//計算均值
private def vectorMean（x： IndexedSeq［BV［Double］］）： BDV［Double］ = {
val v = BDV。zeros［Double］（x（0）。length）
x。foreach（xi => v += xi）
v / x。length。toDouble
}
//計算方差
private def initCovariance（x： IndexedSeq［BV［Double］］）： BreezeMatrix［Double］ = {
val mu = vectorMean（x）
val ss = BDV。zeros［Double］（x（0）。length）
x。foreach（xi => ss += （xi - mu） ：^ 2。0）
diag（ss / x。length。toDouble）
}
3。2。2 EM演算法求引數
初始化後，就可以使用EM演算法迭代求似然函式中的引數。迭代結束的條件是迭代次數達到了我們設定的次數或者兩次迭代計算的對數似然值之差小於閾值。
while （iter < maxIterations && math。abs（llh-llhp） > convergenceTol）
在迭代內部，就可以按照E-步和M-步來更新引數了。
E-步
：更新引數gamma
val compute = sc。broadcast（ExpectationSum。add（weights， gaussians）_）
val sums = breezeData。aggregate（ExpectationSum。zero（k， d））（compute。value， _ += _）
我們先要了解ExpectationSum以及add方法的實現。
private class ExpectationSum（
var logLikelihood： Double，
val weights： Array［Double］，
val means： Array［BDV［Double］］，
val sigmas： Array［BreezeMatrix［Double］］） extends Serializable
ExpectationSum是一個聚合類，它表示部分期望結果：主要包含對數似然值，權重值（第二章中介紹的pi），均值，方差。add方法的實現如下：
def add（ weights： Array［Double］，dists： Array［MultivariateGaussian］）
（sums： ExpectationSum， x： BV［Double］）： ExpectationSum = {
val p = weights。zip（dists）。map {
//計算pi_i * N（x）
case （weight， dist） => MLUtils。EPSILON + weight * dist。pdf（x）
}
val pSum = p。sum
sums。logLikelihood += math。log（pSum）
var i = 0
while （i < sums。k） {
p（i） /= pSum
sums。weights（i） += p（i）
sums。means（i） += x * p（i）
//A ：= alpha * x * x^T^ + A
BLAS。syr（p（i）， Vectors。fromBreeze（x），
Matrices。fromBreeze（sums。sigmas（i））。asInstanceOf［DenseMatrix］）
i = i + 1
}
sums
}
從上面的實現我們可以看出，最終，logLikelihood表示公式**（2）
中的對數似然。p和weights分別表示公式
（3）
中的gamma和pi，means表示公式
（6）
中的求和部分，sigmas表示公式
（7）**中的求和部分。
呼叫RDD的aggregate方法，我們可以基於所有給定資料計算上面的值。利用計算的這些新值，我們可以在M-步中更新mu和sigma。
M-步
：更新引數mu和sigma
var i = 0
while （i < k） {
val （weight， gaussian） =
updateWeightsAndGaussians（sums。means（i）， sums。sigmas（i）， sums。weights（i）， sumWeights）
weights（i） = weight
gaussians（i） = gaussian
i = i + 1
}
private def updateWeightsAndGaussians（
mean： BDV［Double］，
sigma： BreezeMatrix［Double］，
weight： Double，
sumWeights： Double）： （Double， MultivariateGaussian） = {
// mean/weight
val mu = （mean /= weight）
// -weight * mu * mut +sigma
BLAS。syr（-weight， Vectors。fromBreeze（mu），
Matrices。fromBreeze（sigma）。asInstanceOf［DenseMatrix］）
val newWeight = weight / sumWeights
val newGaussian = new MultivariateGaussian（mu， sigma / weight）
（newWeight， newGaussian）
}
基於
E-步
計算出來的值，根據公式**（6）
，我們可以透過(mean /= weight)來更新mu；根據公式
（7）
，我們可以透過BLAS.syr()來更新sigma；同時，根據公式
（5）**， 我們可以透過weight / sumWeights來計算pi。
迭代執行以上的
E-步
和
M-步
，到達一定的迭代數或者對數似然值變化較小後，我們停止迭代。這時就可以獲得聚類後的引數了。
3。3 多元高斯模型中相關方法介紹
在上面的求參程式碼中，我們用到了MultivariateGaussian以及MultivariateGaussian中的部分方法，如pdf。MultivariateGaussian定義如下：
class MultivariateGaussian @Since（“1。3。0”） （
@Since（“1。3。0”） val mu： Vector，
@Since（“1。3。0”） val sigma： Matrix） extends Serializable
MultivariateGaussian包含一個向量mu和一個矩陣sigma，分別表示期望和協方差。MultivariateGaussian最重要的方法是pdf，顧名思義就是計算給定資料的機率密度函式。它的實現如下：
private［mllib］ def pdf（x： BV［Double］）： Double = {
math。exp（logpdf（x））
}
private［mllib］ def logpdf（x： BV［Double］）： Double = {
val delta = x - breezeMu
val v = rootSigmaInv * delta
u + v。t * v * -0。5
}
上面的rootSigmaInv和u透過方法calculateCovarianceConstants計算。根據公式**（1）**，這個機率密度函式的計算需要計算sigma的行列式以及逆。
sigma = U * D * U。t
inv（Sigma） = U * inv（D） * U。t = （D^{-1/2}^ * U。t）。t * （D^{-1/2}^ * U。t）
-0。5 * （x-mu）。t * inv（Sigma） * （x-mu） = -0。5 * norm（D^{-1/2}^ * U。t * （x-mu））^2^
這裡，U和D是奇異值分解得到的子矩陣。calculateCovarianceConstants具體的實現程式碼如下：
private def calculateCovarianceConstants： （DBM［Double］， Double） = {
val eigSym。EigSym（d， u） = eigSym（sigma。toBreeze。toDenseMatrix） // sigma = u * diag（d） * u。t
val tol = MLUtils。EPSILON * max（d） * d。length
try {
//所有非0奇異值的對數和
val logPseudoDetSigma = d。activeValuesIterator。filter（_ > tol）。map（math。log）。sum
//透過求非負值的倒數平方根，計算奇異值對角矩陣的根偽逆矩陣
val pinvS = diag（new DBV（d。map（v => if （v > tol） math。sqrt（1。0 / v） else 0。0）。toArray））
（pinvS * u。t， -0。5 * （mu。size * math。log（2。0 * math。Pi） + logPseudoDetSigma））
} catch {
case uex： UnsupportedOperationException =>
throw new IllegalArgumentException（“Covariance matrix has no non-zero singular values”）
}
}
上面的程式碼中，eigSym用於分解sigma矩陣。
4 參考文獻
【1】漫談 Clustering （3）： Gaussian Mixture Model
