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

最新推荐文章于 2022-10-24 14:27:42 发布

tzc_fly

最新推荐文章于 2022-10-24 14:27:42 发布

阅读量666

点赞数 1

分类专栏：随机过程与概率图模型

本文链接：https://blog.csdn.net/qq_40943760/article/details/115660070

版权

随机过程与概率图模型专栏收录该内容

23 篇文章 17 订阅

订阅专栏

参数估计问题

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