python EM算法2

1硬币问题

先看一个抛硬币问题，如果我们有A和B两个不均匀硬币，选择任意一个硬币抛10次(这里我们知道选择是的哪一个硬币)，共计选择5次。正面记为H，背面记为T。记录实验结果，求A和B再抛正面向上的概率？

使用极大似然估计(Maximum likelihood)来算：

• 统计出每次实验，正反面的次数
• 多次实验结果相加
• 相除得到结果，P(A)=0.8,P(B)=0.45

但是在实际过程中，很有可能我们只知道有两个硬币，不知道每次选择的哪一个硬币，问是否能求出每个硬币抛出正面的概率？

是不是第一感觉这个问题很无解，我都不知道选择是哪个硬币，我怎么求解啊？
事实上是可以求出的。这里使用的时EM算法。

使用期望最大值(Expectation maximization,EM)算法来算:

• 假设$\hat{\theta}_A^{(0)}=0.6,\hat{\theta}_B^{(0)}=0.5$
• 统计每次的实验结果，记录正反面
• 通过贝叶斯公式，估计每次实验选择的A硬币或是B硬币的概率
• 依据计算出的选择硬币概率得到该概率下的正反面结果
• 5次平均得到$\hat{\theta}_A^{(1)}\approx 0.71,\hat{\theta}_B^{(1)}\approx 0.58$
• 重复上面的过程，例如迭代10次后，得到$\hat{\theta}_A^{(2)}\approx 0.8,\hat{\theta}_B^{(2)}\approx 0.52$
• $${\hat{\theta }}_A^{(10)},{\hat{\theta }}_B^{(10)}$`就是使用EM算法计算出的概率值

这里比较难理解的是：如何利用贝叶斯公式计算每次实验选择的A硬币或是B硬币的概率？

---

那么先看下面这个例子：

假如现在有射击运动员甲和乙，甲乙射击靶心的概率为$\hat{\theta}_{\text{甲}}=0.9$，$${\hat{\theta }}_{\text{乙}}=0.2$`，如果现在有一组实验结果为

中，不中，中，中，中 问这次是谁射击的？

**直观上的来看，非常大的概率是甲射击的，但是也有可能是乙走狗屎运了。那么该如何从概率的角度计算是谁射击的？ **

首先我们知道选择甲和乙的概率为p(甲)=p(乙)=0.5(先验概率),本次实验记为E。通过贝叶斯公式

$$p(\text{甲}|\text{击中})=\frac{p(\text{击中}|\text{甲})p(\text{甲})}{p(\text{击中})}$$

$$=\frac{p(\text{击中}|\text{甲})p(\text{甲})}{p(\text{击中}|\text{甲})p(\text{甲})+p(\text{击中}|\text{乙})p(\text{乙})}$$

$$\approx 0.98$`

故本次实验有98%可能是A射击的，2%的可能是B射击的。

---

有了上面的案例，我们再回到抛硬币的问题上，由贝叶斯公式:

p(A\vert E)=\frac{p(E|A)P(A)}{p(E)}

A 为选用硬币A，E为本次实验。而选择两个硬币的概率是相同的：$p(A)=p(B)=0.5$,且$p(E|A)=0.6^5\times 0.4^5\approx 0.0008,P(E|B)=0.5^5\times 0.5^5\approx 0.001$

\begin{array}{lll}
p(A|E)&=&\frac{p(E|A)p(A)}{p(E)}\\
&=&\frac{p(E|A)p(A)}{p(E|A)p(A)+p(E|B)p(B)}\\
&\approx & 0.45
\end{array}

有$0.45\times 10\times 0.5\approx 2.2$,故第一次实验结果平均下来，有2.22.2个A硬币正面的可能。

最后将5次结果平均得到新的A,B抛硬币正面估计值$\hat{\theta}_A^{(1)}\approx 0.71,\hat{\theta}_B^{(1)}\approx 0.58$,，这是我们第一次迭代的值(这就是一次学习过程)，照着这个流程迭代多次，得到最后的估测值。

上面我们计算出每次实验中是抛A或抛B的概率值就是隐变量.这个过程就是EM算法的简单案例。

2形式化EM算法

2.1 算法理论

`$p(x)=\sum _{j=1}^k{\varphi }_{j}N(x|{\mu }_j,{\Sigma }_j)$$

`$N(x|{\mu }_j,{\Sigma }_j)=\frac{1}{(2\pi {)}^{(}n/2)|{\Sigma }_j{|}^{(}1/2)}\exp (-\frac{1}{2}(x-{\mu }_j{)}^T{\Sigma }_j^{-1}(x-{\mu }_j))$$

$$\sum _{j=1}^k{\varphi }_j=1$`

这是一个非凸问题，只能求局部极值，一般采用EM算法。

步骤如下:

(1).随机化` { ( ϕ j , μ j , Σ j ) ∣ 1 ≤ j ≤ k } \{(\phi_j,\mu_j,\Sigma_j)\vert 1\leq j \leq k\} {(ϕj​,μj​,Σj​)∣1≤j≤k}$

(2).E步：计算

$${\omega }_j^{(i)}=Q_i(z^{(i)}=j)=P(z^{(i)}=j|x^{(i)};\varphi ,\mu ,\Sigma )$$

$$=\frac{{\varphi }_{j}N(x^{(i)}|{\mu }_j,{\Sigma }_j)}{\sum _{j=1}^{k}N(x^{(i)}|{\mu }_j,{\Sigma }_j)}$$

上式表示：第i个样本落在第j个高斯的概率.

(3). M步：

$${\mu }_j:=\frac{{\sum }_{i=1}^m{\omega }_{\ell }^{(i)}{x}^{(i)}}{{\sum }_{i=1}^m{\omega }_{\ell }^{(i)}}$$

$${\Sigma }_j=\frac{\sum _{i=1}^m{\omega }_j^{(i)}(x^{(i)}-{\mu }_j)(x^{(i)}-{\mu }_j{)}^T}{\sum _{i=1}^m{\omega }_j^{(i)}}$$

$${\varphi }_j:=\frac{1}{m}\sum _{i=1}^m{\omega }_j^{(i)}$$

(4).回到(2)直至收敛.

2.2 python实现
以西瓜数据集4.0为例:

编号	密度	含糖率
1	0.697	0.46
2	0.774	0.376
3	0.634	0.264
4	0.608	0.318
5	0.556	0.215
6	0.403	0.237
7	0.481	0.149
8	0.437	0.211
9	0.666	0.091
10	0.243	0.267
11	0.245	0.057
12	0.343	0.099
13	0.639	0.161
14	0.657	0.198
15	0.36	0.37
16	0.593	0.042
17	0.719	0.103
18	0.359	0.188
19	0.339	0.241
20	0.282	0.257
21	0.748	0.232
22	0.714	0.346
23	0.483	0.312
24	0.478	0.437
25	0.525	0.369
26	0.751	0.489
27	0.532	0.472
28	0.473	0.376
29	0.725	0.445
30	0.446	0.459

import numpy as np
# 预处理数据
def loadData(filename):
     dataSet = []
     fr = open(filename)
     for line in fr.readlines():
         curLine = line.strip().split(' ')
         fltLine = list(map(float, curLine))
         dataSet.append(fltLine)
     return dataSet
# 计算高斯函数
def Gaussian(data,mean,cov):
    dim = np.shape(cov)[0]   # 计算维度
    covdet = np.linalg.det(cov) # 计算|cov|
    covinv = np.linalg.inv(cov) # 计算cov的逆
    if covdet==0:              # 以防行列式为0
        covdet = np.linalg.det(cov+np.eye(dim)*0.01)
        covinv = np.linalg.inv(cov+np.eye(dim)*0.01)
    m = data - mean
    z = -0.5 * np.dot(np.dot(m, covinv),m)    # 计算exp()里的值
    return 1.0/(np.power(np.power(2*np.pi,dim)*abs(covdet),0.5))*np.exp(z)  # 返回概率密度值

# 获取最初的聚类中心
def GetInitialMeans(data,K,criterion):
    dim = data.shape[1]  # 数据的维度
    means = [[] for k in range(K)] # 存储均值
    minmax=[]
    for i in range(dim):
        minmax.append(np.array([min(data[:,i]),max(data[:,i])]))  # 存储每一维的最大最小值
    minmax=np.array(minmax)
    while True:
        for i in range(K):
            means[i]=[]
            for j in range(dim):
                 means[i].append(np.random.random()*(minmax[j][1]-minmax[j][0])+minmax[j][0] ) #随机产生means
            means[i]=np.array(means[i])

        if isdistance(means,criterion):
            break
    return means
# 用于判断初始聚类簇中的means是否距离离得比较近
def isdistance(means,criterion=0.03):
     K=len(means)
     for i in range(K):
         for j in range(i+1,K):
             if criterion>np.linalg.norm(means[i]-means[j]):
                 return False
     return True

dataSet=loadData('d:/watermelon4.txt')
means=GetInitialMeans(np.array(dataSet),3,0.03)

# K均值算法，估计大约几个样本属于一个GMM
def Kmeans(data,K):
    N = data.shape[0]  # 样本数量
    dim = data.shape[1]  # 样本维度
    means = GetInitialMeans(data,K,15)
    means_old = [np.zeros(dim) for k in range(K)]
    # 收敛条件
    while np.sum([np.linalg.norm(means_old[k] - means[k]) for k in range(K)]) > 0.01:
        means_old = cp.deepcopy(means)
        numlog = [0] * K  # 存储属于某类的个数
        sumlog = [np.zeros(dim) for k in range(K)]
        # E步
        for i in range(N):
            dislog = [np.linalg.norm(data[i]-means[k]) for k in range(K)]
            tok = dislog.index(np.min(dislog))
            numlog[tok]+=1         # 属于该类的样本数量加1
            sumlog[tok]+=data[i]   # 存储属于该类的样本取值

        # M步
        for k in range(K):
            means[k]=1.0 / numlog[k] * sumlog[k]
    return means

def GMM(data,K):
    N = data.shape[0]
    dim = data.shape[1]
    means= Kmeans(data,K)
    convs=[0]*K
    # 初始方差等于整体data的方差
    for i in range(K):
        convs[i]=np.cov(data.T)
    pis = [1.0/K] * K
    gammas = [np.zeros(K) for i in range(N)]
    loglikelyhood = 0
    oldloglikelyhood = 1

    while np.abs(loglikelyhood - oldloglikelyhood) > 0.0001:
        oldloglikelyhood = loglikelyhood

        # E步
        for i in range(N):
            res = [pis[k] * Gaussian(data[i],means[k],convs[k]) for k in range(K)]
            sumres = np.sum(res)
            for k in range(K):           # gamma表示第n个样本属于第k个混合高斯的概率
                gammas[i][k] = res[k] / sumres
        # M步
        for k in range(K):
            Nk = np.sum([gammas[n][k] for n in range(N)])  # N[k] 表示N个样本中有多少属于第k个高斯
            pis[k] = 1.0 * Nk/N
            means[k] = (1.0/Nk)*np.sum([gammas[n][k] * data[n] for n in range(N)],axis=0)
            xdiffs = data - means[k]
            convs[k] = (1.0/ Nk)*np.sum([gammas[n][k]* xdiffs[n].reshape(dim,1) * xdiffs[n] for  n in range(N)],axis=0)
        # 计算最大似然函数
        loglikelyhood = np.sum(
            [np.log(np.sum([pis[k] * Gaussian(data[n], means[k], convs[k]) for k in range(K)])) for n in range(N)])
        print(means)
        print(loglikelyhood)
GMM(np.array(dataSet),3)