Task1-概率PCA

PPCA的模型形式
- 1、极大似然估计
- - 1.1 μ\muμ的ML估计
  - 1.2 W\textrm{W}W的ML估计
  - 1.3 σ2\sigma^{2}σ2的ML估计
- 2、PPCA的EM算法
- - 2.1 计算后验分布
  - 2.2 E Step
  - 2.3 M Step
  - - 2.3.1 W的估计
    - 2.3.2 τ\tauτ的估计
- 3 模拟分析
- - 3.1 极大似然估计R语言代码
  - 3.2 极大似然估计matlab代码
  - 3.3 产生随机数模拟
  - 3.4 EM算法R代码
  - 3.5 EM算法matlab代码
  - 3.6 模拟（matlab代码）
- 4 总结

（以下内容是自己在完成老师任务时经过查阅资料自行整理的笔记，在变量的字母符号设置上前后有所差异，但不影响阅读，如有推导错误，请联系我，此笔记仅用于学习，如需转载，请注明来源，谢谢！）

PPCA的模型形式

假设d维的观测变量t对应q维的潜变量(latent)x。最常见的模型形式为线性模型：
t=Wx+μ+ϵ(1)\textrm{t}=\textrm{W}\textrm{x}+\mu+\epsilon \tag{1}t=Wx+μ+ϵ(1)
其中，W\textrm{W}W是 d×qd \times qd×q的矩阵，μ\muμ表示模型具有非零均值。q<dq<dq<d从而达到降维简化数据的目的。
假设噪声ϵ∼N(0,σ2Id)\epsilon\sim N(\textbf{0},\sigma^{2}\textbf{I}_d)ϵ∼N(0,σ2Id)，潜变量x\textbf{x}x的先验分布也是高斯的，且x∼N(0,Iq)\textbf{x}\sim N(\textbf{0},\textbf{I}_q)x∼N(0,Iq),则在给定x的情况下t的条件分布为：
t∣x∼N(Wx+μ,σ2Id)(2)\textbf{t}\mid \textbf{x} \sim N(\textbf{W}\textbf{x}+\mu,\sigma^{2}\textbf{I}_d) \tag{2}t∣x∼N(Wx+μ,σ2Id)(2)
则观测数据t的边际分布容易通过积分取得，如：
p(t)=∫p(t∣x)p(x)dxp(t)=\int p(t\mid x)p(x)dxp(t)=∫p(t∣x)p(x)dx
很明显，t的观测分布也是高斯的，假设t∼N(μd,Cd)t\sim N(\mu_d,\textbf{C}_d)t∼N(μd,Cd),则
E(x)=E(Wx+μ+ϵ)=μE(\textbf{x})=E(\textrm{W}\textrm{x}+\mu+\epsilon)=\mu E(x)=E(Wx+μ+ϵ)=μ
C=Cov(t)=cov(Wx+μ+ϵ,Wx+μ+ϵ)=cov(Wx+ϵ,Wx+ϵ)=WWT+σ2Id(3)\begin{aligned} C=Cov(t)&=cov(\textrm{W}\textrm{x}+\mu+\epsilon ,\textrm{W}\textrm{x}+\mu+\epsilon ) \\ &=cov(\textrm{W}\textrm{x}+\epsilon ,\textrm{W}\textrm{x}+\epsilon )\\ &=\textbf{W}\textbf{W}^{T}+\sigma^{2}\textbf{I}_d \end{aligned} \tag{3}C=Cov(t)=cov(Wx+μ+ϵ,Wx+μ+ϵ)=cov(Wx+ϵ,Wx+ϵ)=WWT+σ2Id(3)
则观测数据t相关的似然函数为：
l=∏n=1Np(tn∣C)=∏n=1N∣2πC∣−12e−12(tn−μ)′C−1(tn−μ)=(2π)−Nd2∣C∣−N2e−12∑n=1N(tn−μ)′C−1(tn−μ)\begin{aligned} l&=\prod_{n=1}^Np(t_n|C)\\ &=\prod_{n=1}^N|2\pi C|^{-\frac{1}{2}}e^{-\frac{1}{2}(t_n-\mu)'C^{-1}(t_n-\mu)} \\ &=(2\pi)^{-\frac{Nd}{2}}|C|^{-\frac{N}{2}}e^{-\frac{1}{2}\sum_{n=1}^{N}(t_n-\mu)'C^{-1}(t_n-\mu)} \end{aligned}l=n=1∏Np(tn∣C)=n=1∏N∣2πC∣−21e−21(tn−μ)′C−1(tn−μ)=(2π)−2Nd∣C∣−2Ne−21∑n=1N(tn−μ)′C−1(tn−μ)
则t的对数似然函数为：
L=log(l)=log(∏n=1Np(tn∣C))=−N2{dlog⁡(2π)+log∣C∣+1N∑n=1N(tn−μ)′C−1(tn−μ)}=−N2{dlog⁡(2π)+log∣C∣+tr(C−1S)}(4)\begin{aligned} L&=log(l)=log(\prod_{n=1}^Np(t_n|C))\\ &=-\frac{N}{2} \{d\log(2\pi)+log|C|+\frac{1}{N}\sum_{n=1}^{N}(t_n-\mu)'C^{-1}(t_n-\mu)\} \\ &=-\frac{N}{2} \{d\log(2\pi)+log|C|+tr(C^{-1}S)\} \\ \end{aligned} \tag{4}L=log(l)=log(n=1∏Np(tn∣C))=−2N{dlog(2π)+log∣C∣+N1n=1∑N(tn−μ)′C−1(tn−μ)}=−2N{dlog(2π)+log∣C∣+tr(C−1S)}(4)
其中,
S=1N∑n=1N(tn−μ)(tn−μ)′S=\frac{1}{N}\sum_{n=1}^{N}(t_n-\mu)(t_n-\mu)' S=N1n=1∑N(tn−μ)(tn−μ)′
1N∑n=1N(tn−μ)′C−1(tn−μ)=1N∑n=1Ntr(tn−μ)′C−1(tn−μ)=1N∑n=1Ntr(C−1(tn−μ)(tn−μ)′)=tr(C−11N∑n=1N(tn−μ)(tn−μ)′)=tr(C−1S)\begin{aligned} \frac{1}{N}\sum_{n=1}^{N}(t_n-\mu)'C^{-1}(t_n-\mu)&=\frac{1}{N}\sum_{n=1}^{N}tr(t_n-\mu)'C^{-1}(t_n-\mu) \\ &=\frac{1}{N}\sum_{n=1}^{N}tr(C^{-1}(t_n-\mu)(t_n-\mu)') \\ &=tr(C^{-1}\frac{1}{N}\sum_{n=1}^{N}(t_n-\mu)(t_n-\mu)') \\ &=tr(C^{-1}S) \end{aligned}N1n=1∑N(tn−μ)′C−1(tn−μ)=N1n=1∑Ntr(tn−μ)′C−1(tn−μ)=N1n=1∑Ntr(C−1(tn−μ)(tn−μ)′)=tr(C−1N1n=1∑N(tn−μ)(tn−μ)′)=tr(C−1S)

通过以上计算发现，观测数据t的对数似然函数与C和μ\muμ有关，而C与因子载荷矩阵W和噪声方差σ2\sigma^2σ2有关。因此，对数似然函数中的参数有μ,W和σ2\mu,\textbf{W}和\sigma^2μ,W和σ2,可通过极大似然估计进行求解。

1、极大似然估计

1.1 μ\muμ的ML估计

将对数似然函数省略与 μ\muμ无关的项，得到
L=log(l)=log(∏n=1Np(tn∣C))=−N2{dlog⁡(2π)+log∣C∣+1N∑n=1N(tn−μ)′C−1(tn−μ)}∝∑n=1N(tn−μ)′C−1(tn−μ)∝∑n=1N(−tn′C−1μ−μ′C−1tn+μ′C−1μ})(5)\begin{aligned} L&=log(l)=log(\prod_{n=1}^Np(t_n|C))\\ &=-\frac{N}{2} \{d\log(2\pi)+log|C|+\frac{1}{N}\sum_{n=1}^{N}(t_n-\mu)'C^{-1}(t_n-\mu)\} \\ &\propto \sum_{n=1}^{N}(t_n-\mu)'C^{-1}(t_n-\mu)\\ &\propto \sum_{n=1}^{N}(-t_{n}^{'}C^{-1}\mu-\mu^{'} C^{-1}t_n+\mu'C^{-1}\mu\})\\ \end{aligned} \tag{5} L=log(l)=log(n=1∏Np(tn∣C))=−2N{dlog(2π)+log∣C∣+N1n=1∑N(tn−μ)′C−1(tn−μ)}∝n=1∑N(tn−μ)′C−1(tn−μ)∝n=1∑N(−tn′C−1μ−μ′C−1tn+μ′C−1μ})(5)
基于两个公式：
∂xTa∂x=∂aTx∂x=a\frac{\partial \textrm{x}^{T}\textrm{a}}{\partial \textrm{x}}=\frac{\partial \textrm{a}^{T}\textrm{x}}{\partial \textrm{x}}=\textrm{a} ∂x∂xTa=∂x∂aTx=a
∂xTAx∂x=∂xT∂xAx+∂(Ax)T∂xx=(A+A′)x=2Ax\frac{\partial \textrm{x}^{T}\textrm{A}\textrm{x}}{\partial \textrm{x}}=\frac{\partial \textrm{x}^{T}}{\partial \textrm{x}}\textrm{A}\textrm{x} +\frac{\partial(\textrm{Ax})^{T}}{\partial\textrm{x}}\textrm{x}=(\textrm{A}+\textrm{A}')\textrm{x}=2\textrm{A}\textrm{x} ∂x∂xTAx=∂x∂xTAx+∂x∂(Ax)Tx=(A+A′)x=2Ax
因此，对对数似然函数关于μ\muμ求偏导得到如下结果：
∂L∂μ=∑n=1N(−2C−1tn+2C−1μ)=0C−1∑n=1N(2μ−2tn)=0∑n=1N(μ−tn)=0μ^=1N∑n=1Ntn=tˉ(6)\frac{\partial L}{\partial \mu}=\sum_{n=1}^{N}(-2C^{-1}t_n+2C^{-1}\mu)=0 \\ C^{-1}\sum_{n=1}^{N}(2\mu-2t_n)=0 \\ \sum_{n=1}^{N}(\mu-t_n)=0 \\ \hat{\mu} = \frac{1}{N}\sum_{n=1}^{N}t_n=\bar{t} \tag{6} ∂μ∂L=n=1∑N(−2C−1tn+2C−1μ)=0C−1n=1∑N(2μ−2tn)=0n=1∑N(μ−tn)=0μ^=N1n=1∑Ntn=tˉ(6)
因此，参数μ\muμ的ML估计即是观测数据的均值向量tˉ\bar{t}tˉ。

1.2 W\textrm{W}W的ML估计

对公式(4)关于W\textrm{W}W求导，但是(4)主要与C有关，因此先对(4)式关于C求导。注意：此处要使几个矩阵微分公式和矩阵公式：
dlog⁡(∣C∣)=1∣C∣d∣C∣=1∣C∣∣C∣tr(C−1dC)=tr(C−1dC)(7)d\log(|C|)=\frac{1}{|C|}d|C|=\frac{1}{|C|}|C|tr(C^{-1}dC)=tr(C^{-1}dC) \tag{7} dlog(∣C∣)=∣C∣1d∣C∣=∣C∣1∣C∣tr(C−1dC)=tr(C−1dC)(7)
dC−1=−C−1dCC−1(8)dC^{-1}=-C^{-1}dCC^{-1} \tag{8} dC−1=−C−1dCC−1(8)
dC=(dW)W′+W(dW′)(9)dC=(d\textrm{W})\textrm{W}'+\textrm{W}(d\textrm{W}') \tag{9} dC=(dW)W′+W(dW′)(9)
tr(AB)=tr(BA)tr(AB)=tr(BA)tr(AB)=tr(BA)
故根据以上结论，对W\textrm{W}W微分得：
dL=dlog⁡(∣C∣)+dtr(C−1S)=tr(C−1dC)+tr(d(C−1)S)=tr(C−1((dW)W′+W(dW′)))−tr(C−1((dW)W′+W(dW′))C−1S)=2tr(C−1WdW′)−2tr(C−1SC−1WdW′)=2tr(C−1WdW′−C−1SC−1WdW′)=0(10)\begin{aligned} dL&=d\log(|C|) +dtr(C^{-1}S) \\ &=tr(C^{-1}dC)+tr(d(C^{-1})S) \\ &=tr(C^{-1}((d\textrm{W})\textrm{W}'+\textrm{W}(d\textrm{W}')))-tr(C^{-1}((d\textrm{W})\textrm{W}'+\textrm{W}(d\textrm{W}') )C^{-1}S) \\ &=2tr(C^{-1}WdW')-2tr(C^{-1}SC^{-1}WdW') \\ &=2tr(C^{-1}WdW'-C^{-1}SC^{-1}WdW')=0 \end{aligned} \tag{10} dL=dlog(∣C∣)+dtr(C−1S)=tr(C−1dC)+tr(d(C−1)S)=tr(C−1((dW)W′+W(dW′)))−tr(C−1((dW)W′+W(dW′))C−1S)=2tr(C−1WdW′)−2tr(C−1SC−1WdW′)=2tr(C−1WdW′−C−1SC−1WdW′)=0(10)
由(10)式可知，
C−1W=C−1SC−1WC^{-1}W=C^{-1}SC^{-1}WC−1W=C−1SC−1W
W=SC−1W(11)W=SC^{-1}W \tag{11}W=SC−1W(11)
由
(I+AB)−1A=A(I+BA)−1(I+AB)^{-1}A=A(I+BA)^{-1} (I+AB)−1A=A(I+BA)−1
将C带入（11）式，得
W=S(WW′+τId)−1W=SW(W′W+τIq)−1(12)W=S(WW'+\tau I_{d})^{-1}W=SW(W'W+\tau I_{q})^{-1} \tag{12}W=S(WW′+τId)−1W=SW(W′W+τIq)−1(12)
对W进行奇异值分解，W=ULV′W=ULV'W=ULV′,其中，U是d×qd\times qd×q的列正交矩阵，L是q×qq \times qq×q的奇异值对角矩阵，V是q×qq \times qq×q的正交矩阵。将此式带入上式得：
ULV′=SULV′((ULV′)′ULV′+τIq)−1=SULV′(VL2V′+τIq)−1\begin{aligned} ULV'&=SULV'((ULV')'ULV'+\tau I_q)^{-1} \\ &=SULV'(VL^2V'+\tau I_q)^{-1}\\ \end{aligned} ULV′=SULV′((ULV′)′ULV′+τIq)−1=SULV′(VL2V′+τIq)−1
两边同时乘以VL2V′+τIqVL^2V'+\tau I_qVL2V′+τIq,得到
ULV′(VL2V′+τIq)=SULV′UL(L2V′+τV′)=SULV′UL(L2+τIq)=SUL(两边同时乘以V)U(L2+τIq)L=SUL(13)\begin{aligned} ULV'(VL^2V'+\tau I_q)&=SULV' \\ UL(L^2V'+\tau V')&=SULV' \\ UL(L^2+\tau I_q)&=SUL(两边同时乘以V) \\ U(L^2+\tau I_q)L&=SUL \tag{13} \end{aligned} ULV′(VL2V′+τIq)UL(L2V′+τV′)UL(L2+τIq)U(L2+τIq)L=SULV′=SULV′=SUL(两边同时乘以V)=SUL(13)
L是W的奇异值，当ljl_jlj不为0时，上式表示Suj=(τ+lj2)ujSu_j=(\tau+l_j^{2})u_jSuj=(τ+lj2)uj，因此，u的每一列是S的特征向量，所对应的特征值λj=τ+lj2\lambda_j=\tau+l_j^{2}λj=τ+lj2,因此，lj=(λj−τ)1/2l_j=(\lambda_j-\tau)^{1/2}lj=(λj−τ)1/2。当lj=0l_j=0lj=0时，uju_juj是任意的。
故W所有的潜在解为：
W=Uq(Kq−τIq)1/2R(14)W=U_q(K_q-\tau I_q)^{1/2}R \tag{14}W=Uq(Kq−τIq)1/2R(14)
R为任意的正交矩阵，一般取单位阵I即可。KqK_qKq是对角矩阵，是
kj={λjuj对应的特征值τ其他k_j=\begin{cases} \lambda_j & u_j对应的特征值 \\ \tau &其他 \end{cases} kj={λjτuj对应的特征值其他
但是上式中依然包含参数τ\tauτ,需要将τ\tauτ估计出来，再将其带入即可。

1.3 σ2\sigma^{2}σ2的ML估计

令σ2=τ\sigma^{2}=\tauσ2=τ,则对（4）式关于τ\tauτ微分可得：
dLτ=dlog⁡(∣C∣)+dtr(C−1S)=tr(C−1dτ)−tr(C−1dτC−1S)=tr((C−1−C−1SC−1)dτ)=tr(C−1−C−1SC−1)dτ(因为dτ是个标量，可以提取出来)(15)\begin{aligned} dL_{\tau}&=d\log(|C|) +dtr(C^{-1}S) \\ &=tr(C^{-1}d\tau)-tr(C^{-1}d\tau C^{-1}S) \\ &=tr((C^{-1}-C^{-1}SC^{-1})d\tau)\\ &=tr(C^{-1}-C^{-1}SC^{-1})d\tau(因为d\tau是个标量，可以提取出来) \end{aligned} \tag{15} dLτ=dlog(∣C∣)+dtr(C−1S)=tr(C−1dτ)−tr(C−1dτC−1S)=tr((C−1−C−1SC−1)dτ)=tr(C−1−C−1SC−1)dτ(因为dτ是个标量，可以提取出来)(15)
则有(15)推出，
dLτdτ=tr(C−1−C−1SC−1)=0(16)\frac{dL_{\tau}}{d\tau}=tr(C^{-1}-C^{-1}SC^{-1})=0 \tag{16} dτdLτ=tr(C−1−C−1SC−1)=0(16)
注意下面公式：
(A+BD−1C)−1=A−1−A−1B(D+CA−1B)−1CA−1(17)(A+BD^{-1}C)^{-1}=A^{-1}-A^{-1}B(D+CA^{-1}B)^{-1}CA^{-1} \tag{17}(A+BD−1C)−1=A−1−A−1B(D+CA−1B)−1CA−1(17)
(I+AB)−1A=A(I+BA)−1(18)(I+AB)^{-1}A=A(I+BA)^{-1} \tag{18}(I+AB)−1A=A(I+BA)−1(18)
则
C−1=τ−1Id−τ−1WM−1W′=τ−1(Id−WM−1W′)(19)C^{-1}=\tau^{-1}I_d-\tau^{-1}WM^{-1}W'=\tau^{-1}(I_d-WM^{-1}W') \tag{19}C−1=τ−1Id−τ−1WM−1W′=τ−1(Id−WM−1W′)(19),其中矩阵M=W′W+τIqM=W'W+\tau I_qM=W′W+τIq
由上文知道，SW=WMSW=WMSW=WM,所以
C−1−C−1SC−1=C−1(Id−SC−1)=C−1(Id−τ−1S−τ−1WW′)=τ−1C−1(τId−S−WW′)=τ−1C−1(C−S)=τ−1(Id−C−1S)\begin{aligned} C^{-1}-C^{-1}SC^{-1}&=C^{-1}(I_d-S C^{-1}) \\ &=C^{-1}(I_d -\tau^{-1} S -\tau^{-1}WW') \\ &=\tau^{-1}C^{-1}(\tau I_d - S -WW') \\ &=\tau^{-1}C^{-1}(C-S) \\ &=\tau^{-1}(I_d -C^{-1}S) \end{aligned} C−1−C−1SC−1=C−1(Id−SC−1)=C−1(Id−τ−1S−τ−1WW′)=τ−1C−1(τId−S−WW′)=τ−1C−1(C−S)=τ−1(Id−C−1S)
tr(C−1−C−1SC−1)=tr(τ−1(Id−C−1S))=τ−1tr(Id−τ−1S+τ−1WM−1W′S)=τ−1tr(Id−τ−1UdKdUd+τ−1SWM−1W′)=τ−2[tr(τId)−tr(UdKdUd′)+tr(WW′)](令W=UqLV′)=τ−2[τd−tr(Kd)+tr(Lq2)]\begin{aligned} tr(C^{-1}-C^{-1}SC^{-1})&=tr(\tau^{-1}(I_d -C^{-1}S)) \\ &=\tau^{-1}tr(I_d - \tau^{-1}S +\tau^{-1}WM^{-1}W'S) \\ &=\tau^{-1}tr(I_d - \tau^{-1}U_d K_d U_d +\tau^{-1}SWM^{-1}W') \\ &=\tau^{-2}[tr(\tau I_d )- tr(U_d K_d U_d')+tr(WW') ] (令W=U_qLV')\\ &=\tau^{-2}[\tau d-tr(K_d)+tr(L_q^2)] \end{aligned} tr(C−1−C−1SC−1)=tr(τ−1(Id−C−1S))=τ−1tr(Id−τ−1S+τ−1WM−1W′S)=τ−1tr(Id−τ−1UdKdUd+τ−1SWM−1W′)=τ−2[tr(τId)−tr(UdKdUd′)+tr(WW′)](令W=UqLV′)=τ−2[τd−tr(Kd)+tr(Lq2)]
由上述式子可得，
tr(C−1−C−1SC−1)=τ−2[τd−Σi=1dλi+Σi=1q(λi−τ)=0tr(C^{-1}-C^{-1}SC^{-1})=\tau^{-2}[\tau d - \Sigma_{i=1}^{d}\lambda_i +\Sigma_{i=1}^{q}(\lambda_i - \tau)=0tr(C−1−C−1SC−1)=τ−2[τd−Σi=1dλi+Σi=1q(λi−τ)=0
τ=Σi=q+1dλid−q\tau = \frac{\Sigma_{i=q+1}^{d}\lambda_i}{d-q}τ=d−qΣi=q+1dλi

当τ\tauτ被估计出来后，则W也可以估计出来了，至此极大似然估计完成！

2、PPCA的EM算法

用EM算法迭代最大化PPCA的似然函数求解参数值，μ\muμ由上文已知容易求得，不用迭代，因此只对W和σ2\sigma^2σ2进行迭代，为了符号上的方便，假设参数为θ=(W,σ2)\theta=(W,\sigma^2)θ=(W,σ2)（这样写不太准确，只是为了符号的方便，因为W是一个矩阵，σ2是一个标量\sigma^2是一个标量σ2是一个标量）。
将xnx_nxn看作为缺失数据（missing data）,与观测数据一起构成完全数据(complete data),则有观测数据的对数似然函数为
log⁡(p(tn∣θ))=log⁡(∫p(tn,xn∣θ)dxn)=log⁡(∫p(tn,xn∣θ)q(xn)q(xn)dxn)≥∫log⁡(p(tn,xn)q(xn))q(xn)dxn(通过Jenson不等式求得)=∫(log⁡p(tn,xn∣θ))q(xn)dxn−∫(log⁡q(xn))q(xn)dxn=(1)−(2)(20)\begin{aligned} \log(p(t_n|\theta))&=\log(\int p(t_n,x_n|\theta)dx_n) \\ &=\log(\int \frac{p(t_n,x_n|\theta)}{q(x_n)}q(x_n)dx_n) \\ &\geq\int\log(\frac{p(t_n,x_n)}{q(x_n)})q(x_n)dx_n (通过Jenson不等式求得)\\ &=\int (\log p(t_n,x_n|\theta))q(x_n)dx_n-\int(\log q(x_n))q(x_n)dx_n\\ &=(1)-(2) \tag{20} \end{aligned} log(p(tn∣θ))=log(∫p(tn,xn∣θ)dxn)=log(∫q(xn)p(tn,xn∣θ)q(xn)dxn)≥∫log(q(xn)p(tn,xn))q(xn)dxn(通过Jenson不等式求得)=∫(logp(tn,xn∣θ))q(xn)dxn−∫(logq(xn))q(xn)dxn=(1)−(2)(20)
q(xn)q(x_n)q(xn)是关于xnx_nxn的函数，在此将其定义为在观测变量tnt_ntn给定时xnx_nxn的后验分布p(xn∣tn,θ)p(x_n|t_n,\theta)p(xn∣tn,θ)，由于时迭代算法，当初始值θ\thetaθ给定时，p(xn∣tn,θ)p(x_n|t_n,\theta)p(xn∣tn,θ)与参数θ\thetaθ无关。因此对于log⁡(p(tn∣θ))\log(p(t_n|\theta))log(p(tn∣θ))，与参数θ\thetaθ有关的仅有第（1）部分，则我们将（1）定义为Q函数：
Qn(θ∣θ(t))=∫(log⁡p(tn,xn∣θ))p(xn∣tn,θ(t)dxn=E[log⁡p(tn,xn∣θ)](21)Q_n(\theta|\theta^{(t)})=\int (\log p(t_n,x_n|\theta))p(x_n|t_n,\theta^{(t)}dx_n=E[\log p(t_n,x_n|\theta)] \tag{21}Qn(θ∣θ(t))=∫(logp(tn,xn∣θ))p(xn∣tn,θ(t)dxn=E[logp(tn,xn∣θ)](21)

2.1 计算后验分布

在此，我们先计算p(xn∣tn,θ)p(x_n|t_n,\theta)p(xn∣tn,θ)，也即潜变量得后验分布。
潜变量xnx_nxn先验分布假设为x∼N(0,Iq)x \sim N(0,I_q)x∼N(0,Iq),x的后验分布可有贝叶斯公式精确计算得知，为了简便，我们只根据分布函数的核进行判断：
p(xn∣tn)∝p(tn∣xn)p(xn)∝exp⁡{−(tn−Wxn−μ)′(tn−Wxn−μ)2τ}exp⁡{−xn′xn2}=exp⁡{−[xn′(W′W+τIq)xn−2xn′W′(tn−μ)]2τ}=exp⁡{−[xn′(W′W+τIq)xn]2τ+xn′W′(tn−μ)τ}=exp⁡{−12[xn′W′W+τIqτxn−2xn′W′(tn−μ)τ]}(22)\begin{aligned} p(x_n|t_n)&\propto p(t_n|x_n)p(x_n) \\ &\propto \exp\{-\frac{(t_n-Wx_n-\mu)'(t_n-Wx_n-\mu)}{2\tau}\}\exp\{-\frac{x_n'x_n}{2}\} \\ &=\exp\{-\frac{[x_n'(W'W+\tau I_q)x_n-2x_n'W'(t_n-\mu)]}{2\tau}\} \\ &=\exp\{-\frac{[x_n'(W'W+\tau I_q)x_n]}{2\tau}+\frac{x_n'W'(t_n-\mu)}{\tau}\} \\ &=\exp\{-\frac{1}{2}[x_n'\frac{W'W+\tau I_q}{\tau}x_n-2\frac{x_n'W'(t_n-\mu)}{\tau}]\} \tag{22} \end{aligned}p(xn∣tn)∝p(tn∣xn)p(xn)∝exp{−2τ(tn−Wxn−μ)′(tn−Wxn−μ)}exp{−2xn′xn}=exp{−2τ[xn′(W′W+τIq)xn−2xn′W′(tn−μ)]}=exp{−2τ[xn′(W′W+τIq)xn]+τxn′W′(tn−μ)}=exp{−21[xn′τW′W+τIqxn−2τxn′W′(tn−μ)]}(22)
假设xn∣tn∼N(α,Σ)x_n|t_n\sim N(\alpha,\Sigma)xn∣tn∼N(α,Σ),则其核函数可得：
p(xn∣tn)=exp⁡{−12(xn−α)′Σ−1(xn−α)}∝exp⁡{−12[xn′Σ−1xn−2xn′Σ−1α]}p(x_n|t_n)=\exp\{-\frac{1}{2}(x_n-\alpha)'\Sigma^{-1}(x_n-\alpha)\}\propto \exp\{-\frac{1}{2}[x_n'\Sigma^{-1}x_n-2x_n'\Sigma^{-1}\alpha]\}p(xn∣tn)=exp{−21(xn−α)′Σ−1(xn−α)}∝exp{−21[xn′Σ−1xn−2xn′Σ−1α]}

故根据上式得到Σ−1=W′W+τIqτ\Sigma^{-1}=\frac{W'W+\tau I_q}{\tau}Σ−1=τW′W+τIq,则Σ=τM−1,其中M=W′W+Iq\Sigma=\tau M^{-1},其中M=W'W+I_qΣ=τM−1,其中M=W′W+Iq,然而，W′(tn−μ)=Σ−1α,求得α=M−1W′(tn−μ)W'(t_n-\mu)=\Sigma^{-1}\alpha,求得\alpha=M^{-1}W'(t_n-\mu)W′(tn−μ)=Σ−1α,求得α=M−1W′(tn−μ),故给定观测数据tnt_ntn,潜变量xnx_nxn的分布为：
xn∣tn∼N(M−1W′(tn−μ),τM−1)(23)x_n|t_n \sim N(M^{-1}W'(t_n-\mu),\tau M^{-1}) \tag{23}xn∣tn∼N(M−1W′(tn−μ),τM−1)(23)

2.2 E Step

对Q函数求解就是EM算法的第一步（E step）,之所以称为E步，是因为我们可以将Q函数看作是完全数据的似然函数log⁡p(tn,xn∣θ)\log p(t_n,x_n|\theta)logp(tn,xn∣θ)关于xnx_nxn的分布p(xn∣tn,θ)p(x_n|t_n,\theta)p(xn∣tn,θ)的期望。
由于n=1,……,Nn=1,……,Nn=1,……,N，所以

Q=∑n=1NQn(θ∣θ(t))=∑n=1NE[log⁡p(tn,xn∣θ)]=∑n=1NE{−d2log⁡(2πτ)−q2log⁡(2π)−(tn−μ−Wxn)T(tn−μ−Wxn)2τ−xnTxn2}∝∑n=1NE{−d2log⁡τ−(tn−μ)T(tn−μ)2τ+xnTWT(tn−μ)τ−xnTWTWxn2τ−xnTxn2}(去掉与参数无关的部分)=−∑n=1N{d2log⁡τ+(tn−μ)T(tn−μ)2τ−E[xnTWT(tn−μ)]τ+E[xnTWTWxn]2τ+E[xnTxn]2}=−∑n=1N((1)+(2)−(3)+(4)+(5))(24)\begin{aligned} Q&=\sum_{n=1}^{N}Q_n(\theta|\theta^{(t)})=\sum_{n=1}^{N}E[\log p(t_n,x_n|\theta)] \\ &=\sum_{n=1}^{N}E\{-\frac{d}{2}\log(2\pi\tau)-\frac{q}{2}\log(2\pi)-\frac{(t_n-\mu-Wx_n)^T(t_n-\mu-Wx_n)}{2\tau}-\frac{x_n^{T}x_n}{2}\} \\ &\propto\sum_{n=1}^{N}E\{-\frac{d}{2}\log\tau-\frac{(t_n-\mu)^T(t_n-\mu)}{2\tau}+\frac{x_n^{T}W^T(t_n-\mu)}{\tau}-\frac{x_{n}^{T}W^{T}Wx_n}{2\tau}-\frac{x_n^{T}x_n}{2}\} (去掉与参数无关的部分)\\ &=-\sum_{n=1}^{N}\{\frac{d}{2}\log\tau+\frac{(t_n-\mu)^T(t_n-\mu)}{2\tau}-\frac{E[x_n^{T}W^T(t_n-\mu)]}{\tau}+\frac{E[x_{n}^{T}W^{T}Wx_n]}{2\tau}+\frac{E[x_n^{T}x_n]}{2}\} \\ &=-\sum_{n=1}^{N}((1)+(2)-(3)+(4)+(5)) \tag{24} \end{aligned}Q=n=1∑NQn(θ∣θ(t))=n=1∑NE[logp(tn,xn∣θ)]=n=1∑NE{−2dlog(2πτ)−2qlog(2π)−2τ(tn−μ−Wxn)T(tn−μ−Wxn)−2xnTxn}∝n=1∑NE{−2dlogτ−2τ(tn−μ)T(tn−μ)+τxnTWT(tn−μ)−2τxnTWTWxn−2xnTxn}(去掉与参数无关的部分)=−n=1∑N{2dlogτ+2τ(tn−μ)T(tn−μ)−τE[xnTWT(tn−μ)]+2τE[xnTWTWxn]+2E[xnTxn]}=−n=1∑N((1)+(2)−(3)+(4)+(5))(24)
因此根据以上计算结果，要计算Q函数，仅需要计算(3)(4)(5)式即可，也就是说，分别计算这三个式子关于分布p(xn∣tn,θ)p(x_n|t_n,\theta)p(xn∣tn,θ)的期望。

(3)式易得，直接得到结果：
E[xnTWT(tn−μ)]τ=[Exn]TWT(tn−μ)τ=⟨xn⟩TWT(tn−μ)τ\frac{E[x_n^{T}W^T(t_n-\mu)]}{\tau}=\frac{[Ex_n]^{T}W^T(t_n-\mu)}{\tau}=\frac{\langle x_n\rangle^{T}W^T(t_n-\mu)}{\tau}τE[xnTWT(tn−μ)]=τ[Exn]TWT(tn−μ)=τ⟨xn⟩TWT(tn−μ)
（4）和(5)的计算要使用二次型求期望公式：
E[X′AX]=μ′Aμ+tr(AΣ),其中μ=E[X],Σ=cov(X)E[X'AX]=\mu'A\mu+tr(A\Sigma),其中\mu=E[X],\Sigma=cov(X) E[X′AX]=μ′Aμ+tr(AΣ),其中μ=E[X],Σ=cov(X)
令A=W′W,E[xn]=αA=W'W,E[x_n]=\alphaA=W′W,E[xn]=α,则(4)式得：
E[xnTWTWxn]2τ=α′Aα+tr(AτM−1)2τ=tr(α′Aα)+tr(AτM−1)2τ=tr(Aαα′)+tr(AτM−1)2τ=tr(A(αα′+τM−1))2τ=tr(WTW⟨xnxnT⟩)2τ\begin{aligned} \frac{E[x_{n}^{T}W^{T}Wx_n]}{2\tau}&=\frac{\alpha'A\alpha+tr(A\tau M^{-1})}{2\tau} \\ &=\frac{tr(\alpha'A\alpha)+tr(A\tau M^{-1})}{2\tau} \\ &=\frac{tr(A\alpha\alpha')+tr(A\tau M^{-1})}{2\tau} \\ &=\frac{tr(A(\alpha\alpha'+\tau M^{-1}))}{2\tau} \\ &=\frac{tr(W^{T}W\langle x_nx_n^{T}\rangle)}{2\tau} \end{aligned}2τE[xnTWTWxn]=2τα′Aα+tr(AτM−1)=2τtr(α′Aα)+tr(AτM−1)=2τtr(Aαα′)+tr(AτM−1)=2τtr(A(αα′+τM−1))=2τtr(WTW⟨xnxnT⟩)
同理，(5)式可得，
E[xnTxn]2=α′α+tr(τM−1)2=tr(α′α)+tr(τM−1)2=tr(αα′)+tr(τM−1)2=tr(αα′+τM−1)2=tr(αα′+τM−1)2=tr(⟨xnxnT⟩)2\begin{aligned} \frac{E[x_n^{T}x_n]}{2}&=\frac{\alpha'\alpha+tr(\tau M^{-1})}{2} \\ &=\frac{tr(\alpha'\alpha)+tr(\tau M^{-1})}{2} \\ &=\frac{tr(\alpha\alpha')+tr(\tau M^{-1})}{2} \\ &=\frac{tr(\alpha\alpha'+\tau M^{-1})}{2} \\ &=\frac{tr(\alpha\alpha'+\tau M^{-1})}{2} \\ &=\frac{tr(\langle x_nx_n^{T}\rangle)}{2} \end{aligned}2E[xnTxn]=2α′α+tr(τM−1)=2tr(α′α)+tr(τM−1)=2tr(αα′)+tr(τM−1)=2tr(αα′+τM−1)=2tr(αα′+τM−1)=2tr(⟨xnxnT⟩)

故Q函数得到最终结果为：
Q=−∑n=1N{d2log⁡τ+(tn−μ)T(tn−μ)2τ−⟨xn⟩TWT(tn−μ)τ+tr(WTW⟨xnxnT⟩)2τ+tr(⟨xnxnT⟩)2}(25)Q=-\sum_{n=1}^{N}\{\frac{d}{2}\log\tau+\frac{(t_n-\mu)^T(t_n-\mu)}{2\tau}-\frac{\langle x_n\rangle^{T}W^T(t_n-\mu)}{\tau}+\frac{tr(W^{T}W\langle x_nx_n^{T}\rangle)}{2\tau}+\frac{tr(\langle x_nx_n^{T}\rangle)}{2}\} \tag{25} Q=−n=1∑N{2dlogτ+2τ(tn−μ)T(tn−μ)−τ⟨xn⟩TWT(tn−μ)+2τtr(WTW⟨xnxnT⟩)+2tr(⟨xnxnT⟩)}(25)

2.3 M Step

2.3.1 W的估计

对Q函数关于参数θ\thetaθ最大化，分别对θ\thetaθ的分量求导数即可。

d(Q(θ))w=−∑n=1N{−⟨xn⟩TdWT(tn−μ)τ+tr(d(WTW)⟨xnxnT⟩)2τ}=∑n=1N{tr[⟨xn⟩TdWT(tn−μ)]τ−tr((dWT∗W+WTdW)⟨xnxnT⟩)2τ}=∑n=1N{tr[(tn−μ)⟨xn⟩TdWT]τ−tr(W⟨xnxnT⟩)dWTτ}∝=∑n=1N{tr[((tn−μ)⟨xn⟩T−W⟨xnxnT⟩)dWT]}\begin{aligned} d(Q(\theta))_{w}&=-\sum_{n=1}^{N}\{-\frac{\langle x_n\rangle^{T}dW^T(t_n-\mu)}{\tau}+\frac{tr(d(W^{T}W)\langle x_nx_n^{T}\rangle)}{2\tau}\} \\ &=\sum_{n=1}^{N}\{\frac{tr[\langle x_n\rangle^{T}dW^T(t_n-\mu)]}{\tau}-\frac{tr((dW^{T}*W+W^{T}dW)\langle x_nx_n^{T}\rangle)}{2\tau}\} \\ &=\sum_{n=1}^{N}\{\frac{tr[(t_n-\mu)\langle x_n\rangle^{T}dW^T]}{\tau}-\frac{tr(W\langle x_nx_n^{T}\rangle)dW^{T}}{\tau}\} \\ &\propto=\sum_{n=1}^{N}\{tr[((t_n-\mu)\langle x_n\rangle^{T}-W\langle x_nx_n^{T}\rangle)dW^{T}]\} \end{aligned}d(Q(θ))w=−n=1∑N{−τ⟨xn⟩TdWT(tn−μ)+2τtr(d(WTW)⟨xnxnT⟩)}=n=1∑N{τtr[⟨xn⟩TdWT(tn−μ)]−2τtr((dWT∗W+WTdW)⟨xnxnT⟩)}=n=1∑N{τtr[(tn−μ)⟨xn⟩TdWT]−τtr(W⟨xnxnT⟩)dWT}∝=n=1∑N{tr[((tn−μ)⟨xn⟩T−W⟨xnxnT⟩)dWT]}

其中，
d(tr(W′W⟨xnxnT⟩))=tr(dW′∗W⟨xnxnT⟩)+tr(W′dW⟨xnxnT⟩)=tr(W⟨xnxnT⟩dW′)+tr(⟨xnxnT⟩dW′W)=tr(W⟨xnxnT⟩dW′)+tr(W⟨xnxnT⟩dW′)=2tr(W⟨xnxnT⟩dW′)\begin{aligned} d(tr(W'W\langle x_nx_n^{T}\rangle))&=tr(dW'*W\langle x_nx_n^{T}\rangle)+tr(W'dW\langle x_nx_n^{T}\rangle) \\ &=tr(W\langle x_nx_n^{T}\rangle dW')+tr(\langle x_nx_n^{T}\rangle dW'W) \\ &=tr(W\langle x_nx_n^{T}\rangle dW')+tr(W\langle x_nx_n^{T}\rangle dW') \\ &=2tr(W\langle x_nx_n^{T}\rangle dW') \\ \end{aligned}d(tr(W′W⟨xnxnT⟩))=tr(dW′∗W⟨xnxnT⟩)+tr(W′dW⟨xnxnT⟩)=tr(W⟨xnxnT⟩dW′)+tr(⟨xnxnT⟩dW′W)=tr(W⟨xnxnT⟩dW′)+tr(W⟨xnxnT⟩dW′)=2tr(W⟨xnxnT⟩dW′)
故
dQdW=∑n=1N{(tn−μ)⟨xn⟩T−W⟨xnxnT⟩}=0\frac{dQ}{dW}=\sum_{n=1}^{N}\{(t_n-\mu)\langle x_n\rangle^{T}-W\langle x_nx_n^{T}\rangle\}=0 dWdQ=n=1∑N{(tn−μ)⟨xn⟩T−W⟨xnxnT⟩}=0
则有，
W~=(∑n=1N(tn−μ)⟨xn⟩T)(∑n=1N⟨xnxnT⟩)−1(26)\tilde{W}=(\sum_{n=1}^{N}(t_n-\mu)\langle x_n\rangle^{T})(\sum_{n=1}^{N}\langle x_nx_n^{T}\rangle)^{-1} \tag{26} W~=(n=1∑N(tn−μ)⟨xn⟩T)(n=1∑N⟨xnxnT⟩)−1(26)

2.3.2 τ\tauτ的估计

∂Q(θ)∂τ=−∑n=1N{d2τ−12(tn−μ)′(tn−μ)τ2+⟨xn⟩TWT(tn−μ)τ2−12tr(WTW⟨xnxnT⟩)τ2}=0\begin{aligned} \frac{\partial Q(\theta)}{\partial \tau} &=-\sum_{n=1}^{N}\{\frac{d}{2\tau}-\frac{1}{2}\frac{(t_n-\mu)'(t_n-\mu)}{\tau^2}+\frac{\langle x_n\rangle^{T}W^T(t_n-\mu)}{\tau^2}-\frac{1}{2}\frac{tr(W^{T}W\langle x_nx_n^{T}\rangle)}{\tau^2}\}=0 \end{aligned}∂τ∂Q(θ)=−n=1∑N{2τd−21τ2(tn−μ)′(tn−μ)+τ2⟨xn⟩TWT(tn−μ)−21τ2tr(WTW⟨xnxnT⟩)}=0
则
τ~=1Nd∑n=1N{(tn−μ)′(tn−μ)−2⟨xn⟩TW~T(tn−μ)+tr(W~TW~⟨xnxnT⟩)}(27)\tilde{\tau}=\frac{1}{Nd}\sum_{n=1}^{N}\{(t_n-\mu)'(t_n-\mu)-2\langle x_n\rangle^{T}\tilde{W}^T(t_n-\mu)+tr(\tilde{W}^{T}\tilde{W}\langle x_nx_n^{T}\rangle)\} \tag{27}τ~=Nd1n=1∑N{(tn−μ)′(tn−μ)−2⟨xn⟩TW~T(tn−μ)+tr(W~TW~⟨xnxnT⟩)}(27)

由上文已知，
xn∣tn∼N(M−1W′(tn−μ),τM−1)(28)x_n|t_n \sim N(M^{-1}W'(t_n-\mu),\tau M^{-1}) \tag{28}xn∣tn∼N(M−1W′(tn−μ),τM−1)(28)
给定初始值W和τ\tauτ后，则xn∣tnx_n|t_nxn∣tn的分布已知，因此基于公式（26）（27）（28）则迭代过程如下：

给定初始值：W，τ=σ2W，\tau=\sigma^2W，τ=σ2,ttt代表第t次迭代

(1)⟨xn⟩(t+1)=M−1(t)WT(t)(tn−μ)\langle x_n\rangle^{(t+1)}=M^{-1(t)}W^{T(t)}(t_n-\mu)⟨xn⟩(t+1)=M−1(t)WT(t)(tn−μ)

(2)⟨xnxnT⟩)(t+1)=τM−1(t)+⟨xn⟩(t+1)⟨xn⟩(t+1)T\langle x_nx_n^{T}\rangle)^{(t+1)}=\tau M^{-1(t)}+\langle x_n\rangle^{(t+1)}\langle x_n\rangle^{(t+1)T}⟨xnxnT⟩)(t+1)=τM−1(t)+⟨xn⟩(t+1)⟨xn⟩(t+1)T

(3)W~t+1=(∑n=1N(tn−μ)⟨xn⟩T(t+1))(∑n=1N⟨xnxnT⟩(t+1))−1\tilde{W}^{t+1}=(\sum_{n=1}^{N}(t_n-\mu)\langle x_n\rangle^{T(t+1)})(\sum_{n=1}^{N}\langle x_nx_n^{T}\rangle^{(t+1)})^{-1}W~t+1=(∑n=1N(tn−μ)⟨xn⟩T(t+1))(∑n=1N⟨xnxnT⟩(t+1))−1

(4)τ~(t+1)=1Nd∑n=1N{(tn−μ)′(tn−μ)−2⟨xn⟩T(t+1)W~T(t+1)(tn−μ)+tr(W~T(t+1)W~(t+1)⟨xnxnT⟩)(t+1)}\tilde{\tau}^{(t+1)}=\frac{1}{Nd}\sum_{n=1}^{N}\{(t_n-\mu)'(t_n-\mu)-2\langle x_n\rangle^{T(t+1)}\tilde{W}^{T(t+1)}(t_n-\mu)+tr(\tilde{W}^{T(t+1)}\tilde{W}^{(t+1)}\langle x_nx_n^{T}\rangle)^{(t+1)}\}τ~(t+1)=Nd1∑n=1N{(tn−μ)′(tn−μ)−2⟨xn⟩T(t+1)W~T(t+1)(tn−μ)+tr(W~T(t+1)W~(t+1)⟨xnxnT⟩)(t+1)}

按照（1）-（4）的顺序进行迭代，直到算法收敛为止。

有时为了方便，将（28）式中的期望和方差代入（26）和（27），经过一系列计算得到：
W~t+1=SW(t)(τ(t)I+M−1(t)W(t)TSW(t))−1(29)\tilde{W}^{t+1}=SW^{(t)}(\tau^{(t)} I+M^{-1(t)}W^{(t)T}SW^{(t)})^{-1} \tag{29} W~t+1=SW(t)(τ(t)I+M−1(t)W(t)TSW(t))−1(29)
τ~(t+1)=1dtr{S−SW(t)M−1(t)W~t+1}(30)\tilde \tau^{(t+1)}=\frac{1}{d}tr\{S-SW^{(t)}M^{-1{(t)}}\tilde{W}^{t+1}\} \tag{30} τ~(t+1)=d1tr{S−SW(t)M−1(t)W~t+1}(30)

给定初始值后，直接对(29)和(30)进行迭代也可以。

3 模拟分析

3.1 极大似然估计R语言代码

PPCA.MLE= function(data,q){N = nrow(data) d = ncol(data) S = cov(data)*(N-1)/N S.e=eigen(S,symmetric = TRUE) S.values = S.e$valuesS.vectors = S.e$vectorsmu.hat=apply(data, 2, mean)sigma.hat = (1/(d-q))*sum(S.values[(q+1):d])U.q = S.vectors[,1:q] W.hat = U.q%*%sqrt(diag(S.values[1:q])-sigma.hat*diag(q))C = W.hat%*%t(W.hat)+sigma.hat*diag(d)L = -N*(d*log(2*pi)+log(det(C))+sum(diag(solve(C)%*%S)))/2return(list(mu.hat = mu.hat,sigma.hat = sigma.hat,W.hat = W.hat,L = L))
}

3.2 极大似然估计matlab代码

function [ mu,w,sigma,L ] = myppcaMle( data,q )
%   myppcaMle:          perform the ML estimation for ppca
%   data:               the input data
%   q:                  the dimension of latent space
%   mu:                 the mean of data
%   w:                  the factor loading matrix
%   sigma:              the noise variance
%   L:                  the log likelihood
[N,d] = size(data);
mu = mean(data,1);
T = data-repmat(mu,N,1); %Tn-mu   N x d matrix
S = T'*T/N; % sample covirance matrix
[V,D] = eig(S); % Eigenvalue decomposition
[D,ind] = sort(diag(D),'descend');
V = V(:,ind);
sigma = sum(D((q+1):d))/(d-q);
Uq = V(:,1:q);
lambda = D(1:q);
w = Uq*sqrt(diag(lambda)-sigma*eye(q));
C = w*w'+sigma*eye(d);
L = -N*(d*log(2*pi)+log(det(C))+trace(C\S))/2;
end

3.3 产生随机数模拟

产生500x50的多元正态分布随机数，调用myppcaMle.m文件求解

%matlab代码
N=500;d=50;
%C=[1 2 3;2 5 -1;3 4 0.5; 5 2 4; 8 4 3];
C=rand(d,5);
Sigma=C*C'+eye(d);
mu=zeros(1, d);
data=mvnrnd(mu,Sigma,N);
[ mu,w,sigma,L ] = myppcaMle( data,2 ) %q = 2

结果：
sigma =

0.9973

L =

-3.8156e+04
w =

1.5720    0.0438    0.0156   -0.0540   -0.0530
1.2005    0.1058    0.3620   -0.0587    0.2661
0.9459    0.2878   -0.1874   -0.1412    0.2027
1.0019    0.1230    0.0857    0.2362    0.0516
1.6361    0.2430    0.3402    0.1778   -0.2621
0.7754    0.1628   -0.1843   -0.3780   -0.0606
1.4058   -0.3866   -0.1356   -0.6775    0.0570
1.2503    0.7431    0.0505    0.0292   -0.0347
1.1661   -0.4553    0.5875    0.1110    0.1197
1.2999   -0.4238   -0.0797   -0.2630    0.3962
1.1997   -0.2367    0.3248   -0.4927    0.1052
1.7639   -0.3564    0.0873   -0.0554   -0.2347
1.8527   -0.2764   -0.1938    0.0309    0.3184
1.4023    0.2942    0.0876   -0.4853    0.3365
1.6301   -0.1035   -0.1867    0.2207    0.0515
0.7464   -0.1425    0.3904    0.0734    0.1356
1.1553    0.2479   -0.5237    0.2115    0.3468
1.4609   -0.1241    0.0029    0.3209    0.4468
1.2326    0.3062   -0.2927    0.0507   -0.0743
1.0559   -0.0836   -0.2309   -0.0832   -0.2286
1.4525   -0.0833    0.2430   -0.0228   -0.2699
1.1152    0.5510    0.4083    0.2175   -0.2681
1.2620   -0.1765   -0.2791    0.1263   -0.4167
1.1804    0.1549   -0.1592    0.1618    0.1028
0.8164    0.1464    0.1776   -0.3515   -0.1701
1.6320    0.4115   -0.2098   -0.1775   -0.2799
0.7181   -0.6356    0.3240   -0.1115    0.2148
0.5191    0.2463   -0.0182    0.0771    0.1920
0.6614    0.2289   -0.2178   -0.1585   -0.2960
0.8900   -0.0316   -0.0713    0.4507   -0.0870
1.4090   -0.1145   -0.1788    0.1291    0.3101
0.8474    0.7222    0.3297   -0.1545    0.2102
1.4167    0.1059   -0.7243    0.0525   -0.1818
1.1168   -0.4581   -0.3452    0.3257   -0.2130
1.3321   -0.2459    0.2864   -0.2336   -0.6533
1.1172   -0.0403    0.3609    0.4312    0.2921
0.9191    0.0462    0.3226    0.4194    0.2352
1.4706   -0.3615    0.3804    0.1931    0.0590
1.5295    0.1815    0.1910    0.1588   -0.2919
1.3108   -0.5129   -0.0599    0.1237   -0.0287
0.8889   -0.4877    0.2955   -0.3663   -0.0079
1.1872    0.2181   -0.3542    0.2937    0.2705
1.0097    0.1335   -0.3662   -0.0299    0.1271
0.7747    0.0458   -0.3526    0.0113   -0.1185
1.3889    0.7790    0.0525   -0.5751    0.2995
1.3609    0.5182    0.3476    0.0922   -0.1800
1.3972   -0.5621   -0.2327    0.2446   -0.0955
1.3566    0.3017    0.2629    0.4965   -0.2197
1.5782   -0.2135   -0.0658   -0.4633   -0.1950
1.6668   -0.2806   -0.3327   -0.1161    0.0951

3.4 EM算法R代码

ppcaEm = function(data,q){#data:  待输入的数据#q：待输入的潜变量维度data = as.matrix(data)N = nrow(data) #样本???d = ncol(data) #样本维度mu = colMeans(data) #样本均???I.q = diag(q)  #q维单位阵I.d = diag(d)  #d维单位阵w = matrix(rnorm(d*q,mean = 0,sd = 1),nrow = d,ncol = q) # W: 待输入的因子载荷矩阵初始???sigma = 1e-3 # 初始的噪声方???llh = c() #记录对数似然???muMat = t(matrix(rep(mu,N),ncol  = N)) #生成N*d的均值矩???S = t(data-muMat)%*%(data-muMat)/N#S = cov(data)*(N-1)/Nn = 100for (i in 1:n){C = w%*%t(w) + sigma*I.dM = t(w)%*%w + sigma*I.qinvC = solve(C)invM = solve(M)llh[i] = -(N/2)*(d*log(2*pi)+log(det(C))+sum(diag(invC%*%S))) #对数似然???postMeanX = invM%*%t(w)%*%t(data-muMat) #潜变量的后验均SumpostMeanX2 = N*sigma*invM + postMeanX%*%t(postMeanX) #w = t(data - muMat)%*%t(postMeanX)%*%solve(SumpostMeanX2)sigma = (1/(N*d))*(sum(diag((data-muMat)%*%t(data-muMat)))-2*sum(diag(t(postMeanX)%*%t(w)%*%t(data-muMat)))+sum(diag(SumpostMeanX2%*%t(w)%*%w)))}return(list(w = w,sigma = sigma,llh = llh))
}

3.5 EM算法matlab代码

function [ mu,w,sigma,L ] = myppca(data,q)
%myppca EM algorithm for ppca
%   data
%   q      factor numbers
%   mu     the estimated mean vector
%   w      weighted matrix
%   sigma  noise variance
%   L    likelihood
[N,d] = size(data);
mu = mean(data,1);
% q = 2
T = data-repmat(mu,N,1); %Tn-mu   N x d matrix
Iq = eye(q);% identity matrix q x q
Id = eye(d);
w = randn(d,q); % Initial value for w, d x q normal distribution random numbers
sigma = 1/randg; % Initial value for sigma
S = T'*T/N; % sample covirance matrix
niter = 500;
L = zeros(1,niter);
error = 1e-4;
for i = 2:niterM = w'*w+sigma*Iq;C = w*w'+ sigma*Id;detC = det(C);invC = inv(C);invM = inv(M);L(i) = -N*(d*log(2*pi)+log(detC)+trace(invC*S))/2;if abs(L(i)-L(i-1))<= error*abs(L(i-1));break;end% E stepEx = invM*w'*T'; %latent space mean of xExx = N*sigma*invM + Ex*Ex';% M stepw = T'*Ex'*inv(Exx);sigma = (dot(T(:),T(:))-2*trace(Ex'*w'*T')+trace(Exx*w'*w))/(N*d);
end
L = L(2:i);
end

3.6 模拟（matlab代码）

%matlab代码
N=500;d=50;
%C=[1 2 3;2 5 -1;3 4 0.5; 5 2 4; 8 4 3];
C=rand(d,5);
Sigma=C*C'+eye(d);
mu=zeros(1, d);
data=mvnrnd(mu,Sigma,N);
[ mu,w,sigma,L ] = myppca(data,5)

sigma =

0.9974

L =

1.0e+04 *

-5.1462 -4.0888 -3.8880 -3.8535 -3.8323 -3.8258 -3.8240 -3.8233 -3.8229 -3.8226

4 总结

通过自己编写的代码和示例代码进行对比，主要总结如下几点：

log⁡∣C∣\log |C|log∣C∣的计算
我的代码：直接计算；
示例代码：利用cholesky分解以及公式det(eA)=etrAdet(e^A)=e^{trA}det(eA)=etrA
设eA=B，则log⁡∣eA∣=log⁡∣B∣=log⁡trA=tr(log⁡B)设e^A=B，则\log|e^A|=\log |B|=\log trA=tr(\log B)设eA=B，则log∣eA∣=log∣B∣=logtrA=tr(logB)
对B矩阵进行cholesky分解，B=UTUB=U^T UB=UTU,故得：
log⁡∣B∣=tr(log⁡B)=tr(log⁡(UTU))=2tr(log⁡(U))\log |B|=tr(\log B)=tr(\log (U^T U))=2tr(\log(U))log∣B∣=tr(logB)=tr(log(UTU))=2tr(log(U))
C=WWT+σ2IC=W W^T+\sigma^2IC=WWT+σ2I和M=WTW+σ2IM=W^T W+\sigma^2IM=WTW+σ2I 的转换
为了计算方便，由于M的维度更低，所以利用M求解C
∣Id+ABT∣=∣Iq+ATB∣|I_d+AB^T|=|I_q+A^TB|∣Id+ABT∣=∣Iq+ATB∣
∣σ2Id+WWT∣=σ2d∣Id+WWTσ2∣|\sigma^2 I_d+WW^T|=\sigma^{2d}|I_d+\frac{WW^T}{\sigma^2}|∣σ2Id+WWT∣=σ2d∣Id+σ2WWT∣
∣σ2Iq+WTW∣=σ2q∣Id+WTWσ2∣|\sigma^2 I_q+W^TW|=\sigma^{2q}|I_d+\frac{W^TW}{\sigma^2}|∣σ2Iq+WTW∣=σ2q∣Id+σ2WTW∣

∣C∣=σ2(d−q)∣M∣|C|=\sigma^{2(d-q)}|M|∣C∣=σ2(d−q)∣M∣
C−1=σ−2(I−WM−1WT)C^{-1} = \sigma^{-2} (I-W M^{-1} W^T)C−1=σ−2(I−WM−1WT)

充分利用dot()函数
如矩阵范数∣∣A∣∣2||A||^2∣∣A∣∣2
示例代码：dot(A(:),(:))dot(A(:), (:))dot(A(:),(:))
我的代码：trace(A′A)trace(A'A)trace(A′A)

入门任务1-PPCA相关推荐

用Construct 2制作入门小游戏~
今天在软导课上了解到了Construct 2这个神器,本零基础菜鸟决定尝试做一个简单的小游戏(实际上是入门的教程啊= = 首先呢,肯定是到官网下载软件啊,点击我下载~ 等安装完毕后我便按照新手教程开始 ...
Docker入门六部曲——Swarm
原文链接:http://www.dubby.cn/detail.html?id=8738 准备工作安装Docker(版本最低1.13). 安装好Docker Compose,上一篇文章介绍过的. 安 ...
Docker入门六部曲——Stack
原文链接:http://www.dubby.cn/detail.html?id=8739 准备知识安装Docker(版本最低1.13). 阅读完Docker入门六部曲--Swarm,并且完成其中介绍 ...
Docker入门六部曲——服务
原文链接:http://www.dubby.cn/detail.html?id=8735 准备已经安装好Docker 1.13或者以上的版本. 安装好Docker Compose.如果你是用的是Do ...
【springboot】入门
简介: springBoot是spring团队为了整合spring全家桶中的系列框架做研究出来的一个轻量级框架.随着spring4.0推出而推出,springBoot可以説是J2SEE的一站式解决方案 ...
SpringBoot (一) ：入门篇 Hello World
什么是SpringBoot Spring Boot是由Pivotal团队提供的全新框架,其设计目的是用来简化新Spring应用的初始搭建以及开发过程.该框架使用了特定的方式来进行配置,从而使开发人员不 ...
入门指南目录页 -PaddlePaddle 飞桨入门指南 FAQ合集-深度学习问题
入门指南目录页 -PaddlePaddle 飞桨入门指南 FAQ合集 GT_Zhang关注 0.1012019.08.01 18:43:34字数 1,874阅读 795 Hi,欢迎各位来自Paddl ...
5 分钟入门 Google 最强NLP模型：BERT
BERT (Bidirectional Encoder Representations from Transformers) 10月11日,Google AI Language 发布了论文 BERT: ...
命名实体识别入门教程（必看）
关于开发自己的命名实体识别先期思路: 虽然网上有很多相关代码,但实际如何入门材料较少,故整理下: CRF:先期可以用人民日报语料库去做,步骤如下: https://blog.csdn.net/hude ...
Shiro第一个程序:官方快速入门程序Qucickstart详解教程
目录一.下载解压二.第一个Shiro程序 1. 导入依赖 2. 配置shiro配置文件 3. Quickstart.java 4. 启动测试三.shiro.ini分析四.Quickstart. ...

入门任务1-PPCA