EM算法的仿真（三硬币+混合高斯）python 程序

1. 求解抛硬币的参数

(1)首先给出样本数据及说明

• 现在有A,B,C三个硬币，正面出现的概率分别为 $\pi$,$p$,$q$。
• 抛硬币的步骤：首先抛硬币A，硬币A正面朝上接下来抛硬币B，否则抛硬币C。
• 现在得到的样本数据为: 1,1,0,1,0,0,1,0,1,1。

(2)建立参数估计模型

• 假定未知的模型参数 $\mathbf{\theta }=(\pi ,p,q)$，每次抛掷硬币A（隐变量）的数据结果为 $z$，硬币B、C（观测变量）正反面的结果为 $y$。
  $$\begin{array}{r}P(y|\mathbf{\theta })=\sum _{z}P(y,z|\mathbf{\theta })=\sum _{z}P(y|z,\mathbf{\theta })P(z|\mathbf{\theta })\\ =\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$
• 假定观测数据 $\mathbf{y}=(y_1,y_2,...,y_n{)}^T$，观测次数为 $n$，则最大似然估计的公式为
  $$\mathbf{\theta }=\underset{\mathbf{\theta }}{\mathrm{argmax}}\log P(\mathbf{y}|\mathbf{\theta })$$
  其中 $$P(\mathbf{y}|\mathbf{\theta })=\prod _{i=0}^n[\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}]$$

(3)采用EM的步骤

• 计算隐变量
  假定此时估计的 ${\mathbf{\theta }}^{(i)}=({\pi }^{(i)},p^{(i)},q^{(i)})$。

  观测数据来自B的概率为：

  $${\mu }^{(i+1)}=\frac{{\pi }^{(i)}(p^{(i)}{)}^{y_i}{\left(1-p^{(i)}\right)}^{1-y_i}}{{\pi }^{(i)}(p^{(i)}{)}^{y_i}{\left(1-p^{(i)}\right)}^{1-y_i}+(1-{\pi }^{(i)})(q^{(i)}{)}^{y_i}{\left(1-q^{(i)}\right)}^{1-y_i}}$$
  观测数据来自C的概率为：

  $$1-{\mu }^{(i+1)}$$

• 此时对观测数据来自 B 的概率取平均值就可以得到参数 $${\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)})}$$

• 这里面仿真的是对于不同的初值最终 EM 算法会收敛到不同的概率预测值

sample = [1,1,0,1,0,0,1,0,1,1]
times = 20 #设定迭代参数为 20次
mu = []
def coin_EM(pi,p,q):
    pi_copy = pi
    q_copy = q
    p_copy = p
    for i in range(times):
        sum_p = 0
        sum_q = 0
        for flap in sample:
            if flap == 0:
                mu.append(pi*(1-p)/(pi*(1-p)+(1-pi)*(1-q)))
            else:
                mu.append(pi*p/(pi*p+(1-pi)*q))
                sum_p = sum_p + pi*p/(pi*p+(1-pi)*q)
                sum_q = sum_q + 1 - pi*p/(pi*p+(1-pi)*q)
        pi = sum(mu)/len(sample)
        p = sum_p/sum(mu)
        q = sum_q/(len(sample) - sum(mu))
        mu.clear()
    print("给定的初值分别为：")
    print("pi = {:.4f} , p = {:.4f} , q = {:.4f}".format(pi_copy,p_copy,q_copy))
    print("迭代得到的终值分别为：")
    print("pi = {:.4f} , p = {:.4f} , q = {:.4f}".format(pi,p,q))
    
    
def MLE(pi,p,q):
    index_one = 0
    index_zero = 0
    for clap in sample:
        if clap == 1:
            index_one = index_one + 1
    index_zero = len(sample) - index_one
    proability = pow(pi * p + (1-pi) * q,index_one) \
    * pow(pi * (1-p) + (1-pi) * (1-q),index_zero)
    print("对这一组的最大似然结果为：{:.5f}".format(proability))
    
print("第一组：")
coin_EM(0.46,0.55,0.67)
MLE(0.46,0.55,0.67)
print("###############################################")
print("第二组：")
coin_EM(0.5,0.5,0.5)
MLE(0.5,0.5,0.5)
print("###############################################")
print("第三组：")
coin_EM(0.4,0.6,0.7)
MLE(0.4,0.6,0.7)

#结果
第一组：
给定的初值分别为：
pi = 0.4600 , p = 0.5500 , q = 0.6700
迭代得到的终值分别为：
pi = 0.4619 , p = 0.5346 , q = 0.6561
对这一组的最大似然结果为：0.00119
###############################################
第二组：
给定的初值分别为：
pi = 0.5000 , p = 0.5000 , q = 0.5000
迭代得到的终值分别为：
pi = 0.5000 , p = 0.6000 , q = 0.6000
对这一组的最大似然结果为：0.00098
###############################################
第三组：
给定的初值分别为：
pi = 0.4000 , p = 0.6000 , q = 0.7000
迭代得到的终值分别为：
pi = 0.4064 , p = 0.5368 , q = 0.6432
对这一组的最大似然结果为：0.00110

2. 求解混合高斯分布的参数

(1)首先生成并展示出混合高斯分布的数据

• 高斯混合模型的概率分布如下：
  $$P(y|\mathbf{\theta })=\sum _{k=1}^K{\alpha }_k\varphi (y|{\mathbf{\theta }}_{\mathbf{k}})$$
  参数为 $\mathbf{\theta }=({\alpha }_1,{\alpha }_2,...,{\alpha }_k;{\mathbf{\theta }}_1，{\mathbf{\theta }}_2,...,{\mathbf{\theta }}_{\mathbf{k}})$，且其中 ${\alpha }_k$ 是系数（表示从各模型抽取的概率），$\sum _{k=1}^K{\alpha }_k=1，\varphi (y|{\mathbf{\theta }}_{\mathbf{k}})$ 是高斯分布密度，${\mathbf{\theta }}_{\mathbf{k}}=({\mu }_k,{\sigma }_k^2)$，第 $k$ 个分模型的概率密度函数是
  $$\varphi (y|{\mathbf{\theta }}_{\mathbf{k}})=\frac{1}{\sqrt{2}\pi {\sigma }_k}\exp (-\frac{{\left(y-{\mu }_k\right)}^2}{2{\sigma }_k^2})$$

• 假定样本是从三个一维高斯分布中抽取出来的，抽取的概率分别为：
  $$p_1=\frac{1}{2}，p_2=\frac{1}{3}，p_3=\frac{1}{6}$$

• 这三个高斯函数的参数分别为：(${\mu }_1,{\sigma }_1$) = (1,2) , (${\mu }_2,{\sigma }_2$) = (2,1),(${\mu }_3,{\sigma }_3$) = (4,8)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
#下面完成50000次取样
times = 50000
choice = np.random.choice([1,2,3],size = (1,times),p = [1/2,1/3,1/6]).tolist()
sample_list = []
for i in choice[0]:
    if i == 1:
        sample_list.append(np.random.normal(1,2,1).item())
    elif i == 2:
        sample_list.append(np.random.normal(2,1,1).item())
    else:
        sample_list.append(np.random.normal(4,8,1).item())

sample_list = np.array(sample_list)

bins = np.arange(-5,10,0.1)

plt.hist(sample_list,bins)
plt.show()

(2)EM的理论算法流程

• 明确隐变量，写出完全数据的对数似然函数
  观测数据 $y_1,y_2,...,y_N$，反映观测数据 $y_k$ 来自第 $k$ 个分模型的隐变量 ${\gamma }_{jk}$
  $${\gamma }_{jk}=\{\begin{array}{ll}1,& 第j个观测变量来自第k个分模型\\ 0,& 其余情况\end{array}$$
  综合观测数据和未观测数据可以得到完全数据是 $(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jk})$。
  写出完全数据的似然函数
  $$\begin{array}{rl}P(y,\mathbf{\gamma }|\mathbf{\theta })& =\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK}|\mathbf{\theta })\\ & =\prod _{k=1}^K\prod _{j=1}^N({\alpha }_k\varphi (y_j|{\mathbf{\theta }}_{\mathbf{k}}){)}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\varphi (y_j|{\mathbf{\theta }}_{\mathbf{k}}){]}^{{\gamma }_{jk}}\\ & =\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\frac{1}{\sqrt{2}\pi {\sigma }_k}\exp (-\frac{{\left(y_j-{\mu }_k\right)}^2}{2{\sigma }_k^2}){]}^{{\gamma }_{jk}}\end{array}$$
  其中 $n_k=\sum _{j=1}^N{\gamma }_{jk}$，$\sum _{k=1}^K{n}_k=N$
  变成对数形式 $$\log P(y,\gamma |\mathbf{\theta })=\sum _{k=1}^K{n}_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2]$$

• E步：确定 Q 函数
  $$\begin{array}{rl}Q(\theta ,{\theta }^{(i)})& =E[\log P(y,\gamma |\mathbf{\theta })]\\ & =\sum _{k=1}^K{n}_k\log {\alpha }_k+\sum _{j=1}^N{\hat{\gamma }}_{jk}[\log (\frac{1}{\sqrt{2\pi }})-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2]\end{array}$$
  其中 $${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi \left(y_j|{\theta }_k\right)}{\sum _{k=1}^K{\alpha }_k\varphi \left(y_j|{\theta }_k\right)}\left(j=1,2,..,N\right)\left(k=1,2,...,K\right)$$

• M步：确定 ${\theta }^{(i+1)}$
  $${\theta }^{(i+1)}=\underset{\mathbf{\theta }}{\mathrm{argmax}}Q(\mathbf{\theta },{\mathbf{\theta }}^{(i)})$$
  通过求偏导可得
  $${\hat{\mu }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{\hat{\sigma }}_k^2=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{\left(y_j-{\mu }_k\right)}^2}{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{\hat{\alpha }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N}\left(k=1,2,...,K\right)$$

(3)EM的编程算法流程

• EM 的步骤：

  输入：观测数据 $y_1,y_2,...,y_N$，高斯混合模型；

  输出：高斯混合模型参数

• 取参数的初始值并开始迭代。

• E步：依据当前模型参数，计算分模型$k$对观测数据$y_j$的响应度

$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi \left(y_j|{\theta }_k\right)}{\sum _{k=1}^K{\alpha }_k\varphi \left(y_j|{\theta }_k\right)}\left(j=1,2,..,N\right)\left(k=1,2,...,K\right)$$

• M步：计算新一轮的迭代参数

$${\hat{\mu }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{\hat{\sigma }}_k^2=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}{\left(y_j-{\mu }_k\right)}^2}{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{\hat{\alpha }}_k=\frac{\sum _{j=1}^N{\hat{\gamma }}_{jk}}{N}\left(k=1,2,...,K\right)$$

(3)在这个例子中的EM

• 在这里总共有50000个样本，因此$N=50000$。总共有3个高斯分布，因此$K=3$。
• 给出初始迭代值
  $$p_1=\frac{1}{4}，p_2=\frac{1}{2}，p_3=\frac{1}{4}$$
  (${\mu }_1,{\sigma }_1$) = (1,1) , (${\mu }_2,{\sigma }_2$) = (2,2), (${\mu }_3,{\sigma }_3$) = (3,3)
• 得到的迭代收敛值
  得到的概率值为：$p_1=0.17,p_2=0.35,p_3=0.47$
  得到的分布为：(${\mu }_1,{\sigma }_1$) = (1,1) , (${\mu }_2,{\sigma }_2$) = (1,1), (${\mu }_3,{\sigma }_3$) = (2,5)

import scipy.stats
gamma = np.empty([times,3])
p = np.array([1 / 4, 1 / 2, 1 / 4])
mu = np.array([1, 2, 3])
std = np.array([1, 2, 3])
#这里设定迭代次数为1000次
iteration = 1000

# 这里进行 E-step的操作
def E_step():
    calculate_temp = np.empty((times,3))
    for k in range(3):
        calculate_temp[:, k] = p[k] * scipy.stats.norm(mu[k], std[k]).pdf(sample_list)
    divide_temp = np.sum(calculate_temp, axis=1).reshape(times, )
    for k in range(3):
        gamma[:, k] = calculate_temp[:, k] / divide_temp

#这里进行M-step的操作
def M_step():
    for k in range(3):
        # 这里更新mu的值
        mu[k] = np.divide(np.sum(gamma[:,k] * sample_list) ,np.sum(gamma[:,k]))
        # 这里更新std的值
        std[k] = np.sqrt(np.sum(gamma[:,k] * np.square(sample_list - mu[k]))/np.sum(gamma[:,k]))
        # 这里更新p的值
        p[k] = np.sum(gamma[:,k]) / times

        
for cnt in  np.arange(iteration):   
    E_step()
    M_step()
print("得到的概率值为：p1 = {:.2f},p2 = {:.2f},p3 = {:.2f}".format(p[0],p[1],p[2]))
print("得到的分布为：(mu1,std1) = ({},{}) (mu2,std2) = ({},{}) (mu3,std3) = ({},{})".format(mu[0],std[0],mu[1],std[1],mu[2],std[2]))
    

(4)画出EM计算出的结果

import math
EM_lis = np.arange(-5,10,0.1)
EM_draw = np.zeros(EM_lis.shape)

for j in range(3):
    EM_draw = EM_draw + scipy.stats.norm(mu[j],std[j]).pdf(EM_lis)* p[j]
#下面乘以4500只是做一个scale，并没有其他的含义
plt.plot(EM_lis,EM_draw*4500)
plt.hist(sample_list,bins)
plt.show()