《神经网络与机器学习》第10章 EM算法

神经网络与机器学习

第10章 EM算法

§10.1 EM算法引入

EM算法是一种迭代算法，1977年Dempster等提出的，用于隐变量的概率模型参数的极大似然估计，或者极大后验概率估计。

EM算法两步：E步求期望，M步求极大，所以叫做期望最大算法，再高斯混合模型的学习中应用非常成功，后来又推广为广义EM算法。

隐变量：如果数据给定，模型变量都是可以观测到的，那么最大似然估计或者贝叶斯估计可直接进行，但是如果含有未知的观测过程，那么需要EM算法。

EM算法，实际上是含有 隐变量 ( h i d d e n   v a r i a b l e ) (hidden variable)(hidden variable)或说 潜在变量概率模型参数的最大似然估计法，构造EM算法的关键之一在于选择合适的隐变量。

An elegant and powerful method for finding maximum likelihood solutions for models with latent variables is called the expectation-maximization algorithm. —— From 《Pattern Recognition and Machine Learning》 § 9.2.2

例子：三硬币模型，3枚硬币分别记为A、B、C，单独抛下正面出现的概率分别是$\pi$，p，q，实验记录如下，抛A出现正面，则选择B，出现反面则选C，然后再抛硬币B或C，出现正面记录1，反面则记录0。独立进行了n=10次实验，结果为

1 1 0 1 0 0 1 0 1 1

对于观察者来说，只知道上述记录结果，不能观测到抛的过程，请估计A、B、C单独抛下正面出现的概率$\pi$，p，q？

https://zfoox.blog.csdn.net/article/details/103483863

https://www.cnblogs.com/pinard/p/6912636.html

https://blog.csdn.net/howard789/article/details/104903937/

这个模型参数估计没有解析解，设是观测变量结果，不是1就是0，设随机变量z为隐变量，是未观测到的A的结果，参数向量为$\theta = [\pi,p,q]$

1 建模：

考虑硬币正面（反面）出现的概率$p(y|\theta)$，y=1表示掷出硬币的结果为正面，y=0表示掷出硬币的结果为反面。那么，上面10次试验的观察结果就是

\[\left \{y_1,y_2,\cdots,y_{10}\right \}\]

以第一次投掷为例，y = 1表示我们观测到了"硬币出现正面"的结果，但是我们只能知道掷出了硬币的正面，并不知道到底是硬币B的正面（记为 z = 1），还是 硬币C的正面（记为 z = 0）。

因为选择投掷硬币B、还是硬币C来产生 "观测 y"的结果，是由投掷硬币A的结果z来确定，但这个中间过程是无法观测的。

 

如上面表格所示，如果单独考虑投掷硬币B ( z = 1 )，关于观测变量 y的概率是一个条件概率，即：

 \[Prob(y|z=1,\theta)=p^y(1-p)^{1-y},y\in\left \{0,1\right \}\]

如果单独考虑投掷硬币c ( z = 0)，关于观测变量 y的概率是一个条件概率，即：

  \[Prob(y|z=0,\theta)=q^y(1-q)^{1-y},y\in \left \{0,1\right \}\]

2. 引入隐藏变量

考虑第n次观测 $y_n$的产生过程：

如果投掷的是硬币B（必然伴随着事件 z = 1同时发生），关于的概率是联合概率

 \[Prob(y_n,z=1|\theta)=Prob(z=1|\theta)Prob(y_n|z=1,\theta)=\pi p^{y_n}(1-p)^{1-y_n}\]

如果投掷的是硬币C（必然伴随着事件 z = 0同时发生） ，关于的概率是联合概率

  \[Prob(y_n,z=0|\theta)=Prob(z=0|\theta)Prob(y_n|z=0,\theta)=(1-\pi） q^{y_n}(1-q)^{1-y_n}\]

上述过程实际上是通过引入隐藏变量 z来表示掷出硬币A的结果（该过程无法被观测）。引入隐藏变量 z之后， 虽然我们观测到的 N次独立试验的结果为

\[\left \{y_1,y_2,\cdots,y_{10}\right \},y_n\in \left \{0,1\right \}\]

 

实际上是${y_n,z_1}$或者${y_n,z_0}$。

最大似然解： 综上所述：每观测到一枚硬币的投掷结果（y = 1表示出现正面，y = 0表示出现反面），观测变量 y的概率表示为

 

\[Prob(y|\theta)=\sum_{z\in \left \{0,1\right \}}Prob(y|z,\theta)Prob(z|\theta)\\
\quad \quad \quad \quad = Prob(z=1|\theta)Prob(y|z=1,\theta)+Prob(z=0|\theta)Prob(y|z=0,\theta)\\
\quad \quad \quad \quad = \pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y}\]

设N次实验的观测数据用随机向量$\textbf{y}=\{y_1,y_2,\cdots,y_{10}\}$，$y_n \in \{0,1\}$,隐变量为$z=\{z_1,z_2,\cdots,z_{10}\}$，$z_n\in \{0,1\},$那么观测数据用随机向量y的似然函数为

 

\[P(\textrm{y}|\theta)=Prob(y_1,y_2,\cdots,y_{10}|\theta)=\prod_{n=1}^NProb(y_n|\theta)\\
\quad\quad\quad=\prod_{n=1}^N\left \{\sum_{z\in \left \{0,1\right \}}Prob(y_n,z|\theta)\right \}\\
\quad\quad\quad=\prod_{n=1}^N\left \{\pi p^{y_n}(1-p)^{1-y_n}+(1-\pi)q^{y_n}(1-q)^{1-y_n}\right \}\]

 

 

 

我们希望得到最大似然解

 \[\hat{\theta}=\underset{\theta}{argmax}P(\textrm{y}|\theta)\]

可惜没有解析解，就是求对数似然函数也不行。

3. 三硬币问题EM公式的推导过程

三硬币问题概率模型的参数为$\theta=[\pi,\textrm{p},\textrm{q}]$

(1)三硬币问题的 1-of-K 表示

隐藏变量$z_1=1$，表示A硬币正面，概率$Prob(z_1|\theta)=\pi_1$

隐藏变量$z_2=1$，表示A硬币反面，概率$Prob(z_2|\theta)=\pi_2=1-\pi_1$

则投掷硬币 A 的概率可以统一描述为：$Prob(z_k|\theta)=\pi_k$

合并都写成向量

 

\[Prob(z|\theta)=\prod_{k=1}Prob(z_k|\theta)=\pi_1^{z_1}\pi_2^{z_2},z=\begin{Bmatrix}
z_1\\
z_2
\end{Bmatrix}=\left \{\begin{bmatrix}
1\\
0
\end{bmatrix},\begin{bmatrix}
0\\
1
\end{bmatrix}
\right \}\]

 

(2)引入隐藏向量 表示观测值 y的概率

单独考虑硬币 B 和硬币 C 关于观测值$\textrm{y}_n \in \{0,1\}$的概率，为了数学上的描述方便，$\alpha_1=p$，$\alpha_2=q$，于是，硬币 B 和C 的概率可以统一描述为： 

\[Prob(\textrm{y}|z_k=1,\theta)=\pi_k\alpha_k^y(1-\alpha_k)^{1-y}\]

用隐藏向量z 表示则写为

\[Prob(\textrm{y}|z,\theta)\prod_{k=1}^{2}Prob(\textrm{y}|z_k=1,\theta)=\prod_{k=1}^{2}[\pi_k\alpha_k^y(1-\alpha_k)^{1-y}]^{z_k}\]

\[z=\begin{Bmatrix}
z_1\\
z_2
\end{Bmatrix}=\left \{\begin{bmatrix}
1\\
0
\end{bmatrix},\begin{bmatrix}
0\\
1
\end{bmatrix}\right \}\]

(3)非完整数据集的似然函数

对于样本集$Y=\{y_1,y_2,\cdots,y_N\}$和隐藏变量$Z=\{z_1,z_2,\cdots,z_N\}$，

\[z_n=\begin{Bmatrix}
z_{n1}\\
z_{n2}
\end{Bmatrix}=\left \{\begin{bmatrix}
1\\
0
\end{bmatrix},\begin{bmatrix}
0\\
1
\end{bmatrix}\right \}\]

\[P(Y|Z,\theta)=P(\left \{y_1,y_2,\cdots,y_N\right \}|\left \{z_1,z_2,\cdots,z_N\right \},\theta)\\
\quad\quad\quad=\prod_{n=1}^{N}P(y_n|Z,\theta)=\prod_{n=1}^{N}\left \{\prod_{k=1}^{2}[P(y_n|z_{nk}=1,\theta)]^{z_{nk}}\right \}\\
\quad\quad\quad=\prod_{n=1}^{N}\left \{ \prod_{k=1}^{2}[\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]^{z_{nk}}\right \}\]

条件概率

 \[P(Z|\theta)=P(\left \{z_1,z_2,\cdots,z_n\right \}|\theta)=\prod_{n=1}^{N}P(z_n|\theta)=\prod_{n=1}^{N}\left \{\prod_{k=1}^{2}\pi_k^{z_{nk}}\right \}\]

隐藏变量$Z=\{z_1,z_2,\cdots,z_N\}$无法观测，因此EM算法先对隐藏变量Z进行猜测，使其变成完备数据，然后再对参数向量$\theta$进行估计。使得观测数据从 imcomplete 形式 的Y变成 complete 形式 的(Y,Z)。如果Z=$\{z_1,z_2,\cdots,z_N\}$已知，那么完备数据(Y,Z)的概率分布

 

\[P(Y,Z|\theta)=P(Y|Z,\theta)P(Z|\theta)\\
\quad\quad=\prod_{n=1}^{N}\left \{ \prod_{k=1}^{2}[\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]^{z_{nk}}\right \}\prod_{n=1}^{N}\left \{\prod_{k=1}^{2}\pi_k^{z_{nk}}\right \}\\
\quad\quad=\prod_{n=1}^{N}\left \{ \prod_{k=1}^{2}\pi_k^{z_{nk}}[\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]^{z_{nk}}\right \}\]

 

对数似然函数

 \[lnP(Y,Z|\theta)=\sum_{n=1}^{N}\sum_{k=1}^{2}z_{nk}\left \{ln\pi_k+ln[\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]\right \}\]

我们不清楚每一次投掷用的是硬币B还是硬币C，也就是不知道 $z_n$的取值到底是下面何种向量

\[\left \{\begin{bmatrix}
1\\
0
\end{bmatrix},\begin{bmatrix}
0\\
1
\end{bmatrix}\right \}\]

那么对数似然函数的值$lnP(Y,Z|\theta)$无法确定。

(4) 计算对数似然函数$lnP(Y,Z|\theta)$ 的期望

由于隐藏向量的值无法观测，导致无法采用最大似然估计计算参数$\theta=[\pi,p,q]$的值。为了解决这个问题，E M算法对所有的隐藏变量进行了合理猜测，实际上就是计算其期望E[Z]。

a）当假定一个初始值$\theta^i=[\pi^i,p^i,q^i]$，就可以计算猜测值给定情况下的对数似然函数$lnP(Y,Z|\theta)$。

b）那么什么样子的猜测值是比较可靠的？则通过对于似然函数$lnP(Y,Z|\theta)$的期望进行评价，就是完全数据的对数似然函数似然函数$lnP(Y,Z|\theta)$关于给定观测数据Y和当前参数$\theta^i$下对于隐藏数据变量条件分布概率P$(Z|Y,\theta)$的期望，此期望函数为

\[Q(\theta,\theta^i)=\mathbb{E}_{P(Z|Y,\theta)}[lnP(Y,Z|\theta)]\\ \quad\quad=\sum_{n=1}^{N}\sum_{k=1}^{2}\mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]\{ln\pi_k+ln[\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]\}\\ \quad\quad \mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]=\sum_{z_{nk}\in \{0,1\}}z_{nk}P(z_n|y_n,\theta^i)\]

期望需要用到后验概率

\[P(\textbf{z}_n|y_n,\theta)=\frac{P(y_n,\textbf{z}_n|\theta)}{P(y_n|\theta)}=\frac{P(y_n|\textbf{z}_n,\theta)P(\textbf{z}_n|\theta)}{\sum_{z_n}P(y_n|\textbf{z}_n,\theta)P(\textbf{z}_n|\theta)}\\
\quad\quad\quad\quad=\frac{P(y_n|z_{nk}=1,\theta)P(z_{nk=1}|\theta)}{\sum_{k=1}^2P(y_n|z_{nk}=1,\theta)P(z_{nk=1}|\theta)}\\
\quad\quad\quad\quad=\frac{P(y_n|z_{nk}=1,\theta)P(z_{nk=1}|\theta)}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\]

那么期望

\[\mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]=\sum_{z_{nk}\in\{0,1\}}z_{nk}P(\mathbf{z}_n|\mathbf{y}_n,\mathbf{\theta}^i)\\
\quad\quad\quad\quad=\frac{1*P(y_n|z_{nk}=1,\theta^i)P(z_{nk}=1|\theta^i)+0*P(y_n|z_{nk}=0,\theta^i)P(z_{nk}=0|\theta^i)}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\\
\quad\quad\quad\quad=\frac{P(y_n|z_{nk}=1,\theta^i)P(z_{nk}=1|\theta^i)}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\\
\quad\quad\quad\quad=\frac{\prod_{k=1}^2 [\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]^{z_{nk}}\prod_{k=1}^2\pi_k^{z_{nk}}}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\]

因此，上面就是在初始条件$\mathbf{\theta}^i=[\pi^i,p^i,q^i]$，对于未知的隐变量的期望$\mathbb{E}_{P(Z|Y,\theta)}[z_nk]$代替函数$G(\mathbf{\theta},\mathbf{\theta}^i)$中$z_{nk}$值，从而完成似然函数的计算。

(5) 总结EM算法的期望

 \[\mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]=\frac{\prod_{k=1}^2 [\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]^{z_{nk}}\prod_{k=1}^2\pi_k^{z_{nk}}}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\]

如果观测数据是B硬币产生的，那么,$\mathbf{z}_n=[0,1]^{\mathrm{T}}$ ，$z_{n1}=1$，$z_{n2}=0$,则

 

\[\mathbb{E}_{P(Z|Y,\theta)}[z_{n1}=1]=\frac{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\\
\quad\quad\quad\quad\quad=\frac{\pi^{(i)}(p^{(i)})^{y_n}(1-p^{(i)})^{1-y_n}}{\pi^{(i)}(p^{(i)})^{y_n}(1-p^{(i)})^{1-y_n}+(1-\pi^{(i)}(q^{(i)}))^{y_n}(1-q^{(i)})^{1-y_n}}\]

如果观测数据是C硬币产生的，那么$\mathbf{z}_n=[0,1]^{\mathrm{T}}$,$z_{n1}=0$ ，$z_{n2}=1$，则

 

\[\mathbb{E}_{P(Z|Y,\theta)}[z_{n2}=1]=\frac{\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}{\pi_1\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}+\pi_2\alpha_2^{y_n}(1-\alpha_2)^{1-y_n}}\\
\quad\quad\quad\quad\quad=\frac{\pi^{(i)}(q^{(i)})^{y_n}(1-q^{(i)})^{1-y_n}}{\pi^{(i)}(p^{(i)})^{y_n}(1-p^{(i)})^{1-y_n}+(1-\pi^{(i)}(q^{(i)}))^{y_n}(1-q^{(i)})^{1-y_n}}\]

 

简单记为 

\[\mu_n^{(i+1)}=\mathbb{E}_{P(Z|Y,\theta)}[z_{n1}=1]\\
1-\mu_n^{(i+1)}=\mathbb{E}_{P(Z|Y,\theta)}[z_{n2}=1]\]

 

(6) EM算法的似然最大化M步

a）在上面求得可靠猜测值$\mathbf{Z}$，然后求出初始化参数值$\mathbf{\theta}^i=[\pi^i,p^i,q^i]$的似然值$Q(\theta,\theta^i)$

b）此时将$\theta$当成未知数，对于函数$Q(\theta,\theta^i)$进行似然估计，得到新的参数值$\theta^{i+1}=[p^{i+1},q^{i+1}]$，这就是最大化步（M步）

 

\[Q(\theta,\theta^i)=\mathbb{E}_{P(Z|Y,\theta)}[lnP(Y,Z|\theta)]\\
\quad\quad=\sum_{n=1}^{N}\sum_{k=1}^{2}\mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]\{ln\pi_k+ln[\alpha_k^{y_n}(1-\alpha_k)^{1-y_n}]\}\\
\quad\quad=\sum_{n=1}^{N}\mu_n^{(i+1)}\{ln\pi_1+ln[\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}]\}+\sum_{n=1}^{N}(1-\mu_n^{(i+1)})\{ln\pi_1+ln[\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}]\}\]

\[\frac{\partial Q(\theta,\theta^i)}{\partial \alpha_1}=\sum_{n=1}^N\mu_n^{(i+1)}\frac{\partial\{ln\pi_1+ln[\alpha_1^{y_n}(1-\alpha_1)^{1-y_n}]\}}{\partial \alpha_1}=\sum_{n=1}^{N}\mu_n^{(i+1)}\frac{y_n-\alpha_1}{\alpha_1(1-\alpha_1)}=0\]

\[p^{(i+1)}=\alpha_1^{(i+1)}=\frac{\sum_{n=1}^{N}\mu_n^{(i+1)}y_n}{\sum_{n=1}^{N}\mu_n^{(i+1)}}\]

同理

 \[\frac{\partial Q(\theta,\theta^i)}{\partial \alpha_2}=\sum_{n=1}^{N}(1-\mu_n^{(i+1)})\frac{y_n-\alpha_2}{\alpha_2(1-\alpha_2)}=0\]

 \[q^{(i+1)}=\alpha_2^{(i+1)}=\frac{\sum_{n=1}^{N}(1-\mu_n^{(i+1)})y_n}{\sum_{n=1}^{N}(1-\mu_n^{(i+1)})}\]

注意到我们有个约束问题那就是

 \[\pi_1+\pi_2=1\]

所以我们采用拉格朗日乘子法

 \[\underset{\pi_k,\alpha_k}{max}Q(\theta,\theta^i)+\lambda(\sum_{k}\pi_k-1)\]

对于$\pi_k$求导

 \[\sum_{n=1}^{N}\frac{\mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]}{\pi_k}+\lambda=0\]

\[\sum_{n=1}^{N} \mathbb{E}_{P(Z|Y,\theta)}[z_{nk}]+\pi_k\lambda=0\]

两边同时求和

 \[\sum_{k=1}^{2}\big\{\sum_{n=1}^{N}\mathbb{E}_{P(Z|Y,\theta}[z_{nk}]+\pi_k\lambda \big\}=0\]

注意

 \[\mathbb{E}_{P(Z|Y,\theta}[z_{nk}]=\frac{P(y_n|z_{nk}=1,\theta)P(z_{nk}=1|\theta)}{\sum_{k=1}^{2}P(y_n|z_{nk=1},\theta)P(z_{nk}=1|\theta)}\]

\[\sum_{k=1}^{2}\big\{\sum_{n=1}^{N}\mathbb{E}_{P(Z|Y,\theta)}
[z_{nk}]\big\}=\sum_{n=1}^{N}\sum_{k=1}^2\frac{P(y_n|z_{nk}=1,\theta)P(z_{nk}=1|\theta)}{\sum_{k=1}^{2}P(y_n|z_{nk=1},\theta)P(z_{nk}=1|\theta)}=N\]

那么

 \[\lambda=N\]

\[\pi_1^{i+1}=\frac{1}{N}\sum_{n=1}^{N}\mu_n^{(i+1)}\]

因此

 \[\pi_2^{i+1}=1-\pi_1^{i+1}=1-\frac{1}{N}\sum_{n=1}^{N}\mu_n^{(i+1)}=\frac{1}{N}\sum_{n=1}^{N}(1-\mu_n^{(i+1)})\]

上面是对于三硬币模型的EM算法，对于一般的问题描述如下：

Y表示观测数据，又称为不完全观测数据，因为含有隐藏变量，其概率分布$P(Y|\theta)$；Z表示隐藏变量的数据，不可观测，其概率分布$P(Z|\theta)$；设二者的联合概率分布$P(Y,Z|\theta)=P(Y|Z,\theta)P(Z|\theta)$。

EM算法就是通过迭代求似然函数$L(\theta)=logP(Y|\theta)$的极大估计：分成期望和最大化两步

1)选择初始参数$\theta^{(0)}$

2)E步：记$\theta^{(i)}$为第i步得到的参数迭代值，第i+1步计算期望

\[Q(\theta,\theta^{(i)})=\mathbb{E}_{P(Z|Y,\theta^{(i)})}[logP(Y,Z|\theta)]=\sum_z[lnP(Y,Z|\theta)]P(Z|Y,\theta^{(i)})\]

$P(Z|Y,\theta^{(i)})$是给定当前观测数据Y和第i步得到的参数迭代值$\theta^{(i)}$下的隐藏变量Z的条件分布。

3) M步：求使得$Q(\theta,\theta^{(i)})$最大化的参数，得到第i+1步得到的参数估计值$\theta^{(i+1)}$

 \[\theta^{(i+1)}=\underset{\theta}{argmax}Q(\theta,\theta^{(i)})\]

4)重复2)-3)直至收敛。

§10.2 EM算法理论基础

EM算法为何能够从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。

先理解Jensen不等式，以丹麦技术大学数学家约翰·延森（Johan Jensen）命名。它给出积分的凸函数值和函数的积分值间的关系。

凸函数定义：函数，，对于任意和，那么如果

则函数在区间内是（下）凸的。从一元函数图像可知，就是两点之间函数值总是在弦的下方。如，。

Jensen不等式：给定凸函数和随机变量，那么

当不等式取等号时,X=E[X],表示X是个常向量。

证明：我们只证明离散分布情况，如果两点分布，那么就符合前面的凸函数定义，有$\lambda=p_1$，$1-\lambda=p_2$，那么

 \[f(p_1x_1+p_2x_2)\geq p_1f(x_1)+p_2f(x_2)\]

假设离散分布点k-1个成立，

\[p_i=\frac{p_i}{1-p_k}\]

 

那么归纳法

\[\sum_{i=1}^{k}p_if(x_i)=p_kf(x_k)+(1-p_k)\sum_{i=1}^{k-1}p_if(x_i)\\
\geq p_kf(x_k)+(1-p_k)f(\sum_{i=1}^{k-1}p_ix_i)\\
\geq f\Big(p_kf(x_k)+(1-p_k)\sum_{i=1}^{k-1}p_ix_i\Big)=f\Big(\sum_{i=1}^{k}p_ix_i\Big)\]

  

因此离散随机变量分布是成立的，可以推广到连续随机变量。

 

EM算法的理论推导

对于一个隐变量的概率参数估计模型，我们试图极大化观测不完全数据Y关于参数$\theta$的对数似然函数

 \[L(\theta)=logP(Y|\theta)=log[\sum_{z}P(Y,Z|\theta)]=log[\sum_{z}P(Y|Z,\theta)P(Z|\theta)]\]

这样的最大化没有办法计算，因为隐变量的概率分布不知道以及log函数中加和。

EM算法是利用第i步得到的参数迭代值$\theta^{(i)}$近似最大化似然函数$L(\theta)$，希望新的估计$\theta$使得

 \[L(\theta)>  L(\theta^{(i)})\]

并且逐步达到极大。考虑

 

\[L(\theta)-L(\theta^{(i)})=log[\sum_zP(Y|Z,\theta)P(Z|\theta)]-logP(Y|\theta^{(i)}\\
\quad\quad\quad\quad=log\Big[\sum_zP(Z|Y,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}\Big]-logP(Y|\theta^{(i)})\\
\quad\quad\quad\quad\geq \sum_zP(Z|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-logP(Y|\theta^{(i)})\\
\quad\quad\quad\quad=\sum_zP(Z|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-\sum_zP(Z|Y,\theta^{(i)})logP(Y|\theta^{(i)})\\
\quad\quad\quad\quad=\sum_zP(Z|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\]  

因此我们得到

 \[L(\theta)\geq L(\theta^{(i)})+\sum_zP(Z|Z,\theta^{(i)})log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}=B(\theta,\theta^{(i)})\]

$B(\theta,\theta^{(i)})$为似然函数的一个下界，且$L(\theta^{(i)})=B(\theta,\theta^{(i)})$。因此任何使得$B(\theta,\theta^{(i)})$增大的参数估计都能提高似然函数$L(\theta)$，自然地我们选择极大化

 

\[\theta^{(i+1)}=\underset{\theta}{argmax}B(\theta,\theta^{(i)})\\
\quad\quad=\underset{\theta}{argmax}\Big\{L(\theta^{(i)})+\sum_zP(Z|Z,\theta^{(i)})
log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\Big\}\\
\quad\quad=\underset{\theta}{argmax}\Big\{\sum_zP(Z|Z,\theta^{(i)})logP(Y|Z,\theta)P(Z|\theta)\Big\}(L(\theta^{(i)})\quad is \quad subtracted)\\
\quad\quad=\underset{\theta}{argmax}\Big\{\sum_zP(Z|Z,\theta^{(i)})logP(Y,Z|\theta)\Big\}\\
\quad\quad=\underset{\theta}{argmax}Q(\theta,\theta^{(i)})\] 

EM算法就是不断极大化其下界逼近最大的似然函数。

 

如图，上面的曲线表示$L(\theta)$，下面的是下界曲线$B(\theta,\theta^{(i)})$，二者在当前参数$\theta^{(i)}$下相等，下一步寻找$\theta^{(i+1)}$使得$B(\theta,\theta^{(i)})$极大化，也是使得$\underset{\theta}{max}Q(\theta,\theta^{(i)})$，由于$L(\theta)\geq B(\theta,\theta^{(i)})$，所以也使得似然函数进一步增大，这一过程不断地进行下去，但是不能保证得到全局最大值。

补充定理1：为何构造的$log\Big[\sum_zP(Z|Y,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}\Big]$使用了隐变量的后验概率？这里可以设任意一个关于隐变量的函数$q(Z|\theta)$，那么

 

\[L(\theta)-L(\theta^{(i)})=log\Big[\sum_zP(Y,Z|\theta)\Big]-logP(Y|\theta^{(i)})\\
\quad\quad\quad=log\Big[\sum_zq(Z|\theta)\frac{P(Y,Z|\theta)}{q(Z|\theta)}\Big]-logP(Y|\theta^{(i)})\\
\quad\quad\quad\geq \sum_zq(Z|\theta)log\frac{P(Y,Z|\theta)}{q(Z|\theta)}-logP(Y|\theta^{(i)})\]

 

为使得$B(\theta,\theta^{(i)})$增大到与 $L(\theta)$相等，那么依据Jensen不等式等号条件

 \[\frac{P(Y,Z|\theta)}{q(Z|\theta)}=C\Rightarrow P(Y,Z|\theta)=Cq(Z|\theta)\Rightarrow \sum_zP(Y,Z|\theta)=C\sum_zq(Z|\theta)=C\]

再代回去得到

 \[\frac{P(Y,Z|\theta)}{\sum_zP(Y,Z|\theta)}=\frac{P(Y,Z|\theta)}{P(Y|\theta)}=q(Z|\theta)=P(Z|Y,\theta)\]

就是后验概率，此时

 \[L(\theta)=B(\theta,\theta^{(i)})\]

也就是说完整数据对数似然在隐变量分布下的期望是观测数据对数似然的下界，所以最大化观测数据对数似然转化为最大化这个下界。

补充定理2：收敛性。$\theta^{(i)}$是由EM算法得到的参数估计序列，那么似然函数序列$P(Y|\theta^{(i)})$是单调递增的，即

 \[P(Y|\theta^{(i+1)})\geq P(Y|\theta^{(i)})\]

证明：由于

 

\[P(Y|\theta)=\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}\\
logP(Y|\theta)=logP(Y,Z|\theta)-logP(Z|Y,\theta)\\
Q(\theta,\theta^{(i)})=\mathbb{E}_{P(Z|Y,\theta^{(i)})}[logP(Y,Z|\theta)]=\sum_z[logP(Y,Z|\theta)]P(Z|Y,\theta^{(i)})\]

设

 \[H(\theta,\theta^{(i)})=\mathbb{E}_{P(Z|Y,\theta^{(i)})}[logP(Z|Y,\theta)]=\sum_z[logP(Z|Y,\theta)]P(Z|Y,\theta^{(i)})\]

则似然函数

 \[logP(Y|\theta)=Q(\theta,\theta^{(i)})-H(\theta,\theta^{(i)})\]

那么

\[logP(Y|\theta^{(i+1)})-logP(Y|\theta^{(i)})\\=Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})-[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)})]\]

由于$\theta^{(i+1)}=\underset{\theta}{argmax}Q(\theta,\theta^{(i)})$，那么$Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})\geq 0$。而

 

\[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)})\\
=\sum_z[logP(Z|Y,\theta^{(i+1)})]P(Z|Y,\theta^{(i)})-\sum_z[logP(Z|Y,\theta^{(i)})]P(Z|Y,\theta^{(i)})\\
=\sum_zlog\Big[\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}\Big]P(Z|Y,\theta^{(i)})\\
\leq log\sum_z\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}P(Z|Y,\theta^{(i)})=log\sum_zP(Z|Y,\theta^{(i+1)})=log1=0\]

 

  

因此

 \[logP(Y|\theta^{(i+1)})-logP(Y|\theta^{(i)})\geq 0\]

补充定理3：如果观测数据分布$(Y|\theta)$存在上界，在似然函数$L(\theta^{(i)})=logP(Y|\theta^{(i)})$收敛到某一值。EM算法得到参数序列$\theta^{(i)}$收敛值$\theta^{(*)}$是对应着似然函数的局部极大值$L(\theta^{(*)})$。

 

三硬币模型Python程序

import numpy as np

 

y = np.array([1,1,0,1,0,0,1,0,1,1])

N = len(y)

pi_n = 0.4

p_n = 0.6

q_n = 0.7

flag = 1

iter = 0

while flag:

pi = pi_n

p = p_n

q = q_n

mu = np.zeros(N)

i = 0

for n in y:

t1 = pi*np.power(p,n)*np.power(1-p,1-n)

t2 = (1-pi)*np.power(q,n)*np.power(1-q,1-n)

mu[i] = t1/(t1+t2)

i = i + 1

 

pi_n = np.sum(mu)/N

p_n = np.sum(y*mu)/np.sum(mu)

q_n = np.sum(y*(1-mu))/np.sum(1-mu)

print(('%1.4f %5.4f %5.4f') % (pi_n,p_n,q_n))

iter = iter + 1

if iter==2:

flag = 0

Matlab版本程序如下

Y=[1 1 0 1 0 0 1 0 1 1];%观测数据

N=length(Y);

pi_n = 0.4;%硬币A正面的初始概率

p_n = 0.6;%硬币B正面的初始概率

q_n = 0.7;%硬币C正面的初始概率

flag = 1;%迭代判定

cond=1e-4;

while flag>cond

pi = pi_n

p = p_n

q = q_n

mu = zeros(1,N);

%E步李航统计学习公式9.5

t1 = pi*p.^Y.*(1-p).^(1-Y);

t2 = (1-pi)*q.^Y.*(1-q).^(1-Y);

mu= t1./(t1+t2);

%M步李航统计学习公式9.6-9.8

pi_n = sum(mu)/N

p_n = sum(Y.*mu)/sum(mu)

q_n = sum(Y.*(1-mu))/sum(1-mu)

flag=abs(p-p_n);

end

 

 

 

§10.3 EM算法与混合高斯

Kevin P. Murphy, Machine Learning: A probabilistic Perspective, MIT Press，London，2012.

使用最广泛的混合模型就是混合高斯模型（mixture of Gaussian MOG），也叫做高斯混合模型（Gaussian mixture model GMM）。在这个模型中，混合分布中的每个基分布都是多元高斯分布，其中均值为$\mu_k$，协方差矩阵为$\sum$。因此模型具有如下的形式：

 \[P(x_i|\boldsymbol{\theta})=\sum_{k=1}^K\alpha_k\phi(x_i|\boldsymbol{\theta}_k)\]

系数满足$\sum_{k=1}^K\alpha_k=1$，$\alpha_k \geq 0$,参数向量$\boldsymbol{\theta}_k=[\mu_k,\sum_k]$，$\phi$是高斯分布，也可以任意分布模型替代高斯分布。下图给出了3个高斯混合的分布模型示意图。

高斯混合模型数据分布图

EM算法：假设观测数据

\[\mathbf{y}=\{y_1,y_2,\cdots,y_N\}\]

 

高斯混合模型简写为

 \[P(\mathbf{y}|\boldsymbol{\theta})=\sum_{k=1}^K\alpha_k\phi(\mathbf{y}|\boldsymbol{\theta}_k)\]

参数包含了全部$\theta_k=\alpha_1,\cdots,\alpha_K,\mu_1,\cdots,\mu_K,\sigma_1,\cdots,\sigma_K$.

1) 隐变量与对数似然函数

假设观测数据产生过程$\mathbf{y}=\{y_1,y_2,\cdots,y_N\}$，$y_j$是搜先依据概率$\alpha_k$选择的第k个$\phi(\mathbf{y}|\boldsymbol{\theta}_k)$高斯分布生成的，虽然生成数据$y_j$已知，但是反映$y_j$来自于第k个分模型的数据知识未知，以隐变量$\gamma_{jk}$

 

\[\gamma_{jk}=\left\{\begin{matrix}
1，The\quad j-th\quad data\quad comes\quad from\quad the\quad k-th\quad submodel\\
0
\end{matrix}\right.\]

表示，j表示数据$y_j$下标，k表示第k个模型来源。那么完全数据应该表示为

 \[y_j,\gamma_{j1},\cdots,\gamma_{jk},j=1,2,\cdots,N\]

那么完全数据的似然函数

 

\[P(\mathbf{y},\boldsymbol{\gamma}|\boldsymbol{\theta})=\prod_{j=1}^{N}P(y_j,\gamma_{j1},\cdots,\gamma_{jk}|\boldsymbol{\theta})\\
=\prod_{j=1}^{N}\prod_{k=1}^{K}[\alpha_k\phi(y_j|\boldsymbol{\theta}_k)]^{\gamma_{jk}}\\
=\prod_{k=1}^{K}[\alpha_k]^{n_k}\prod_{j=1}^{N}[\phi(y_j|\boldsymbol{\theta}_k)]^{\gamma_{jk}},n_k=\sum_{j=1}^{N}\gamma_{jk},N=\sum_{k=1}^{K}n_k=\sum_{k=1}^{K}\sum_{j=1}^{N}\gamma_{jk}\\
=\prod_{k=1}^{K}[\alpha_k]^{n_k}\prod_{j=1}^{N}\Big[\frac{1}{\sqrt{2\pi}\sigma_k}exp\Big(\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\Big)\Big]^{\gamma_{jk}}\]

 

对数似然函数

\[logP(\mathbf{y},\boldsymbol{\gamma}|\boldsymbol{\theta})=\sum_{k=1}^K\Big\{{n_klog\alpha_k+\sum_{j=1}^N}\gamma_{jk}\Big[log\frac{1}{\sqrt{2\pi}}-log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\Big]\Big\}\]

2) E步

确定Q函数

 

\[Q(\theta,\theta^{(i)})=E[logP(\mathbf{y},\boldsymbol{\gamma}|\boldsymbol{\theta})|\mathbf{y},\theta^{(i)}]\\
=\sum_{k=1}^{K}\Big\{\sum_{j=1}^{N}E[\gamma_{jk}|\mathbf{y},\theta^{(i)}]log\alpha_k+\sum_{j=1}^{N}E[\gamma_{jk}|\mathbf{y},\theta^{(i)}]\Big[log\frac{1}{\sqrt{2\pi}}-log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\Big]\Big\}\]

 

那么注意就是计算

 \[E[\gamma_{jk}|\mathbf{y},\theta]=\widehat{\gamma}_{jk}=P(\gamma_{jk}=1|\mathbf{y},\theta)\]

利用后验概率

 

\[=\frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^{K}P(\gamma_{jk}=1,y_j|\theta)}\\
=\frac{P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^{K}P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}\\
=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^{K}\alpha_k\phi(y_j|\theta_k)}\]

 

$\widehat{\gamma}_{jk}$就是当前模型参数下观测数据$y_j$来自第k个分模型的概率，称为分模型对观测数据的响应度。那么将$\widehat{\gamma}_{jk}$，以及$n_k=\sum_{j=1}^{N}\widehat{\gamma}_{jk}$，代入则Q函数

\[Q(\theta,\theta^{(i)})==\sum_{k=1}^{K}\Big\{n_klog\alpha_k+\sum_{j=1}^{N}\widehat{\gamma}_{jk}\Big[log\frac{1}{\sqrt{2\pi}}-log\sigma_k-\frac{(y_j-\mu_k)^2}{2\sigma_k^2}\Big]\Big\}\]

3) M步

对于$Q(\theta,\theta^{(i)})$极大化，得到新的估计参数

 

\[\theta^{(i+1)}=\underset{\theta}{argmax}Q(\theta,\theta^{(i)})\\
\frac{\partial Q(\theta,\theta^{(i)})}{\partial \mu_k}=0\Rightarrow \widehat{\mu}_k=\frac{\sum_{j=1}^{N}\widehat{\gamma}_{jk}y_j}{\sum_{j=1}^{N}\widehat{\gamma}_{jk}}=\frac{\sum_{j=1}^{N}\widehat{\gamma}_{jk}y_j}{n_k}\\
\frac{\partial Q(\theta,\theta^{(i)})}{\partial \sigma_k^2}=0\Rightarrow \widehat{\sigma}_k^2=\frac{\sum_{j=1}^{N}\widehat{\gamma}_{jk}(y_j-\mu_k)^2}{n_k}\\
\frac{\partial\{Q(\theta,\theta^{(i)})+\lambda(\sum_{k=1}^{K}\alpha_k-1)\}}{\partial \alpha_k}=0\Rightarrow \widehat{\alpha}_k=\frac{\sum_{j=1}^{N}\widehat{\gamma}_{jk}}{N}\]

 

 

不断重复，直至收敛。

mixGaussDemoFaithful

mixGaussCreate

mixGaussFit

http://blog.sina.com.cn/yzhmatlab

%% 该程序为混合正态分布参数估计算法

%% 产生拟合数据

k = 3; % 高斯分布的个数

mu = [0 50 100];

sigma = [10 10 10];

nk = [1/3 1/3 1/3];

x = [];

for i = 1:k

x = [x normrnd(mu(i),sigma(i),1,300)];

end

%% 参数初始化：kmeans算法估计

idx = kmeans(x',k);

for i = 1:k

N(i) = length(x(idx == i));

a(i) = N(i)/ length(x);

u(i) = mean(x(idx == i));

o(i) = std(x(idx == i));

end

 

%% 迭代程序

t = 1;

while t < 100

Es =0;

for i =1:k

Es = Es +a(i) * normpdf(x,u(i),o(i));

end

for i =1:k

E = a(i) *normpdf(x,u(i),o(i));

N(i) = sum(E./Es);

a(i) = N(i)/length(x);

u(i) = sum(E./Es.*x)/N(i);

o(i) = sqrt(sum(E.*(x-u(i)).^2./Es)/N(i));

end

t = t +1;

end

disp('理论值')

[nk,mu,sigma]

disp('计算值')

[a,u,o]

理论值

0.3333 0.3333 0.3333 0 50.0000 100.0000 10.0000 10.0000 10.0000

计算值

0.3348 0.3360 0.3293 0.2673 49.7734 100.0348 10.0414 10.2329 9.3430

 

高维的高斯混合模型

EM算法：假设观测数据

\[\mathbf{y}=\{y_1,y_2,\cdots,y_N\}\]

 

高斯混合模型简写为

 \[P(\mathbf{y}|\boldsymbol{\theta})=\sum_{k=1}^K\alpha_k\phi(\mathbf{y}|\boldsymbol{\theta}_k)\]

参数包含了全部$\theta_k=[\alpha_1,\cdots,\alpha_K,\mu_1,\cdots,\mu_K,\sum_1,\cdots,\sum_K]$.

\[\phi(\mathbf{y}|\boldsymbol{\theta}_k)=\frac{1}{(\sqrt{2\pi})^D|\boldsymbol{\sum}_k|^{1/2}}exp\Big(-\frac{1}{2}(\mathbf{y}-\boldsymbol{\mu}_k)^{\mathrm{T}}\boldsymbol{\sum}_k^{-1}(\mathbf{y}-\boldsymbol{\mu}_k)\Big)\]

E步

 \[\widehat{\gamma}_{jk}==\frac{\alpha_k\phi(\mathbf{y}_j|\boldsymbol{\theta}_k)}{\sum_{k=1}^K\alpha_k\phi(\mathbf{y}_j|\boldsymbol{\theta}_k)}\]

M步

\[\widehat{\alpha}_k=\frac{\sum_{j=1}^N\widehat{\gamma}_{jk}}{N}\]

\[\widehat{\mu}_k=\frac{\sum_{j=1}^N\widehat{\gamma}_{jk}\mathbf{y}_j}{\sum_{j=1}^N\widehat{\gamma}_{jk}}=\frac{\sum_{j=1}^N\widehat{\gamma}_{jk}\mathbf{y}_j}{n_k}\]

 \[\sum_k=\frac{\sum_{j=1}^N\widehat{\gamma}_{jk}(\mathbf{y}-\boldsymbol{\mu}_k)(\mathbf{y}-\boldsymbol{\mu}_k)^{\mathrm{T}}}{n_k}\]

%程序中Psi对应二维高斯函数

%Gamma 为隐变量的值，Gamma(i,j)代表第i个样本属于第j个模型的概率。

%Mu为期望，Sigma为协方差矩阵

%Pi为各模型的权值系数

%产生2个二维正态数据

MU1 = [1 2];

SIGMA1 = [1 0; 0 0.5];

MU2 = [-1 -1];

SIGMA2 = [1 0; 0 1];

%生成1000行2列(默认)均值为mu标准差为sigma的正态分布随机数

X = [mvnrnd(MU1, SIGMA1, 1000);mvnrnd(MU2, SIGMA2, 1000)];

%显示

scatter(X(:,1),X(:,2),10,'.');

%====================

K=2;

[N,D]=size(X);

Gamma=zeros(N,K);

Psi=zeros(N,K);

Mu=zeros(K,D);

LM=zeros(K,D);

Sigma =zeros(D, D, K);

Pi=zeros(1,K);

%选择随机的两个样本点作为期望迭代初值

Mu(1,:)=X(randi([1 N/2],1,1),:);

Mu(2,:)=X(randi([1 N/2],1,1),:);

%所有数据的协方差作为协方差初值

for k=1:K

Pi(k)=1/K;

Sigma(:, :, k)=cov(X);

end

LMu=Mu;

LSigma=Sigma;

LPi=Pi;

while true

%Estimation Step

for k = 1:K

Y = X - repmat(Mu(k,:),N,1);

Psi(:,k) = (2*pi)^(-D/2)*det(Sigma(:,:,k))^(-1/2)*diag(exp(-1/2*Y/(Sigma(:,:,k))*(Y'))); %Psi第一列代表第一个高斯分布对所有数据的取值

end

Gamma_SUM=zeros(D,D);

for j = 1:N

for k=1:K

Gamma(j,k) = Pi(1,k)*Psi(j,k)/sum(Psi(j,:)*Pi');

%Psi的第一行分别代表两个高斯分布对第一个数据的取值

end

end

%Maximization Step

for k = 1:K

%update Mu

Mu_SUM= zeros(1,D);

for j=1:N

Mu_SUM=Mu_SUM+Gamma(j,k)*X(j,:);

end

Mu(k,:)= Mu_SUM/sum(Gamma(:,k));

%update Sigma

Sigma_SUM= zeros(D,D);

for j = 1:N

Sigma_SUM = Sigma_SUM+ Gamma(j,k)*(X(j,:)-Mu(k,:))'*(X(j,:)-Mu(k,:));

end

Sigma(:,:,k)= Sigma_SUM/sum(Gamma(:,k));

%update Pi

Pi_SUM=0;

for j=1:N

Pi_SUM=Pi_SUM+Gamma(j,k);

end

Pi(1,k)=Pi_SUM/N;

end

 

R_Mu=sum(sum(abs(LMu- Mu)));

R_Sigma=sum(sum(abs(LSigma- Sigma)));

R_Pi=sum(sum(abs(LPi- Pi)));

R=R_Mu+R_Sigma+R_Pi;

if ( R<1e-10)

disp('期望');

disp(Mu);

disp('协方差矩阵');

disp(Sigma);

disp('权值系数');

disp(Pi);

obj=gmdistribution(Mu,Sigma);

figure,h = ezmesh(@(x,y) pdf(obj,[x,y]),[-8 6], [-8 6]);

break;

end

LMu=Mu;

LSigma=Sigma;

LPi=Pi;

end

 

>> highDGMMEM

期望

-1.0259 -0.9681

1.0232 2.0108

 

协方差矩阵

 

(:,:,1) =

 

1.0229 0.0105

0.0105 1.0224

 

 

(:,:,2) =

 

0.9791 -0.0618

-0.0618 0.4688

 

权值系数

0.5097 0.4903