如果數據 x⁽ⁱ⁾∈Rⁿ是來自混合高斯分布，那么可用 EM 算法來擬合混合模型，但假設前提是有足夠的數據，我們能發現數據中的多高斯結構，這就要訓練數據集大小 m 遠大于數據維度 n。

如果 n>>m 呢？那就很難對數據建模，即使是一個高斯。m 個數據點只能在低維度空間，如果要把數據建模為高斯分布，用極大似然估計來估計均值和協方差：

我們會發現Σ 是奇異矩陣。例如當 m=2，n=5，如下例所示，最后Σ =0。

import numpy as np

x1=np.array([2,3,7,6,5])
x2=np.array([7,9,3,1,4])

x1.shape=(5,1)
x2.shape=(5,1)

mu=(x1+x2)/2.0  # ".0" is request

var=(np.dot(x1-mu,np.transpose(x1-mu))+np.dot(x2-mu,np.transpose(x2-mu)))/2.0

print "The covariance of the matrix is %s." % np.linalg.det(var)

View Code

當 Σ為奇異矩陣時，|∑| = 0，那么Σ^-1 就不存在，1/|Σ|^1/2 = 1/0，但這些在計算多元高斯分布的密度時都是需要的。通常情況下，除非 m 比 n 大得超過一定范圍，不然對均值和協方差的極大似然估計會非常差。當然，我們肯定希望對數據擬合一個合理的高斯模型，來獲得數據中一些有趣的協方差結構信息。那怎么辦呢？

下面我們先來看對 Σ施加兩個限制，可以使我們用少量的數據來擬合Σ，但都不能給出滿意的解決方案。接下來再討論高斯的一些特性，后面會用到，特別是如何尋找高斯的邊緣和條件分布。最后介紹因子分析模型，以及 EM 解法。

1、對Σ 的限制

如果沒有足夠的數據來擬合一個完整的方差矩陣，可以對矩陣Σ 的空間添加一些限制。例如，選擇擬合一個對角矩陣Σ，在這種情況下，協方差矩陣的極大似然估計，即對角矩陣 Σ滿足：

如何通過極大似然估計推出上式，有興趣的讀者可自行驗證。

高斯密度的輪廓是橢圓，當協方差矩陣Σ 是對角矩陣時，橢圓的主軸就與坐標軸平行。

有時我們會對協方差矩陣作進一步的限制，它不僅是對角的，而且對角線上的元素是等值的。這時，Σ=σ²I，其中σ² 是控制參數。σ² 的極大似然估計為：

這個模型所關聯的高斯密度的輪廓在二維情況下是圓，在高維情況下是球體，或超球體。

如何我們要對數據擬合一個完整的、無約束的協方差矩陣Σ，那就需要 m>>n 以使得Σ 的極大似然估計是非奇異的。如果滿足上面兩個限制中的任意一個，只要 m≥2 就能得到一個非奇異的Σ。

但是，限制Σ 為對角矩陣也意味著，建模數據的不同坐標 x_i 和 x_j 為不相關和獨立的。有時就是要獲得數據中的一些有趣的關聯，如果使用了上面對Σ 的限制，就沒法發現這種關聯。所以，下面我們介紹因子分析模型，比對角矩陣Σ 使用更多的參數，但不需要擬合完整的協方差矩陣。

2、高斯邊緣和條件分布

在介紹因子分析之前，先來討論如何尋找聯合多元高斯分布隨機變量的條件和邊緣分布。

假設有一個向量值型隨機變量：

其中 x₁ 屬于 r 維向量空間，x₁∈R^r，x₂∈R^s，x∈R^r+s，假設 x~N(μ,Σ)，其中：

這里 μ₁∈R^r，μ₂∈R^s，Σ₁₁∈R^r×r，Σ₁₂∈R^r×s 等等。因為協方差矩陣是對稱的，所以Σ₁₂=Σ_21^T。

在這種情況下，x₁ 和 x₂ 是聯合多元高斯，那么 x₁ 的邊緣分布是什么？不難看出 E[x₁]=μ₁，Cov(x₁)=E[(x₁-μ₁)(x₁-μ₁)]=Σ，根據 x₁ 和 x₂ 的聯合協方差的定義，有：

因為高斯的邊緣分布本身就是高斯，所以得 x₁ 的邊緣分布為 x₁~N(μ₁,Σ₁₁)。

那么 x₁在給定 x₂ 下的條件分布呢？由多元高斯分布的定義，可得 x₁|x₂~N(μ_1|2,Σ_1|2)，其中

下面介紹因子分析模型時，這些公式對于尋找高斯條件和邊緣分布是很有用的。

3、因子分析模型

在因子分析模型中，我們假設 (x,z) 是如下的聯合分布，其中 z∈R^k 是一個隱隨機變量：

這里模型的參數是，向量μ∈Rⁿ，矩陣Λ∈R^n×k，對角矩陣Ψ∈R^n×n，k 值通常是小于 n 的。

所以每個數據點 x⁽ⁱ⁾ 是通過如下方式生成的：

a. 從 k 維多元高斯分布中采樣得到 z⁽ⁱ⁾；

b. 通過計算μ+Λz⁽ⁱ⁾ 將 z 映射到 n 維空間；

c. 給 μ+Λz⁽ⁱ⁾增加協方差 Ψ噪聲，得到 x⁽ⁱ⁾。

同樣的，我們也可以根據下列方式來定義因子分析模型。

其中ε 和 z 是相互獨立的。

現在來仔細看下這個模型定義了一個什么樣的分布。隨機變量 z 和 x 有一個聯合高斯分布。

現在來尋找μ_zx 和Σ。

我們知道 E[z]=0，因為 z~N(0,I)，同時有

放到一起可得：

下一步，為了尋找 Σ，需要計算Σ_zz=E[(z-E(z))(z-E(z))^T]，Σ_zx=E[(z-E(z))(x-E(x))^T]，Σ_xx=E[(x-E(x))(x-E(x))^T]。

因為 z~N(0,I)，易見Σ_zz=Cov(z)=I，也可得 Σ_zx=

第一步是因為 E[z]=0，x=μ+Λz+ε，E[x]=μ。

第二步消去μ。

第三步是因為 E[zz^T]=Cov(z)=I，E[zε^T]=E[z]E[ε^T]=0（因為 z 和 ε 是相互獨立的，所以乘積的期望也是期望的乘積）。類似地可得Σ_xx：

集中以上結果可得：

所以，可以看到 x 的邊緣分布 x~N(μ,ΛΛ^T+Ψ)。所以，給定訓練集 {x⁽ⁱ⁾;i=1,2,...,m}，我們能寫下參數的 log 似然估計為：

為執行極大似然估計，我們將以這些參數來最大化該式，直接最大化該式很難，沒有封閉形式下的算法可解出該式。所以，我們將使用 EM 算法。

4、EM 解因子分析

E-step 的推導很容易，需要計算 Q_i(z⁽ⁱ⁾)=p(z⁽ⁱ⁾|x⁽ⁱ⁾;μ,Λ,Ψ)。結合上面得到的兩個結果：

可得

其中

用這些定義可得：

現在開始 M-step，這里需要最大化

以μ，Λ 和Ψ 為參數。我們這里只推導以Λ 為參數的最優化，μ 和Ψ 的更新推導留給讀者。

上式可簡化為：

下標 z⁽ⁱ⁾~Q_i 表示期望是指服從 Q_i 分布的 z⁽ⁱ⁾。下面將忽略這個下標，不會存在模糊的風險。去掉與參數Λ 不相關的部分，需要最大化的是：

只有最后一項依賴于Λ，求導，使用事實 tra=a(a∈R)，trAB=trBA，▽_AtrABA^TC=CAB+C^TAB，有：

設為 0 并簡化，得：

解出Λ，得：

從 Q 是均值為 μ_z|x協方差為 Σ_z|x的高斯的定義可得：

最后一式是因為，對于隨機變量 Y，Cov(Y)=E[YY^T]-E[Y]E[Y]^T，所以，E[YY^T]=E[Y]E[Y]^T+Cov(Y)。

代入該式到上面的Λ，可得：

注意到上式右邊有個 Σ_z|x，這是后驗概率 p(z⁽ⁱ⁾|x⁽ⁱ⁾) 的協方差。

我們也能找到μ 和Φ 的最優化參數。

注意到當參數改變時，μ 不是發生改變，不像Λ 依賴于 Q_i(z⁽ⁱ⁾)=p(z⁽ⁱ⁾|x⁽ⁱ⁾;μ,Λ,Φ)，所以，它只需要開始計算一次，算法運行時不需要再更新。同樣，對角Ψ 為

設置Ψ_ii=Φ_ii。對角矩陣Ψ 只包含對角元素Φ。

參考資料：

1、http://cs229.stanford.edu/notes/cs229-notes9.pdf

2、http://blog.csdn.net/stdcoutzyx/article/details/37559995

3、http://www.cnblogs.com/jerrylead/archive/2011/05/11/2043317.html

總結

以上是生活随笔為你收集整理的机器学习笔记—因子分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。