拉普拉斯近似(Laplace Approximation)
对于后验分布,往往是难以解析表达的概率分布形式,进而其性质难以描述,例如期望的计算需要积分,而后验分布常常是难以积分的。
对此,一个简单的想法是利用某种性质良好的概率分布去近似后验分布。利用高斯分布近似后验分布的方法就是拉普拉斯近似。
公式
后验分布来源于贝叶斯公式。
p
(
θ
∣
y
)
=
p
(
y
∣
θ
)
p
(
θ
)
p
(
y
)
p(\theta|y)=\frac{p(y|\theta)p(\theta)}{p(y)}
p(θ∣y)=p(y)p(y∣θ)p(θ)
在此公式中,似然函数
p
(
y
∣
θ
)
p(y|\theta)
p(y∣θ)通常是已知的,先验分布
p
(
θ
)
p(\theta)
p(θ)也是已知的,而分母
p
(
y
)
=
∫
θ
p
(
y
∣
θ
)
p
(
θ
)
d
θ
p(y)=\int_{\theta}p(y|\theta)p(\theta)d\theta
p(y)=∫θp(y∣θ)p(θ)dθ通常是难以计算的。其原因在于需要在
θ
\theta
θ上进行积分,而两个密度函数的乘积会非常难以积分,特别是
θ
\theta
θ常常是多维的。
另一方面,公式左侧是后验分布,是参数
θ
\theta
θ的函数,与变量
y
y
y无关。因此,不妨将公式右侧的分子部分记作:
p
^
(
z
)
=
p
(
y
∣
θ
)
p
(
θ
)
\hat{p}(z)=p(y|\theta)p(\theta)
p^(z)=p(y∣θ)p(θ)
省略
y
y
y,进而将归一化的后验分布记作:
p
(
z
)
=
p
(
z
)
^
Z
p(z)=\frac{\hat{p(z)}}{Z}
p(z)=Zp(z)^其中,
Z
Z
Z是归一化因子。注意这里之所以将
y
y
y省略是因为
y
y
y是某一固定取值,而非变量。这里将后验概率写为
p
(
z
)
p(z)
p(z)是为了和先验分布相区分。
在这样的设定下,拉普拉斯近似首先需要通过某种优化或其他方法,进行最大后验估计,从而得到MAP估计值。
z
M
A
P
=
arg max
z
∈
Θ
p
(
z
)
z_{MAP}=\argmax_{z \in {\Theta}}p(z)
zMAP=z∈Θargmaxp(z)
然后,我们需要在MAP值处对后验概率进行二阶的泰勒展开。进行二阶泰勒展开的原因是高斯分布是二次型,因此二阶泰勒展开会自然而然的导出一个高斯分布。
l
n
p
(
z
)
^
≈
l
n
p
(
z
)
^
−
1
2
(
z
−
z
M
A
P
)
T
H
(
z
−
z
M
A
P
)
ln\hat{p(z)}\approx ln\hat{p(z)}-\frac{1}{2}(z-z_{MAP})^TH(z-z_{MAP})
lnp(z)^≈lnp(z)^−21(z−zMAP)TH(z−zMAP)其中,
H
H
H是Hessian矩阵。那么一次项没有的原因是因为
z
M
A
P
z_{MAP}
zMAP是最大后验处,雅可比是零蛋。
两边同时取指数,仍然成立,有:
p
(
z
)
^
=
p
(
z
M
A
P
)
^
e
x
p
(
−
1
2
(
z
−
z
M
A
P
)
T
H
(
z
−
z
M
A
P
)
)
\hat{p(z)}=\hat{p(z_{MAP)}}exp(-\frac{1}{2}(z-z_{MAP})^TH(z-z_{MAP}))
p(z)^=p(zMAP)^exp(−21(z−zMAP)TH(z−zMAP))
这是对未归一化概率分布的一个二次项,则高斯分布前面的常数就很容易的通过凑高斯的方法获得,协方差矩阵的逆由Hessian矩阵给出,维度已知,则常数可以写出。
p
(
z
)
=
1
(
2
π
)
k
2
∣
H
∣
−
1
2
e
x
p
(
−
1
2
(
z
−
z
M
A
P
)
T
H
(
z
−
z
M
A
P
)
)
p(z)=\frac{1}{(2\pi)^{\frac{k}{2}}|H|^{-\frac{1}{2}}}exp(-\frac{1}{2}(z-z_{MAP})^TH(z-z_{MAP}))
p(z)=(2π)2k∣H∣−211exp(−21(z−zMAP)TH(z−zMAP))
上式就是完整的后验分布的拉普拉斯近似。
证据项的估计
基于对后验分布的拉普拉斯近似,对于贝叶斯公式上难以求得的分母也有了表达:
∫
θ
p
(
θ
∣
y
)
p
(
θ
)
d
θ
≈
∫
z
p
(
z
M
A
P
)
^
e
x
p
(
−
1
2
(
z
−
z
M
A
P
)
T
H
(
z
−
z
M
A
P
)
)
d
z
\int_{\theta}p(\theta|y)p(\theta)d\theta\\\approx \int_{z}\hat{p(z_{MAP)}}exp(-\frac{1}{2}(z-z_{MAP})^TH(z-z_{MAP}))dz
∫θp(θ∣y)p(θ)dθ≈∫zp(zMAP)^exp(−21(z−zMAP)TH(z−zMAP))dz
注意,
θ
\theta
θ和
z
z
z是一个变量,只不过上面为了区分后验和先验分布的表达,写作了两个变量。
把前面的
p
(
z
M
A
P
^
)
\hat{p(z_{MAP}})
p(zMAP^)提出来,然后积里面的高斯分布,实际上就会得到高斯分布前面归一化常数的倒数,也就是
(
2
π
)
k
2
∣
H
∣
−
1
2
(2\pi)^{\frac{k}{2}}|H|^{-\frac{1}{2}}
(2π)2k∣H∣−21这个东西,因此,将积分结果再乘以提出去的
p
(
z
M
A
P
^
)
\hat{p(z_{MAP}})
p(zMAP^),就得到:
(
2
π
)
k
2
∣
H
∣
−
1
2
p
(
z
M
A
P
^
)
(2\pi)^{\frac{k}{2}}|H|^{-\frac{1}{2}}\hat{p(z_{MAP}})
(2π)2k∣H∣−21p(zMAP^)
把
p
(
z
M
A
P
^
)
\hat{p(z_{MAP}})
p(zMAP^)完整写出,在整理,就是:
p
(
y
)
=
(
2
π
)
k
2
∣
H
∣
−
1
2
p
(
y
∣
θ
M
A
P
)
p
(
θ
M
A
P
)
p(y)=(2\pi)^{\frac{k}{2}}|H|^{-\frac{1}{2}}p(y|\theta_{MAP})p(\theta_{MAP})
p(y)=(2π)2k∣H∣−21p(y∣θMAP)p(θMAP)
这在贝叶斯实验设计中可用作采样方法的一个近似。(Ryan, 2003)
Kenneth J Ryan (2003) Estimating Expected Information Gains for Experimental Designs With Application to the Random Fatigue-Limit Model, Journal of Computational and Graphical Statistics, 12:3, 585-603.