参数估计问题
在第一课中,提到使用样本估计模型(比如高斯分布)的参数,并说明了常用的极大似然估计法。假设现在有一枚硬币,但它质地不均匀,导致抛硬币的正面朝上与反面朝上的概率不相等,现在还是想研究正面朝上的概率是多少,所以可以抛硬币 N N N次,正面朝上的次数为 n 1 n_{1} n1,则就使用 n 1 / N n_{1}/N n1/N作为正面朝上概率的估计值。
再举例一个问题,假设已知某个班级内的男同学身高服从正态分布,现在要研究这个正态分布的均值和方差,我们可以随机挑选
N
N
N个男同学的身高数据作为样本,分别统计他们的身高
x
1
,
x
2
,
.
.
.
,
x
N
x_{1},x_{2},...,x_{N}
x1,x2,...,xN,通过极大似然估计法得到均值和方差:
μ
m
l
e
=
1
N
∑
i
=
1
N
x
i
,
σ
m
l
e
2
=
1
N
∑
i
=
1
N
(
x
i
−
μ
m
l
e
)
2
\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}
μmle=N1i=1∑Nxi,σmle2=N1i=1∑N(xi−μmle)2
若考虑无偏性时,需要对方差进行修正:
σ
2
=
1
N
−
1
∑
i
=
1
N
(
x
i
−
μ
m
l
e
)
2
\sigma^{2}=\frac{1}{N-1}\sum_{i=1}^{N}(x_{i}-\mu_{mle})^{2}
σ2=N−11i=1∑N(xi−μmle)2
含隐变量的问题
如果在实际场景中,由于情景的细微改变,极大似然估计将难以解决问题。假设现在有两枚硬币,硬币 A A A和硬币 B B B,它们都是不均匀的,而且两者的正面朝上概率 p A p_{A} pA和 p B p_{B} pB也不相等。每次从这两枚硬币随机摸出一枚投掷(但不知道取出哪一枚),现在依然记录每次投掷的结果是正面还是反面,其中,正面朝上的次数为 N 1 N_{1} N1,反面朝上的次数为 N 2 N_{2} N2,现在想用这些数据直接估计硬币 A A A和硬币 B B B各自正面朝上的概率就变得很困难。
对于统计班级内男同学的身高,如果获得的样本中,既有男同学,又有女同学,当获取 N N N个同学的身高数据后,显然不能直接求出男同学身高的均值和方差,因为男同学,女同学各自会服从不同的正态分布。
以上两个问题,不能直接用极大似然估计参数,因为两个问题都变成了混合模型,不仅仅存在观测变量,还存在隐变量。通常,第 i i i个样本的观测变量记为 x i x_{i} xi,隐变量记为 z i z_{i} zi:
- 对于抛硬币问题,每次抛硬币的正反面记录是观测变量,但某次抛出的硬币是 A A A还是 B B B则属于隐变量;
- 对于身高统计问题,每次得到的身高数据是观测变量,但某个数据是属于男同学还是女同学则为隐变量。
迭代法解决含隐变量的问题
对于混合模型,不能直接用极大似然法估计参数,因此,考虑用迭代法去逐步尝试:即EM算法;
EM算法包含较多计算技巧,较为复杂,因此,本例先不详细展开,而是针对抛硬币问题做简单计算,先感性认识迭代法探索参数的大致过程;
假设有不均匀的硬币 A A A和 B B B,重复5次实验,每次实验都是同样的操作:任取一枚硬币,连续投掷该硬币10次,记录正反面次数。5组实验结果如下:
组别 | 正面次数 | 反面次数 |
---|---|---|
1 | 5 | 5 |
2 | 9 | 1 |
3 | 8 | 2 |
4 | 4 | 6 |
5 | 7 | 3 |
问题为:求出硬币 A A A和硬币 B B B正面朝上的概率 p A p_{A} pA和 p B p_{B} pB。
目前面临的问题是不知道每一组投掷的到底是哪一个硬币,但如果已经知道每组用什么硬币投掷,我们将能快速计算出答案,比如,已知第二组,第三组,第五组使用硬币
A
A
A,而第一组和第四组使用的是硬币
B
B
B,则对于硬币
A
A
A的参数估计可利用所有关于硬币
A
A
A的实验数据求出:
p
A
=
9
+
8
+
7
30
=
0.8
p_{A}=\frac{9+8+7}{30}=0.8
pA=309+8+7=0.8
同理,硬币
B
B
B的参数为:
p
B
=
5
+
4
20
=
0.45
p_{B}=\frac{5+4}{20}=0.45
pB=205+4=0.45
但是并不知道每组实验是用哪个硬币,所以用迭代的步骤反复估计
p
A
p_{A}
pA和
p
B
p_{B}
pB以及
A
,
B
A,B
A,B各自的出现概率
P
A
,
P
B
P_{A},P_{B}
PA,PB。首先随机给定初始值:
p
A
0
=
0.6
,
p
B
0
=
0.5
p_{A}^{0}=0.6,p_{B}^{0}=0.5
pA0=0.6,pB0=0.5
因此,在第一组实验中,记使用硬币
A
A
A的概率为
P
A
1
P_{A1}
PA1,使用硬币
B
B
B的概率为
P
B
1
P_{B1}
PB1,在实验1中,若完全用
A
A
A投掷,出现5正5反的概率为:
(
p
A
0
)
5
(
1
−
p
A
0
)
5
=
0.
6
5
0.
4
5
=
0.000796
(p_{A}^{0})^{5}(1-p_{A}^{0})^{5}=0.6^{5}0.4^{5}=0.000796
(pA0)5(1−pA0)5=0.650.45=0.000796
同理,若完全用
B
B
B投掷,出现5正5反的概率为:
(
p
B
0
)
5
(
1
−
p
B
0
)
5
=
0.
5
5
0.
5
5
=
0.000976
(p_{B}^{0})^{5}(1-p_{B}^{0})^{5}=0.5^{5}0.5^{5}=0.000976
(pB0)5(1−pB0)5=0.550.55=0.000976
所以,第一组实验使用硬币
A
A
A和硬币
B
B
B的概率为:
P
A
1
=
0.
6
5
0.
4
5
0.
6
5
0.
4
5
+
0.
5
5
0.
5
5
=
0.45
P_{A1}=\frac{0.6^{5}0.4^{5}}{0.6^{5}0.4^{5}+0.5^{5}0.5^{5}}=0.45
PA1=0.650.45+0.550.550.650.45=0.45
P
B
1
=
0.
5
5
0.
5
5
0.
6
5
0.
4
5
+
0.
5
5
0.
5
5
=
0.55
P_{B1}=\frac{0.5^{5}0.5^{5}}{0.6^{5}0.4^{5}+0.5^{5}0.5^{5}}=0.55
PB1=0.650.45+0.550.550.550.55=0.55
按照同样的方式,可以分别计算出另外四组实验中,使用硬币
A
A
A和硬币
B
B
B的概率:
组别 | 硬币 A A A | 硬币 B B 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 A 1 = 0.45 P_{A1}=0.45 PA1=0.45和 P B 1 = 0.55 P_{B1}=0.55 PB1=0.55,第一组实验结果为5正5反,则使用硬币 A A A投掷的结果应有:
- 正: 0.45 × 5 = 2.25 0.45\times 5=2.25 0.45×5=2.25
- 反: 0.45 × 5 = 2.25 0.45\times 5=2.25 0.45×5=2.25
使用硬币 B B B投掷的结果为:
- 正: 0.55 × 5 = 2.75 0.55\times 5=2.75 0.55×5=2.75
- 反: 0.55 × 5 = 2.75 0.55\times 5=2.75 0.55×5=2.75
进一步,用同样的流程得到各组实验中,硬币 A A A和硬币 B B B的投掷结果:
组别 | 硬币 A A A | 硬币 B B 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
A
A和硬币
B
B
B正面朝上的概率:
p
A
1
=
21.24
21.24
+
8.56
=
0.71
,
p
B
1
=
11.76
11.76
+
8.44
=
0.58
p_{A}^{1}=\frac{21.24}{21.24+8.56}=0.71,p_{B}^{1}=\frac{11.76}{11.76+8.44}=0.58
pA1=21.24+8.5621.24=0.71,pB1=11.76+8.4411.76=0.58
再重复之前的操作,于是得到
(
p
A
2
,
p
B
2
)
(p_{A}^{2},p_{B}^{2})
(pA2,pB2),依次类推,得到:
(
p
A
3
,
p
B
3
)
,
(
p
A
4
,
p
B
4
)
,
.
.
.
,
(
p
A
T
,
p
B
T
)
(p_{A}^{3},p_{B}^{3}),(p_{A}^{4},p_{B}^{4}),...,(p_{A}^{T},p_{B}^{T})
(pA3,pB3),(pA4,pB4),...,(pAT,pBT)
直到参数
(
p
A
T
,
p
B
T
)
(p_{A}^{T},p_{B}^{T})
(pAT,pBT)收敛;
实例演示
硬币投掷问题的实例演示如下:
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_{A}=0.72225,p_{B}=0.55543 pA=0.72225,pB=0.55543,由于问题较简单,过程容易理解,但是实际问题总是复杂的,最好能总结出一个统一形式的算法,处理这种类似 p t p^{t} pt和 p t + 1 p^{t+1} pt+1的迭代关系,所以人们提出了EM算法,它广泛应用于处理混合模型(包含隐变量的模型),后面的部分内容将会深入分析EM算法。