EM算法 - 笔记

机器学习入门资料

---

1.符号理解：

$L(\theta|X)=P(X;\theta)$:表示$X$的概率密度函数，它是一个以$\theta$为参数的函数。分号左边是随机变量；右边是模型参数。
$P(X;\theta)=\sum_{Z}P(Z;\pi)P(X|Z;\theta)$:前者是Z的以$\pi$为参数的先验概率密度函数，后者是给定Z下的条件概率。

引 Jensen 不等式

如图所示：

1. 如若$f$是凸函数（$f''(x)\geq 0$）:
   有结论：$E(f(X))\geq f(E(X))$
2. 如若$f$是凸函数（$f''(x)\leq 0$）:
   有结论：$E(f(X))\leq f(E(X))$
3. 推论：当且仅当$x$为常数的概率为1时，取到等号
   或者$X=E(X)$

3.最大似然估计（MLE）

$$\theta^*=arg \max \limits_{\theta}\sum_XlogL(\theta|X)$$
当$Z$已知时，对$\theta$求导容易求解$\theta^*$。
$$\theta^*=arg \max \limits_{\theta}\sum_XlogL(\theta|X)\\ =arg \max \limits_{\theta}\sum_XlogP(X;\theta)\\=arg \max \limits_{\theta}\sum_Xlog\sum_ZP(Z;\pi)P(X|Z;\theta)$$
当$Z$未知时，不能通过对$\theta$求导获得。

4.EM algorithm

4.1E-Step

选择合适的关于$Z$的分布$Q(Z)$。

$$Max=\max \limits_{\theta}\sum_Xlog\sum_ZP(X,Z;\theta)\\=\max \limits_{\theta}\sum_Xlog\sum_ZQ(Z;\theta)\cdot\frac{P(X,Z;\theta)}{Q(Z;\theta)}\\=\max \limits_{\theta}\sum_XlogE_Q[\frac{P(X,Z;\theta)}{Q(Z;\theta)}]\\\geq\max \limits_{\theta}\sum_XE_Q[log\frac{P(X,Z;\theta)}{Q(Z;\theta)}]\\=\max \limits_{\theta}\sum_X\sum_ZQ(Z;\theta)log\frac{P(X,Z;\theta)}{Q(Z;\theta)}$$
由Jensen不等式结论可知：
$$\frac{P(X,Z;\theta)}{Q(Z;\theta)}=constant$$
等号才能取到。且$Q(Z;\theta)$是$Z$的某一分布，有$\sum_ZQ(Z;\theta)=1$这个性质：
$$Q(Z;\theta)=\frac{P(X,Z;\theta)}{c}=\frac{P(X,Z;\theta)}{c\cdot\sum_ZQ(Z;\theta)}\\=\frac{P(X,Z;\theta)}{\sum_Zc\cdot Q(Z;\theta)}=\frac{P(X,Z;\theta)}{\sum_ZP(X,Z;\theta)}\\=\frac{P(X,Z;\theta)}{P(X;\theta)}=P(Z|X;\theta)$$
由此可知，$Q$应取为给定$X$下，$Z$的后验概率分布。

4.2M-Step

寻找参数最大化期望似然

$$\theta^*=arg \max \limits_{\theta}\sum_X\sum_ZP(Z|X;\theta)log\frac{P(X,Z;\theta)}{P(Z|X;\theta)}$$

直接求导，依然很麻烦，不过可以用迭代来最大化
(1)由$\theta^{(j)}$，$Q^{(j)}=P(Z|X;\theta^{(j)})$
(2)

$$Q^{(j+1)}=arg\max\limits_{\theta}\sum_X\sum_ZQ^{(j)}log\frac{P(X,Z;\theta)}{Q^{(j)}}\\=arg\max\limits_{\theta}\sum_X\sum_Z(Q^{(j)}logP(X,Z;\theta)-Q^{(j)}logQ^{(j)})\\=arg\max\limits_{\theta}\sum_X\sum_ZQ^{(j)}logP(X,Z;\theta)$$
只需要最大化联合分布$P(X,Z;\theta)$，求出$(\theta^{(j+1)})$

5.应用于高斯混合聚类：

$x$概率密度函数为：

$$p(x)=\frac{1}{(2\pi)^{n/2}|\Sigma|^{\frac 12}}e^{-\frac 12(x-\mu)^T\sum^{-1}(x-\mu)}$$
其中，$\mu$是n维均值向量，$\sum$是$n\times n$的协方差矩阵。
高斯混合分布：
$$p_{M}(x)=\sum_{i=1}^k\alpha_ip(x|\mu_{i},\Sigma_i)$$
分布由$k$个混合成分组成，每个混合部分对应一个高斯分布。其中$\mu_i$与$\sigma_i$是第$i$个高斯混合成分的参数，而$\alpha_i\geq 0$为相应的“混合系数”，且$\sum_{i=1}^k\alpha_i=1$。
E-step
训练集：$D=\{x_1,x_2,\cdots,x_m\}$
随机变量（高斯混合成分）：$z_j\in\{1,2,\cdots,k\}$
$z_j$的先验概率$P(z_j=i)$：$\alpha_i(i=1,2,\cdots,k)$
因此，$z_j$的后验分布：
$${{p}_{M}}({{z}_{j}}=i|{{x}_{j}})=\frac{P({{z}_{j}}=i)\cdot {{p}_{M}}({{x}_{j}}|{{z}_{j}}=i|)}{{{p}_{M}}({{x}_{j}})} \\ =\frac{{{\alpha }_{i}}\cdot p({{x}_{j}}|{{\mu }_{i}},{{\Sigma }_{i}})}{\sum\limits_{l=1}^{k}{{{\alpha }_{j}}\cdot p({{x}_{j}}|{{\mu }_{i}},{{\Sigma }_{i}})}}$$
M-step
模型参数求解，给定样本集$D$,极大似然估计，最大化（对数）似然
$$LL(D)=\ln (\prod\limits_{j=1}^{m}{{{p}_{M}}({{x}_{j}})}) \\=\sum\limits_{j=1}^{m}{\ln (\sum\limits_{i=1}^{k}{{{\alpha }_{i}}\cdot p({{x}_{j}}|{{\mu }_{i}},{{\Sigma }_{i}})})}$$

令$\gamma_{ji}=p_M(z_j=i|x_j)$。
更新$\mu_i^*$:

$$\mu_i=\frac{\sum_{j=1}^m\gamma_{ji}x_j}{\sum_{j=1}^m\gamma_{ji}}$$
各混合成分的均值可通过样本加权平均来获得，样本的权重是每个样本属于该成分的后验概率。
更新$\Sigma_i^*$:
$$\Sigma_i=\frac{\sum_{j=1}^m\gamma_{ji}(x_j-\mu_i)(x_j-\mu_i)^T}{\sum_{j=1}^m\gamma_{ji}}$$
更新混合系数$\alpha_i$
因要满足$\alpha_i\geq 0$,$\sum_{i=1}^k\alpha_i=1$，考虑$LL(D)$拉格朗日形式:
$$LL(D)+\lambda(\sum_{i=1}^k\alpha_i-1)$$
$\lambda=-m$
$$\alpha_i=\frac{1}{m}\sum_{j=1}^m\gamma_{ji}$$
每个高斯成分的混合系数由样本属于改成份的平均后验概率确定。

分簇

$$\lambda_j=arg\max\limits_{i\in\{1,2,\cdots,k\}}\gamma_{ji}$$
簇划分是由原型对应后验概率确定。

MOG-EM Algorithm

伪代码：
输入：样本集 $D=\{x_1,x_2,\cdots,x_m\};$
$\qquad$高斯混合分布个数$k$.
过程：
1:初始化高斯混合分布的模型参数$\{(\alpha_i,\mu_i\Sigma_i)|1\leq i\leq k\}$
2:repeat
3:$\qquad$for $j=1,2,\cdots,m$ do
4:$\qquad$$\qquad$计算$x_j$由各混合成分生成的后验概率，即
$\qquad$$\qquad$$\gamma_{ji}=p_m(z_j=i|x_j)(1\leq i\leq k)$
5:$\qquad$end for
6:$\qquad$for $i=1,2,\cdots,k$ do
7:$\qquad$$\qquad$计算新均值向量：$\mu_i'=\frac{\sum_{j=1}^m\gamma_{ji}x_j}{\sum_{j=1}^m\gamma_{ji}};$
8:$\qquad$$\qquad$计算新协方差矩阵：$\Sigma_{i}'=\frac{\sum_{j=1}^m\gamma_{ji}(x_j-\mu_i')(x_j-\mu_i')^T}{\sum_{j=1}^m\gamma_{ji}};$
9:$\qquad$$\qquad$计算新混合系数：$\alpha_i'=\frac{\sum_{j=1}^m\gamma_{ji}}{m};$
10:$\qquad$end for
11:$\qquad$将模型参数$\{(\alpha_i,\mu_i\Sigma_i)|1\leq i\leq k\}$更新为$\{(\alpha_i’,\mu_i',\Sigma_i’)|1\leq i\leq k\}$
12:until满足停止条件
13:$C_i=\varnothing(1\leq i\leq k)$
14:for $j=1,2,\cdots,m$ do
15:$\qquad$确定$x_j$簇标记$\lambda_j;$
16:$\qquad$将$x_j$划入相应的簇$C_{\lambda_j}=C_{\lambda_j}\bigcup\{x_j\}$
17:end for
输出：簇划分$C=\{C_1,C_2\cdots,C_k\}$

%% Example of the Mixture-of-Gaussian EM algorithm
% 2018.3.28
% by WuQi 
%% Algorithm Start
clear;clc;
tic;                               % Timing start
%% Read data
Data = load('data.txt');
%% Initialization parameters
k = 3; sample = 4;                 % cluster's number is 3
[m,n] = size(Data);                % m = 80, n =2
alpha = [1/3 1/3 1/3];             % parameter alpha = 1/3 
mu = [Data(6,:); Data(22,:); Data(27,:)];           % prior distribution
Sigma(:,:,1)=[0.1,0.0;0.0,0.1];                     % parameter sigma = [0.1,0.0;0.0,0.1]
Sigma(:,:,2)=[0.1,0.0;0.0,0.1];                     
Sigma(:,:,3)=[0.1,0.0;0.0,0.1];
itera = 50;                                         % iteration is 50
count = 0;                                          
Save_mu = zeros(k,n,itera);                         % parameter mu 
Save_index = zeros(m,1,itera);                      % index of cluster
Sample_itera = [5 10 20 50];                        % sampling piont
while count <= itera
%% E-Step :Posterior Distribution
    for j = 1:m
        sum1=zeros(1,3);
        for i = 1:k
            sum1(j,i) = alpha(i) * Gf(Data(j,:),mu(i,:),Sigma(:,:,i));      
        end
        sum2=sum(sum1(j,:));
        gamma(j,:) = sum1(j,:) / sum2;              % gained 
    end
    [max_gamma,index] = max(gamma,[],2);            % classification ; large values are retained
%% M-Step :Update the Parameters
    sum3 = sum(gamma);
    for i = 1:k
        sum4 = 0; sum5 = zeros(2,2);
        for j = 1:m
            sum4 = sum4 + gamma(j,i) * Data(j,:);
        end
        mu1(i,:) = sum4 / sum3(i);                  % update the mu
        for j = 1:m
            Temp = Data(j,:)-mu1(i,:);
            Temp = Temp' * Temp;
            sum5 = sum5 + gamma(j,i) * Temp;
        end
        Sigma1(:,:,i) = sum5 / sum3(i);             % update the sigma
        alpha1(i,:) = sum3(i) / m;                  % update the alpha
    end
    mu = mu1;
    Sigma = Sigma1;
    alpha = alpha1;
    count = count + 1;
    Save_mu(:,:,count) = mu;                        % save the parameter mu
    Save_index(:,:,count) = index;                  % save the parameter index
end
%% Plot and Classification
shape = ['o' 's' '^'];                              % point's shape
color = ['p' 'g' 'k'];                              % point's color
padding = 'filed';
figure(1);
subplot(221);hold on;                               % 221
T_index = Save_index(:,:,Sample_itera(1));
for i = 1:k
    Temp = find(T_index == i);
    scatter(Data(Temp,1),Data(Temp,2),shape(i),color(i),padding);
    Temp = [];
end
scatter(Save_mu(:,1,5),Save_mu(:,2,5),'+','r');
xlabel('密度');ylabel('含糖率');
title('(a)5轮迭代后');
subplot(222);hold on;                               % 222
T_index = Save_index(:,:,Sample_itera(2));
for i = 1:k
    Temp = find(T_index == i);
    scatter(Data(Temp,1),Data(Temp,2),shape(i),color(i),padding);
    Temp = [];
end
scatter(Save_mu(:,1,10),Save_mu(:,2,10),'+','r');
xlabel('密度');ylabel('含糖率');
title('(b)10轮迭代后');
subplot(223);hold on;                               % 223
T_index = Save_index(:,:,Sample_itera(3));
for i = 1:k
    Temp = find(T_index == i);
    scatter(Data(Temp,1),Data(Temp,2),shape(i),color(i),padding);
    Temp = [];
end
scatter(Save_mu(:,1,20),Save_mu(:,2,20),'+','r');
xlabel('密度');ylabel('含糖率');
title('(c)20轮迭代后');
subplot(224);hold on;                               % 224
T_index = Save_index(:,:,Sample_itera(4));
for i = 1:k
    Temp = find(T_index == i);
    scatter(Data(Temp,1),Data(Temp,2),shape(i),color(i),padding);
    Temp = [];
end
scatter(Save_mu(:,1,50),Save_mu(:,2,50),'+','r');
xlabel('密度');ylabel('含糖率');
title('(d)50轮迭代后');
toc;                                % Timing end
---
%% Likehood Function
function f=Gf(x,u,s)
sum1= (-1/2)*(x-u)*(inv(s))*(x-u)';
sum2= 1/(2*pi*det(s)^(1/2));
f=sum2*exp(sum1);
end

效果图

高斯混合聚类($k=3$)在不同轮迭代后的聚类结果。其中样本簇$C_1,C_2$与$C_3$中的样本点分别用“圆形”，“方块”“三角形”表示，各高斯混合成分的均值向量用"+"表示

6.应用于朴素贝叶斯

Example（三硬币模型）
假设有3枚硬币，分别记作A,B,C。这些硬币正面出现的概率分别是$\pi,p$和$q$.
实验:

1. 先投掷硬币A，根据其结果选择硬币B或者硬币C。正面选择硬币B，反面选择硬币C。
2. 然后投掷被选出的硬币，对出现的结果记录。正面记作1；反面计算0

样本：$n=10$
1,1,0,1,0,0,1,0,1,1
只能观测到硬币投掷后的结果，不能观测其过程，如何估计三硬币正面出现的概率？即三硬币的模型参数。
解：
三硬币模型：

$$P(y|\theta)=\sum_{z}P(y,z|\theta)=\sum_{z}P(z|\theta)P(y|z,\theta)\\=\pi p^y(1-p)^{1-y}+(1-\pi)q^y(1-q)^{1-y}$$
这里，随机变量$y$是观测变量，表示一次试验观测的结果是1或0；随机变量$z$是隐变量，表示未观测到的投掷硬币A的结果。$\theta=(\pi,p,q)$是模型参数。

$$P(Y|\theta)=\sum_ZP(Z|\theta)P(Y|Z,\theta)$$
即：
$$P(Y|\theta)=\prod_{j=1}^n[\pi p^{y_j}(1-p)^{1-y_j}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]$$
求模型参数$\theta=(\pi,p,q)$的极大似然估计，即：
$$\widehat{\theta}=arg\max \limits_{\theta}logP(Y|\theta)$$
E-Step
选取初值，记作$\theta^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})$，

$$\mu^{(i+1)}=\frac{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}}{\pi^{(i)}(p^{(i)})^{y_j}(1-p^{(i)})^{1-y_j}+(1-\pi^{(i)})(q^{(i)})^{y_j}(1-q^{(i)})^{1-y_j}}$$

M-Step
计算模型参数的新估计值

$$\pi^{(i+1)}=\frac{1}{n}\sum_{j=1}^n\mu_j^{(i+1)}\\ p^{(i+1)}=\frac{\sum_{j=1}^{n}\mu_j^{(i+1)}y_j}{\sum_{j=1}^n\mu_j^{(i+1)}}\\ q^{(i+1)}=\frac{\sum_{j=1}^n(1-\mu_j^{(i+1)})y_j}{\sum_{j=1}^n(1-\mu_j^{(i+1)})}$$

NB-EM Algorithm

伪代码：
输入：样本集 $D=\{x_1,x_2,\cdots,x_m\};$
$\qquad$硬币A的状态为$k$.
过程：
1:初始化模型参数$\{\pi^{(0)},p^{(0)},q^{(0)}\}$
2:repeat
3:$\qquad$for $j=1,2,\cdots,m$ do
4:$\qquad$$\qquad$计算$x_j$由硬币B,C生成的后验概率，即
$\qquad$$\qquad$$\gamma_{ji}=p_B(z_j=i|x_j)$
5:$\qquad$end for
6:$\qquad$for $i=1,2,\cdots,k$ do
7:$\qquad$$\qquad$计算新的参数$p',q'$：$parameter_i'=\frac{\sum_{j=1}^m\gamma_{ji}x_j}{\sum_{j=1}^m\gamma_{ji}};$
8:$\qquad$$\qquad$计算新的参数$\pi'$：$\pi'=\frac{\sum_{j=1}^m\gamma_{ji}}{m};$
9:$\qquad$end for
10:$\qquad$将模型参数$\{\pi^{(i)},p^{(i)},q^{(i)}\}$更新为$\{\pi^{(i+1)},p^{(i+1)},q^{(i+1)}\}$
11:until满足停止条件
12:$C_i=\varnothing(1\leq i\leq k)$
13:for $j=1,2,\cdots,m$ do
14:$\qquad$确定$x_j$簇标记$\lambda_j;$
15:$\qquad$将$x_j$划入相应的簇$C_{\lambda_j}=C_{\lambda_j}\bigcup\{x_j\}$
16:end for
输出：簇划分$C=\{C_1,C_2\cdots,C_k\}$

%% Example of the Naive-Bayes EM algorithm
% 2018.3.29
% by WuQi 
%% Algorithm Start
clear;clc;
tic;                               % Timing start
%% Read data
Data = [1 1 0 1 0 0 1 0 1 1]';
%% Initialization parameters
k = 2;                             % cluster's number is 2
[m,n] = size(Data);                % m = 10, n = 1
p_pi = 0.4;                        % parameter pi
parameter = [0.6 0.7];             % parameter p,q
itera = 50;                        % iteration is 50
count = 0;          
Save_index = zeros(m,1,itera);                      % index of cluster
Sample_itera = [5 10 20 50];                        % sampling piont
while count <= itera
%% E-Step :Posterior Distribution  
    sum1 = zeros(m,k);
    for j = 1:m
        for i = 1:k
            sum1(j,i) = (p_pi^(2-i)) * ((1-p_pi)^(i-1)) * Bf(Data(j),parameter(i));
        end
        sum2=sum(sum1(j,:));              % gained 
        gamma(j,:) = sum1(j,:) / sum2;
    end
    [max_gamma,index] = max(gamma,[],2);            % classification ; large values are retained
%% M-Step :Update the Parameters
    sum3 = sum(gamma);
    p_pi1 = sum3(1) / m;                            % update the pi
    for i = 1:k
        sum4 = 0; 
        for j = 1:m
            sum4 = sum4 + gamma(j,i) * Data(j,:);
        end
        parameter1(i) = sum4 / sum3(i);             % update the parameter
    end
    p_pi = p_pi1;
    parameter = parameter1;
    count = count + 1;
    Save_index(:,:,count) = index;                  % save the parameter index
end
%% Plot and Classification
shape = ['o' 's' '^'];                              % point's shape
color = ['p' 'g' 'k'];                              % point's color
padding = 'filed';
figure(1);
subplot(221);hold on;                               % 221
T_index = Save_index(:,:,Sample_itera(1));
for i = 1:k
    Temp = find(T_index == i);
    [a,b] = size(Temp);
    scatter(1:a,Data(Temp),shape(i),color(i),padding);
    Temp = [];
end
xlabel('次');ylabel('面');
title('(a)5轮迭代后');
subplot(222);hold on;                               % 222
T_index = Save_index(:,:,Sample_itera(2));
for i = 1:k
    Temp = find(T_index == i);
    [a,b] = size(Temp);
    scatter(1:a,Data(Temp),shape(i),color(i),padding);
    Temp = [];
end
xlabel('次');ylabel('面');
title('(b)10轮迭代后');
subplot(223);hold on;                               % 223
T_index = Save_index(:,:,Sample_itera(3));
for i = 1:k
    Temp = find(T_index == i);
    [a,b] = size(Temp);
    scatter(1:a,Data(Temp),shape(i),color(i),padding);
    Temp = [];
end
xlabel('次');ylabel('面');
title('(c)20轮迭代后');
subplot(224);hold on;                               % 224
T_index = Save_index(:,:,Sample_itera(4));
for i = 1:k
    Temp = find(T_index == i);
    [a,b] = size(Temp);
    scatter(1:a,Data(Temp),shape(i),color(i),padding);
    Temp = [];
end
xlabel('次');ylabel('面');
title('(d)50轮迭代后');
p_pi
parameter
toc;                                % Timing end
---
%% Example of the Naive-Bayes EM algorithm
% Likehood Function
% 2018.3.29
% by WuQi 
function f=Bf(x,parameter)
sum1 = (parameter ^ x);
sum2 = ((1-parameter) ^ (1-x));
f=sum1 * sum2;
end

一般地：
$Y$：表示观测随机变量的数据，
$Z$：表示隐随机变量的数据，
$Y$和$Z$: 连在一起称为完全数据(complete - data)
假设：给定观测数据$Y$,其概率分布为$P(Y|\theta)$，其中$\theta$是需要估计得模型参数，对于不完全数据$Y$的似然函数是$P(Y|\theta)$，对数似然函数$L(\theta)=logP(Y|\theta)$
假设：$Y$和$Z$的联合概率分布是$P(Y,Z|\theta)$，那么完整数据的对数似然函数是$logP(Y,Z|\theta)$。

最后补充：
1.EM算法对初始值很敏感。
2.停止迭代条件是：

$$\left\| {{\theta }^{(i+1)}}-{{\theta }^{(i)}} \right\|<{{\varepsilon }_{1}},\qquad or\qquad\|Q({{\theta }^{(i+1)}},{{\theta }^{(i)}})-Q({{\theta }^{(i)}},{{\theta }^{(i)}}) \|<{{\varepsilon }_{2}}$$
这里的$Q(\theta,\theta^{(i)})$:
$$Q(\theta ,{{\theta }^{(i)}})={{E}_{Z}}[\log P(Y,Z|\theta )|Y,{{\theta }^{(i)}}] \\=\sum\limits_{Z}{\log P(Y,Z|\theta )P(Z|Y,{{\theta }^{(i)}})} \\ $$

---