文章目录
一、EM算法函数
1.1 EM算法问题导入
调查我们学校的男生和女生的身高分布。 假设你在校园里随便找了100个男生和100个女生。调查学生的身高分布。将他们按照性别划分为两组,然后先统计抽样得到的100个男生的身高。假设他们的身高是服从正态分布的。但是这个分布的均值μ和方差
σ
2
σ^2
σ2我们不知道,这两个参数就是我们要估计的。
记作
θ
=
[
μ
,
σ
]
T
。
θ=[μ,σ]^T。
θ=[μ,σ]T。
问题分析:设样本集 X = x 1 , x 2 , … , x N , X=x_1,x_2,…,x_N, X=x1,x2,…,xN,(N=100), p ( x i ∣ θ ) p(x_i|θ) p(xi∣θ)为概率密度函数,表示抽到男生 x i x_i xi的概率。
由于100个样本之间独立同分布,所以我同时抽到这100个男生的概率就是他们各自概率的乘积,也就是样本集X中各个样本的联合概率,用 似然函数
L
(
θ
)
\mathbf{L(θ)}
L(θ) 表示:
L
(
θ
)
=
L
(
x
1
,
x
2
,
…
,
x
n
;
θ
)
=
∏
i
=
1
n
p
(
x
i
;
θ
)
θ
∈
Θ
L(θ) = L(x_1,x_2,…,x_n;θ) = \prod _{i=1}^{n}p(x_i;θ)~~~~~θ \in \Theta
L(θ)=L(x1,x2,…,xn;θ)=i=1∏np(xi;θ) θ∈Θ
找到一个参数
θ
θ
θ,使得抽到
X
X
X 这组样本的概率最大,也就是说需要其对应的似然函数
L
(
θ
)
L(θ)
L(θ)最大。满足条件的
θ
θ
θ 叫做
θ
θ
θ的最大似然估计量:
θ
^
=
a
r
g
max
L
(
θ
)
\hat{\theta } = arg\max L(\theta )
θ^=argmaxL(θ)
1.2 求解最大似然函数 估计值
- 对 似然函数
L
(
θ
)
L(θ)
L(θ) 取对数:
H ( θ ) = ln ( θ ) = ln ∏ i = 1 n p ( x i ; θ ) = ∑ i = 1 n ln p ( x i ; θ ) H(θ) = \ln(θ) = \ln \prod _{i=1}^{n}p(x_i;θ) = \sum _{i=1}^{n}\ln~p(x_i;θ) H(θ)=ln(θ)=lni=1∏np(xi;θ)=i=1∑nln p(xi;θ) - 对上式求导,令导数为 0 ,得到似然方程。
- 求解似然方程,得到的参数 θ θ θ 即为所求。
附 符合高斯分布 概率密度函数
若给定一组样本
x
1
,
x
2
…
x
n
,
x_1,x_2…x_n,
x1,x2…xn,已知它们来自于高斯分布
N
(
μ
,
σ
)
N(μ,σ)
N(μ,σ),试估计参数
μ
,
σ
μ,σ
μ,σ。
f
(
x
)
=
1
2
π
σ
e
−
(
x
−
μ
)
2
2
σ
2
(
高
斯
分
布
的
概
率
密
度
函
数
)
f(x) = \frac{1}{\sqrt{2\pi }\sigma }e^{-\frac{(x-\mu )^2}{2\sigma^2 }}~~~~~~~(高斯分布的概率密度函数)
f(x)=2πσ1e−2σ2(x−μ)2 (高斯分布的概率密度函数)
此时我们已经获得符合高斯分布的概率密度函数:
对
μ
μ
μ 最大化,我们得到最大似然解(即样本均值):
μ
=
1
N
∑
n
=
1
N
x
n
\mu = \frac 1 N \sum_{n=1}^N x_n
μ=N1n=1∑Nxn
对
σ
2
σ^2
σ2 最大化,我们得到(样本方差。):
σ
2
=
1
N
∑
n
=
1
N
(
x
n
−
μ
)
2
\sigma^2 = \frac 1 N \sum_{n=1}^N (x_n-\mu)^2
σ2=N1n=1∑N(xn−μ)2
但是这个解不是无偏的,由
E
[
x
2
]
=
E
[
x
]
2
+
var
[
x
]
=
μ
2
+
σ
2
\mathbb E[x^2] = \mathbb E[x]^2 + \text{var}[x] = \mu^2 + \sigma^2
E[x2]=E[x]2+var[x]=μ2+σ2 我们可以计算它们的期望:
E
[
μ
]
=
1
N
∑
n
=
1
N
E
[
x
n
]
=
μ
E
[
σ
2
]
=
E
[
1
N
∑
n
=
1
N
(
x
n
−
1
N
∑
m
=
1
N
x
m
)
]
=
1
N
∑
i
=
1
N
E
[
x
n
2
−
2
N
x
n
∑
m
=
1
x
m
+
1
N
2
∑
m
=
1
N
∑
l
=
1
N
x
m
x
l
]
=
(
μ
2
+
σ
2
)
−
2
(
μ
2
+
1
N
σ
2
)
+
μ
2
+
1
N
σ
2
=
(
N
−
1
N
)
σ
2
\begin{aligned} \mathbb E[\mu] & = \frac 1 N \sum_{n=1}^N \mathbb E[x_n] = \mu \\ \mathbb E[\sigma^2] & = \mathbb E \left[\frac{1}{N} \sum_{n=1}^N(x_n-\frac 1 N \sum_{m=1}^N x_m)\right] \\ & = \frac 1 N \sum_{i=1}^N \mathbb E \left[x_n^2-\frac 2 N x_n\sum_{m=1} x_m + \frac{1}{N^2} \sum_{m=1}^N\sum_{l=1}^N x_mx_l\right] \\ & = (\mu^2 + \sigma^2) - 2 (\mu^2+ \frac 1 N \sigma^2) + \mu^2+ \frac 1 N \sigma^2 \\ & = \left(\frac{N-1}{N}\right)\sigma^2 \end{aligned}
E[μ]E[σ2]=N1n=1∑NE[xn]=μ=E[N1n=1∑N(xn−N1m=1∑Nxm)]=N1i=1∑NE[xn2−N2xnm=1∑xm+N21m=1∑Nl=1∑Nxmxl]=(μ2+σ2)−2(μ2+N1σ2)+μ2+N1σ2=(NN−1)σ2
用到了:
- 当 m = n m=n m=n 时, E [ x m x n ] = E [ x n 2 ] = μ 2 + σ 2 \mathbb E[x_mx_n]=\mathbb E[x_n^2] = \mu^2 + \sigma^2 E[xmxn]=E[xn2]=μ2+σ2
- 当 m ≠ n m\neq n m=n 时, E [ x m x n ] = E [ x n ] E [ x m ] = μ 2 \mathbb E[x_mx_n]=\mathbb E[x_n]\mathbb E[x_m] = \mu^2 E[xmxn]=E[xn]E[xm]=μ2
因此,方差的一个无偏估计为:
σ
~
2
=
N
N
−
1
σ
2
=
1
N
−
1
∑
n
=
1
N
(
x
n
−
μ
)
2
\tilde \sigma^2 = \frac{N}{N-1}\sigma^2=\frac{1}{N-1}\sum_{n=1}^N(x_n-\mu)^2
σ~2=N−1Nσ2=N−11n=1∑N(xn−μ)2
随着 N N N 的增大,方差估计的误差也随之增大。
附 混合高斯模型
-
直观理解猜测GMM的参数估计:
-
随机变量X是有K个高斯分布混合而成,取各个高斯分布的概率为π1π2… πK,
第i个高斯分布的均值为μi,方差为Σi。
若观测到随机变量X的一系列样本x1,x2,…,xn,试估计参数π,μ,Σ。在学习混合高斯模型的推论时,使我想起蒙特卡洛方法。虽然两者相差很大,但有一点相通。
混合高斯模型:通过数据点的分布来反推公式中的参数
蒙特卡洛方法:通过落入特定形状数据点的概率来反推落入特定形状面积。
关于混合高斯模型(GMM)的参数估计 请点击
二、Jensen不等式
lazy Statistician规则:
设 y 是随机变量 x 的函数, y = g ( x ) ( g 是 连 续 函 数 ) y = g(x)(g是连续函数) y=g(x)(g是连续函数),那么:
(1) x 是离散随机变量,分布为
p
(
x
=
x
k
)
=
p
k
,
k
∈
N
∗
,
p(x = x_k) = p_k,k \in N^*,
p(x=xk)=pk,k∈N∗,若
∑
k
=
1
∞
g
(
x
)
p
k
\sum_{k=1}^{\infty}~g(x)~p_k
∑k=1∞ g(x) pk绝对收敛,则:
E
(
y
)
=
E
∣
g
(
x
)
∣
=
∑
k
=
1
∞
g
(
x
k
)
p
k
E(y) = E|g(x)| = \sum_{k=1}^{\infty }g(x_k)p_k
E(y)=E∣g(x)∣=k=1∑∞g(xk)pk
(2) x 是连续随机变量,分布概率为
f
(
x
)
f(x)
f(x),若
∫
−
∞
∞
g
(
x
)
f
(
x
)
d
x
\int_{-\infty}^{\infty }g(x)f(x)dx
∫−∞∞g(x)f(x)dx绝对收敛,则:
E
(
y
)
=
E
∣
g
(
x
)
∣
=
∫
−
∞
∞
g
(
x
)
f
(
x
)
d
x
E(y) = E|g(x)| =\int_{-\infty}^{\infty }g(x)f(x)dx
E(y)=E∣g(x)∣=∫−∞∞g(x)f(x)dx
三、算法流程
3.1 数据及参数
观测数据:观测到的随机变量X的样本:
X
=
(
x
1
,
.
.
.
,
x
n
)
X = (x_1,..., x_n)
X=(x1,...,xn)
隐含变量:未观测到的随机变量Z的值:
Z
=
(
z
1
,
.
.
.
,
z
n
)
Z = (z_1,..., z_n)
Z=(z1,...,zn)
完整数据:包含观测到的随机变量
X
X
X 和隐含变量
Z
Z
Z 的数据:
Y
=
(
X
,
Z
)
Y = (X, Z)
Y=(X,Z)
Y
=
(
(
x
1
,
z
1
)
,
.
.
.
,
(
x
n
,
z
n
)
)
Y = ((x_1, z_1),..., (x_n, z_n))
Y=((x1,z1),...,(xn,zn))
3.2 EM算法的推导
EM算法是从含有隐含变量的数据(完整数据)中计算极大似然估计。Z为隐含变量,则从可观测数据入手,对参数进行极大似然估计。
-
根据边缘分布列的定义:
∑ j = 1 + ∞ P ( X = x i , Y = y j ) = P ( X = x i ) \sum_{j=1}^{+ \infty}P(X=x_i,Y=y_j) = P(X=x_i) j=1∑+∞P(X=xi,Y=yj)=P(X=xi)改写成 L ( θ ) L(\theta) L(θ) 并取对数似然函数: ℓ ( θ ) = ∑ i = 1 ln p ( x ( i ) ; θ ) = ∑ i ln ∑ z ( i ) p ( x ( i ) , z ( i ) ; θ ) \begin{aligned} \ell (\theta)&=\sum_{i=1} \ln p(x^{(i)};\theta) \\ &=\sum_{i} \ln \sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta) \end{aligned} ℓ(θ)=i=1∑lnp(x(i);θ)=i∑lnz(i)∑p(x(i),z(i);θ)
将 p ( x ( i ) ) p(x^{(i)}) p(x(i)) 用边缘分布列反向拆分为联合分布。
这里, z z z 是隐随机变量,直接找到参数的估计是很困难的。我们建立 l ( θ ) l(\theta) l(θ)的下界,并求该下界的最大值,重复这个过程,直到收敛到局部最大值。
-
定义隐含变量 Z Z Z 的分布 Q i Q_i Qi:
令 Q i Q_i Qi 是z的某一个分布,满足 ∑ z Q i z ( i ) = 1 \mathbf{\sum_zQ_iz^{(i)}=1} ∑zQiz(i)=1, Q i ⩾ 0 \mathbf{Q_i \geqslant 0} Qi⩾0,有:
∑ i ln p ( x ( i ) ; θ ) = ∑ i ln ∑ z ( i ) p ( x ( i ) , z ( i ) ; θ ) = ∑ i ln ∑ z ( i ) Q i ( z ( i ) ) p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) ⩾ ∑ i ∑ z ( i ) Q i ( z ( i ) ) ln p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) \begin{aligned} \sum_{i} \ln p(x^{(i)};\theta)&=\sum_{i}\ln \sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta) \\ &=\sum_{i}\ln \sum_{z^{(i)}} Q_i(z^{(i)}) \frac{p(x^{(i)},z^{(i)};\theta)}{ Q_i(z^{(i)})} \\ &\geqslant \sum_{i}\sum_{z^{(i)}} Q_i(z^{(i)}) \ln \frac{p(x^{(i)},z^{(i)};\theta)}{ Q_i(z^{(i)})} \end{aligned} i∑lnp(x(i);θ)=i∑lnz(i)∑p(x(i),z(i);θ)=i∑lnz(i)∑Qi(z(i))Qi(z(i))p(x(i),z(i);θ)⩾i∑z(i)∑Qi(z(i))lnQi(z(i))p(x(i),z(i);θ)
上述不等式是根据Jensen不等式的来,其中 ln x \ln x lnx为凹函数。为了使上式等号成立:
p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) = C p ( x ( i ) , z ( i ) ; θ ) = C ( Q i ( z ( i ) ) ) ∑ z p ( x ( i ) , z ( i ) ; θ ) = C ( ∑ z Q i ( z ( i ) ) ) \begin{aligned} \frac{p(x^{(i)},z^{(i)};\theta)}{ Q_i(z^{(i)})} &= C \\ &\\ p(x^{(i)},z^{(i)};\theta)&=C(Q_i(z^{(i)})) \\ &\\ \sum_zp(x^{(i)},z^{(i)};\theta)&=C(\sum_zQ_i(z^{(i)})) \end{aligned} Qi(z(i))p(x(i),z(i);θ)p(x(i),z(i);θ)z∑p(x(i),z(i);θ)=C=C(Qi(z(i)))=C(z∑Qi(z(i)))
满足 ∑ z Q i z ( i ) = 1 \mathbf{\sum_zQ_iz^{(i)}=1} ∑zQiz(i)=1, Q i ⩾ 0 \mathbf{Q_i \geqslant 0} Qi⩾0,由上式推导出:
∑ z p ( x ( i ) , z ( i ) ; θ ) = C \sum_zp(x^{(i)},z^{(i)};\theta)=C z∑p(x(i),z(i);θ)=C
进一步推导:
Q i ( z ( i ) ) = p ( x ( i ) , z ( i ) ; θ ) ∑ z p ( x ( i ) , z ( i ) ; θ ) = p ( x ( i ) , z ( i ) ; θ ) p ( x ( i ) ; θ ) = p ( z ( i ) ∣ x ( i ) ; θ ) \begin{aligned} Q_i(z^{(i)}) &= \frac{p(x^{(i)},z^{(i)};\theta)}{\sum_zp(x^{(i)},z^{(i)};\theta)} \\ &\\ &=\frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \\ &\\ &=p(z^{(i)}|x^{(i)};\theta) \end{aligned} Qi(z(i))=∑zp(x(i),z(i);θ)p(x(i),z(i);θ)=p(x(i);θ)p(x(i),z(i);θ)=p(z(i)∣x(i);θ)
最终我们得到了EM算法的一般形式,一般化形式的思想是:
在 E-step:找到对于当前参数 θ \theta θ 使不等式成立的 Q Q Q分布。
在 M-step:对似然函数下界进行极大似然估计,得到新的参数,形式化表述为:
E-step Q i ( z ( i ) ) = p ( z ( i ) ∣ x ( i ) ; θ ) M-step θ = a r g max θ ∑ i = 1 n ∑ z ( i ) Q i ( z ( i ) ) ln p ( x ( i ) ) , z ( i ) ; θ Q i ( z ( i ) ) \begin{aligned} \textbf{E-step} ~~~~~~~& Q_i(z^{(i)}) = p(z^{(i)}|x^{(i)};\theta)\\ \textbf{M-step} ~~~~~~~& \theta =arg \max_{\theta} \sum_{i=1}^n\sum_{z^{(i)}} Q_i(z^{(i)}) \ln \frac{p(x^{(i)}),z^{(i)};\theta}{Q_i(z^{(i)})} \end{aligned} E-step M-step Qi(z(i))=p(z(i)∣x(i);θ)θ=argθmaxi=1∑nz(i)∑Qi(z(i))lnQi(z(i))p(x(i)),z(i);θ
3.1 算法思路
-
直线式迭代优化的路径:
-
可以看到每一步都会向最优值前进一步,而且前进路线是平行于坐标轴的,因为每一步只优化一个变量。
这犹如在x-y坐标系中找一个曲线的极值,然而曲线函数不能直接求导,因此什么梯度下降方法就不适用了。但固定一个变量后,另外一个可以通过求导得到,因此可以使用坐标上升法,一次固定一个变量,
对另外的求极值,最后逐步逼近极值。
对应到EM上, E步:固定θ,优化Q;
M步:固定Q,优化θ;
交替将极值推向最大。 (如图)
3.2 EM 算法示例:
● 假设现在有两枚硬币1和2,随机抛掷后正面朝上概率分别为
P
1
,
P
2
。
P_1,P_2。
P1,P2。
为了估计这两个概率,做实验,每次取一枚硬币,连掷5下,记录下结果,
如下(表一(实验数据)):
表一(实验数据) | 表二(测试数据) | |||||
硬币 | 结果 | 统计 | 硬币 | 结果 | 统计 | |
1 | 正正反正反 | 3正-2反 | UnKnow | 正正反正反 | 3正-2反 | |
2 | 反反正正反 | 2正-3反 | UnKnow | 反反正正反 | 2正-3反 | |
1 | 正反反反反 | 1正-4反 | UnKnow | 正反反反反 | 1正-4反 | |
2 | 正反反正正 | 3正-2反 | UnKnow | 正反反正正 | 3正-2反 | |
1 | 反正正反反 | 2正-3反 | UnKnow | 反正正反反 | 2正-3反 |
表一:可以很容易地估计出P1和P2,如下:
P1 = (3+1+2)/ 15 = 0.4
P2= (2+3)/10 = 0.5
● 现在我们抹去每轮投掷时使用的硬币标记,如上右表(表二(测试数据))。
-
目标没变,还是估计P1和P2:
-
此时多了一个隐变量z,可以把它认为是一个5维的向量(z1,z2,z3,z4,z5),
代表每次投掷时所使用的硬币,比如z1,就代表第一轮投掷时使用的硬币是1还是2。
但是,这个变量z不知道,就无法去估计P1和P2。我们可以先随机初始化一个P1和P2,用它来估计z,然后基于z,还是按照最大似然概率法则去估计新的P1和P2
当与我们初始化的P1和P2一样时,说明是P1和P2很有可能就是真实的值。这里面包含了两个交互的最大似然估计。
概率计算
-
当假设P1 = 0.2,P2 = 0.8时: (看看第一轮抛掷最可能是哪个硬币。)
-
如果是硬币1,得出3正2反的概率为:
0.2×0.2×0.2×0.8×0.8= 0.00512
如果是硬币2,得出3正2反的概率为:
0.8×0.8×0.8×0.2×0.2=0.03087
依次求出其他4轮中的相应概率。(如:表三)
按照最大似然法则:(表三) | 轮数 | 若是硬币1 | 若是硬币2 |
第1轮中最有可能的是硬币2 | 1 | 0.00512 | 0.03087 |
第2轮中最有可能的是硬币1 | 2 | 0.02048 | 0.01323 |
第3轮中最有可能的是硬币1 | 3 | 0.08192 | 0.00567 |
第4轮中最有可能的是硬币2 | 4 | 0.00512 | 0.02048 |
第5轮中最有可能的是硬币1 | 5 | 0.02048 | 0.01323 |
-
我们就把上面的值作为z的估计值。
-
然后按照最大似然概率法则来估计新的P1和P2。(由表三,查到 表二的正反数据计算P值)
P1 = (2+1+2)/15 = 0.33
P2=(3+3)/10 = 0.6
将算的的P1P2反复迭代:
设想我们知道每轮抛掷时的硬币就是如标示的那样,那么,P1和P2的最大似然估计就是0.4和0.5
(下文中将这两个值称为P1和P2的真实值)。那么对比下我们初始化的P1和P2和新估计出的P1和P2:
初始化P1 | 估计出P1 | 真实的P1 | 初始化P2 | 估计出P2 | 真实的P2 | |
0.2 | 0.33 | 0.4 | 0.7 | 0.6 | 0.5 |
3.3 示例优化:
新估计出的 P1和P2 一定会更接近真实的 P1和P2
迭代不一定会收敛到真实的P1和P2。取决于P1和P2的初始化值。
我们是用最大似然概率法则估计出的z值,然后再用z值按照最大似然概率法则估计新的P1和P2。
也就是说,我们使用了一个最可能的z值,而不是所有可能的z值。
所有的z值有多少个呢?显然,有2^5=32种,需要我们进行32次估值???
用期望来简化运算
-
再次利用表三:我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。
-
比如第1轮,使用硬币1的概率是:0.00512/(0.00512+0.03087)=0.14
使用硬币2的概率是:1-0.14=0.86我们按照期望最大似然概率的法则来估计新的P1和P2:
以P1估计为例,(表二)第1轮的3正2反相当于
0.14*3=0.42正
0.14*2=0.28反
依次算出其他四轮,表四如下:new_P1=4.22/(4.22+7.98)=0.35
轮数 | 正面 | 反面 |
---|---|---|
1 | 0.42 | 0.28 |
2 | 1.22 | 1.83 |
3 | 0.94 | 3.76 |
4 | 0.42 | 0.28 |
5 | 1.22 | 1.83 |
总计 | 4.22 | 7.98 |
可以看到,改变了z值的估计方法后,新估计出的P1要更加接近0.4。原因就是我们使用了所有抛掷的数据,而不是之前只使用了部分的数据。我们根据E步中求出的z的概率分布,依据最大似然概率法则去估计P1和P2,被称作M步。
总结