第七课.含隐变量的参数估计

参数估计问题

在第一课中，提到使用样本估计模型（比如高斯分布）的参数，并说明了常用的极大似然估计法。假设现在有一枚硬币，但它质地不均匀，导致抛硬币的正面朝上与反面朝上的概率不相等，现在还是想研究正面朝上的概率是多少，所以可以抛硬币$N$次，正面朝上的次数为$n_1$，则就使用$n_1/N$作为正面朝上概率的估计值。

再举例一个问题，假设已知某个班级内的男同学身高服从正态分布，现在要研究这个正态分布的均值和方差，我们可以随机挑选$N$个男同学的身高数据作为样本，分别统计他们的身高$x_1,x_2,...,x_N$，通过极大似然估计法得到均值和方差：
$${\mu }_{mle}=\frac{1}{N}\sum _{i=1}^N{x}_i,{\sigma }_{mle}^2=\frac{1}{N}\sum _{i=1}^N(x_i-{\mu }_{mle}{)}^2$$
若考虑无偏性时，需要对方差进行修正：
$${\sigma }^2=\frac{1}{N-1}\sum _{i=1}^N(x_i-{\mu }_{mle}{)}^2$$

含隐变量的问题

如果在实际场景中，由于情景的细微改变，极大似然估计将难以解决问题。假设现在有两枚硬币，硬币$A$和硬币$B$，它们都是不均匀的，而且两者的正面朝上概率$p_A$和$p_B$也不相等。每次从这两枚硬币随机摸出一枚投掷（但不知道取出哪一枚），现在依然记录每次投掷的结果是正面还是反面，其中，正面朝上的次数为$N_1$，反面朝上的次数为$N_2$，现在想用这些数据直接估计硬币$A$和硬币$B$各自正面朝上的概率就变得很困难。

对于统计班级内男同学的身高，如果获得的样本中，既有男同学，又有女同学，当获取$N$个同学的身高数据后，显然不能直接求出男同学身高的均值和方差，因为男同学，女同学各自会服从不同的正态分布。

以上两个问题，不能直接用极大似然估计参数，因为两个问题都变成了混合模型，不仅仅存在观测变量，还存在隐变量。通常，第$i$个样本的观测变量记为$x_i$，隐变量记为$z_i$：

• 对于抛硬币问题，每次抛硬币的正反面记录是观测变量，但某次抛出的硬币是$A$还是$B$则属于隐变量；
• 对于身高统计问题，每次得到的身高数据是观测变量，但某个数据是属于男同学还是女同学则为隐变量。

迭代法解决含隐变量的问题

对于混合模型，不能直接用极大似然法估计参数，因此，考虑用迭代法去逐步尝试：即EM算法；

EM算法包含较多计算技巧，较为复杂，因此，本例先不详细展开，而是针对抛硬币问题做简单计算，先感性认识迭代法探索参数的大致过程；

假设有不均匀的硬币$A$和$B$，重复5次实验，每次实验都是同样的操作：任取一枚硬币，连续投掷该硬币10次，记录正反面次数。5组实验结果如下：

组别	正面次数	反面次数
1	5	5
2	9	1
3	8	2
4	4	6
5	7	3

问题为：求出硬币$A$和硬币$B$正面朝上的概率$p_A$和$p_B$。

目前面临的问题是不知道每一组投掷的到底是哪一个硬币，但如果已经知道每组用什么硬币投掷，我们将能快速计算出答案，比如，已知第二组，第三组，第五组使用硬币$A$，而第一组和第四组使用的是硬币$B$，则对于硬币$A$的参数估计可利用所有关于硬币$A$的实验数据求出：
$$p_A=\frac{9+8+7}{30}=0.8$$
同理，硬币$B$的参数为：
$$p_B=\frac{5+4}{20}=0.45$$
但是并不知道每组实验是用哪个硬币，所以用迭代的步骤反复估计$p_A$和$p_B$以及$A,B$各自的出现概率$P_A,P_B$。首先随机给定初始值：
$$p_A^0=0.6,p_B^0=0.5$$
因此，在第一组实验中，记使用硬币$A$的概率为$P_{A1}$，使用硬币$B$的概率为$P_{B1}$，在实验1中，若完全用$A$投掷，出现5正5反的概率为：
$$(p_A^0{)}^5(1-p_A^0{)}^5=0.6^{5}0.4^5=0.000796$$
同理，若完全用$B$投掷，出现5正5反的概率为：
$$(p_B^0{)}^5(1-p_B^0{)}^5=0.5^{5}0.5^5=0.000976$$
所以，第一组实验使用硬币$A$和硬币$B$的概率为：
$$P_{A1}=\frac{0.6^{5}0.4^5}{0.6^{5}0.4^5+0.5^{5}0.5^5}=0.45$$
$$P_{B1}=\frac{0.5^{5}0.5^5}{0.6^{5}0.4^5+0.5^{5}0.5^5}=0.55$$
按照同样的方式，可以分别计算出另外四组实验中，使用硬币$A$和硬币$B$的概率：

组别	硬币$A$	硬币$B$
1	0.45	0.55
2	0.80	0.20
3	0.73	0.27
4	0.35	0.65
5	0.65	0.35

现在考虑第1组实验的$P_{A1}=0.45$和$P_{B1}=0.55$，第一组实验结果为5正5反，则使用硬币$A$投掷的结果应有：

• 正：$0.45\times 5=2.25$
• 反：$0.45\times 5=2.25$

使用硬币$B$投掷的结果为：

• 正：$0.55\times 5=2.75$
• 反：$0.55\times 5=2.75$

进一步，用同样的流程得到各组实验中，硬币$A$和硬币$B$的投掷结果：

组别	硬币$A$	硬币$B$
1	正：2.25，反：2.25	正：2.75，反：2.75
2	正：7.2，反：0.8	正：1.8，反：0.2
3	正：5.84，反：1.46	正：2.16，反：0.54
4	正：1.4，反：2.1	正：2.6，反：3.9
5	正：4.55，反：1.95	正：2.45，反：1.05
合计	正：21.24，反：8.56	正：11.76，反：8.44

基于以上表，重新估计硬币$A$和硬币$B$正面朝上的概率：
$$p_A^1=\frac{21.24}{21.24+8.56}=0.71,p_B^1=\frac{11.76}{11.76+8.44}=0.58$$
再重复之前的操作，于是得到$(p_A^2,p_B^2)$，依次类推，得到：
$$(p_A^3,p_B^3),(p_A^4,p_B^4),...,(p_A^T,p_B^T)$$
直到参数$(p_A^T,p_B^T)$收敛；

实例演示

硬币投掷问题的实例演示如下：

import numpy as np
# 二项分布
from scipy.stats import binom
"""
Bernoulli：伯努利分布，是关于布尔变量的概率分布
Binomial：二项分布，重复多次的独立伯努利实验
"""

#一轮迭代处理
def single_iter(theta_priors, observations):
    """
    param observations：五组试验的观测值
    param theta_priors：上一轮迭代形成的 p_A 和 p_B
    """
    counts = {'A': {'H': 0, 'T': 0}, 'B': {'H': 0, 'T': 0}}
    theta_A = theta_priors[0]
    theta_B = theta_priors[1]
    
    #迭代计算每组试验的数据
    for observation in observations:
        len_observation = len(observation)
        num_heads = observation.sum()
        num_tails = len_observation-num_heads
        
        # binom.pmf(k,n,p)返回n次伯努利实验中结果为目标事件的概率，单次实验目标事件的概率为p，目标事件出现了k次
        P_A = binom.pmf(num_heads, len_observation, theta_A)
        P_B = binom.pmf(num_heads, len_observation, theta_B)
        # 计算出硬币A和硬币B各自出现的概率
        weight_A = P_A / (P_A + P_B)
        weight_B = P_B / (P_A + P_B)
        # 更新在当前硬币A和硬币B各自出现的概率下，硬币A和硬币B各自的正反面次数
        counts['A']['H'] += weight_A * num_heads
        counts['A']['T'] += weight_A * num_tails
        counts['B']['H'] += weight_B * num_heads
        counts['B']['T'] += weight_B * num_tails

    #重新估计新一轮的硬币A和硬币B正面向上的概率
    new_theta_A = counts['A']['H'] / (counts['A']['H'] + counts['A']['T'])
    new_theta_B = counts['B']['H'] / (counts['B']['H'] + counts['B']['T'])
    return [new_theta_A,new_theta_B]


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

prior = [0.6, 0.5] #设定初始的参数值
iteration = 0
iterations = 10000 #最多的迭代次数
tol = 1e-6         #判断参数收敛的阈值
while iteration < iterations:
    new_prior = single_iter(prior,observations)
    print(new_prior)
    delta_change = np.abs(prior[0] - new_prior[0])
    if delta_change < tol:
        break
    else:
        prior = new_prior
        iteration += 1

print('迭代结束，总共迭代轮数{}'.format(iteration))
print('最终估计的参数{}'.format(new_prior))

"""
[0.683267379440028, 0.5794860533160178]
...
[0.7222502854992598, 0.5554380899384829]
迭代结束，总共迭代轮数36
最终估计的参数[0.7222502854992598, 0.5554380899384829]
"""

经过36次迭代，得到收敛的结果为：$p_A=0.72225,p_B=0.55543$，由于问题较简单，过程容易理解，但是实际问题总是复杂的，最好能总结出一个统一形式的算法，处理这种类似$p^t$和$p^{t+1}$的迭代关系，所以人们提出了EM算法，它广泛应用于处理混合模型（包含隐变量的模型），后面的部分内容将会深入分析EM算法。