Saddlepoint Approximation

和Edgeworth展开相比，鞍点近似是更加精细的分布近似方法，鞍点近似同时适用于连续和离散的分布。在许多极端的条件下，鞍点近似都表现出惊人的精度。鞍点近似有许多角度的理解，然而下面的推导展示了鞍点近似和Edgeworth的关系，以及它为什么具有高精度。

Tilted Exponential Families

设 $f(x)$ 是一个密度函数，在区间 $(a,b)$ 具有矩母函数 $M(s) = E e^{sX} $ 和累量母函数 $ K(s) = \log M(s) $。则下列指数分布族称为 tilted 正则指数族 $$ f(x;s) = \exp \{ sx - K(s) \} f(x), \qquad s \in (a,b) $$
$f(x;s)$ 称作 s-tilted 密度，用 $X_s$ 表示服从此密度的一个随机变量。容易推导 $X_s$ 的 均值和方差满足
$$ E(X_s) = K’(s)$$ $$ Var(X_s) = K’’(s)$$
以及 $X_s$ 的高阶标准化累量
$$ \kappa_i(s) = {K^{(i)}(s) \over \{ K’’(s) \}^{i/2} }, \qquad i \geq 3 $$

Saddlepoint Approximation

由 tilted 指数族定义可知
$$ f(x) \equiv \exp \{ K(s) -sx \} f(x;s), \qquad \forall s \in (a,b) $$
注意这里是对任意的 s，两边都是一样的密度函数。所以我们可以针对不同的 x，选择不同的 s，使近似效果达到最佳。

考虑 $ f(x;s) $ 的Edgeworth 展开
$$ f(x;s) = \phi \left( {x -\mu_s \over \sigma_s} \right) \left\{ 1 + {\kappa_{3}(s) \over 6 } \left[ \left( {x -\mu_s \over \sigma_s} \right)^3 - 3 \left( {x -\mu_s \over \sigma_s} \right) \right] + O \left( {1 \over n} \right) \right\} $$
对于固定的 x，若选择s，使得 $\mu_s =x$，则在x点上的近似具有 $O(1/n)$ 的精度。现在对不同的x，我们始终选择 $\hat s$ 满足
$$ K’(\hat s(x)) = x $$
从而可以使 $O(1/n)$ 的精度是全局的。上面的方程被称作鞍点方程。将 $ f(x;s) $ 的近似带入到 $ f(x) $ 中
$$ f(x) \approx \hat f(x) = {1 \over \sqrt{2 \pi K’’(\hat s (x)) }} \exp \{ K(\hat s(x)) -\hat s(x)x \} $$
上式称为 $ f(x) $ 的（一阶）鞍点近似。要注意的是 $ \hat f(x) $ 并不是真正的密度函数，因为
$$ c := \int_{\Bbb R} \hat f(x) dx \neq 1 $$
正规化的鞍点近似密度函数为
$$ \bar f(x) = c^{-1} \hat f(x), \qquad x \in \mathcal X $$

连续分布的鞍点近似

设 $X_i, i \in \Bbb{N}$ 是独立同分布的连续随机变量，具有累量母函数 $K(s)$。则均值 $ \bar X =1/n \sum_{i=1}^n X_i$ 密度函数的鞍点近似为
$$ \hat f _{\bar X}(x) = \sqrt{ { n \over 2 \pi K’’(\hat s (x)) } } \exp \{ K( \hat s(x)) -\hat s(x)x \} $$

$X$ 分布的函数的鞍点近似为
$$ \hat F_{\bar X}(x) = \begin{cases}
\Phi(\hat w) + \phi(\hat w) \left( {1 \over \hat w} - {1 \over \hat u} \right) & \text{if }x \neq \mu \\
{1 \over 2} + { K’’’(0) \over 6 \sqrt{2} \pi K’’(0)^{3/2} } & \text{if }x = \mu
\end{cases} $$
其中
$$\begin{align}
\hat w &= \text{sgn}(\hat s) \sqrt{ 2 \{ \hat s x -K(\hat s) \} } \\
\hat u &= \hat s \sqrt{ K’’(\hat s) }
\end{align} $$
$\Phi(\cdot),\phi(\cdot) $ 分别为标准正态分布的分布函数和密度函数，$\hat s$ 是鞍点方程 $K’(s) =x $ 的解。

离散分布的鞍点近似

设 $X$ 是整数集 $\Bbb{Z}$ 上的随机变量， $p(k)$ 是其质量函数，则 $p(k) $ 具有鞍点近似
$$ \hat p(k) = {1 \over \sqrt{ 2 \pi K’’(\hat s) }} \exp \{ K(\hat s) - \hat s k \} $$
其中 $\hat s$ 是鞍点方程 $K’(\hat s) = k, \quad k \in \mathcal{I_X} $的解， $\mathcal{I_X} $ 是 supp X 延展的内点集。

定义 $ k^- =k-1/2 \in \mathcal{I_X} $ 是 k的偏置，$\hat s$ 是鞍点方程 $K’(\hat s) = k^-$ 的解，则$X$ 生存函数的（第二连续性矫正）鞍点近似为
$$ \hat{Pr}(X \geq k) = \begin{cases}
1 - \Phi(\hat w) - \phi(w)(1/2 - 1/u) & \text{if } k^- \neq \mu \\
1/2 - { K’’’(0) \over 6 \sqrt{2} K’’(0)^{3/2} } & \text{if } k^- = \mu
\end{cases} $$
其中
$$\begin{align}
\hat w &= \text{sgn}(\hat s) \sqrt{2 \{ \hat s k^- -K(\hat s) \}} \\
\hat u &= 2 \text{sinh}(\hat s/2) \sqrt{K’’(\hat s)}
\end{align}$$

Sum of Binomials

设 $X_i,i=1,2,3$ 是独立随机变量，$X_i \sim \text{Binomial}(m_i,\theta_i) $，其中 $ m_i = 8-2i, \theta_i = 2^i/30 $。现在考虑 $ X = \sum_i X_i $ 的质量函数和分布函数。相信这样的事情能让几乎每一个人抓狂。

$X$ 的CGF为
$$ K(s) = 2 \sum_{i=1}^3 (4-i) \log \left\{ {2^i \over 30} (e^s -1) +1\right \}, \quad s \in \Bbb{R} $$
对应的鞍点方程为
$$ K’(s) = {6\hat t \over \hat t +14} + {8\hat t \over 2\hat t +13} +{8\hat t \over 4\hat t +11} $$
其中 $ \hat t = \exp(\hat s) $。下表展示了上述的鞍点近似质量函数 $\hat p(k)$，和对应的归一化鞍点近似 $ \bar p(k) $

在k=1,3时，Poisson近似差强人意，更大的k时，Poisson近似就稀烂了。而（对这个问题）鞍点近似在全局都保持着惊人的精度，注意 $ 0 \leq X \leq 12 $。通常情况下，对密度进行归一化能得到更好的近似，但是归一化常数的计算可能并不直接。

另一方面，对于分布函数，下表展示了鞍点近似的生存函数（右尾概率）

Gumbel(0,1) Distribution

设 $ X \sim \text{Gumbel}(0,1) $，具有分布函数 CDF
$$ F(x) = e^{- e^{-x}}, \qquad x \in (-\infty, \infty) $$
和累量母函数 CGF
$$ K(s) = \log \Gamma (1-s), \qquad s \in (-\infty, 1) $$
鞍点方程为
$$ -\digamma(1- \hat s) = x, \qquad s < 1 $$
其中 $ \digamma(\cdot) $ 是digamma函数 (gamma 函数对数的一阶导数，$ \digamma(z) = { \text d \over \text{dz} } \log \Gamma(z) = {\Gamma’(z) \over \Gamma(z)} $)。
于是 $f(x)$ 鞍点近似为
$$ \hat f(x) = {1 \over \sqrt{2 \pi K’’(\hat s (x)) }} \exp \{ K(\hat s(x)) -\hat s(x)x \} $$
正规化的鞍点近似为
$$ \bar f(x) = c^{-1} \hat f(x) $$
其中
$$ c = \int_{\Bbb R} \hat f(x) dx \approx 0.9792 $$

另一方面，$ EX = \gamma := \lim_{n \rightarrow \infty} (-\log n + \sum_{k=1}^n 1/k) \approx 0.5772 $， $\gamma$ 是 Euler–Mascheroni 常数， $ Var X = \pi^2/6$，高阶标准化累量满足
$$ \kappa_n = { K^{(n)}(0) \over \sigma^n } = {1 \over \sigma^n} (n-1)! \zeta(n) $$
其中 $\zeta (s) = \sum_{n=1}^{\infty} n^{-s} $ 是 Riemann zeta 函数。特别地 $\kappa_3 \approx 1.14, \kappa_4 \approx 2.40 $。$X$ 的Edgeworth近似为
$$ f(x) = \phi(z) \left[ 1+ {\kappa_3 \over 6 } H_3(z) +
\left\{ {\kappa_4 \over 24}H_4(z) + {\kappa_3^2 \over 72}H_6(z) \right\} + O(n^{-3/2}) \right] $$
其中 $\phi(\cdot)$ 是标准正态分布的密度函数，$z=(x-\mu)/\sigma $，$H_k(\cdot)$ 是 Hermite 多项式。

下图展示了鞍点近似和 Edgeworth 近似，鞍点近似的效果非常惊人，正规化的鞍点近似更是（视觉上）与真实的密度无法区分。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

g = function(s, x) {
  -digamma(1 - s) - x
}

# root of saddlepoint equation
saddleEq = function(x) {
  sapply(x, function(t)
    uniroot(g, c(-1e6, 0.9999), tol = 1e-6, x = t)$root)
}

# saddlepoint density
f_saddle = function(x) {
  s = saddleEq(x)
  1 / sqrt(2 * pi * trigamma(1 - s)) * exp(lgamma(1 - s) - s * x)
}

c = integrate(f_saddle,-10, 10)$value

# normalized saddlepoint density
f_norm_saddle = function(x) {
  f_saddle(x) / c
}

# Hermite polynomials
H3 = function(x) (x ^ 3 - 3 * x)
H4 = function(x) (x ^ 4 - 6 * x ^ 2 + 3)
H6 = function(x) (x ^ 6 - 15 * x ^ 4 + 45 * x ^ 2 - 15)

mu = 0.577215665
sigma = sqrt(pi * pi / 6)
k3 = factorial(2) * 1.202057 / sigma ^ 3
k4 = factorial(3) * 1.082323 / sigma ^ 4

# first-order Edgeworth
f_Edge1 = function(x) {
  z = (x - mu) / sigma
  dnorm(z) * (1 + k3 * H3(z) / 6) / sigma
}

# second-order Edgeworth
f_Edge2 = function(x) {
  z = (x - mu) / sigma
  f_Edge1(x) + (k4 * H4(z) / 24 + k3 ^ 2 * H6(z) / 72) * dnorm(z) / sigma
}

# density of Gumbel(0,1)
f_Gumbel = function(x) {
  exp(-x - exp(-x))
}

label = c(
  'Exact',
  'Saddlepoint',
  'Normalized Saddlepoint',
  'First-order Edgeworth',
  'Second-order Edgeworth'
)
x = seq(-3, 5, 1e-2)
data = data.frame(
  density = c(
    f_Gumbel(x),
    f_saddle(x),
    f_norm_saddle(x),
    f_Edge1(x),
    f_Edge2(x)
  ),
  x = rep(x, length(label)),
  Approximations = rep(label, each = length(x))
)

ggplot(data, aes(x, density, color = Approximations)) +
  scale_colour_manual(
    breaks = label,
    values = c('#000000', '#e41a1c', '#377eb8', '#4daf4a', '#984ea3')
  ) +
  theme(legend.position = c(0.85, 0.75)) +
  geom_line(alpha = 0.6)

Edgeworth Expansion

许多时候，我们感兴趣的随机变量 $X$ 具有复杂的结构，我们无法或者很难得到它的精确分布密度PDF或分布函数CDF，对于 $X$ 的研究比如求概率 $P(X \in A)$ 就变得不那么直接。多数情况下，$X$ 可以表示为一些独立甚至同分布随机变量的和（均值）。对于平凡的问题，我们可以利用中心极限定理CLT将 $X$ 的分布近似为正态分布。然而对于一些问题比如

• $\bar X={1 \over n} \sum_{i=1}^n X_i$，n很小只有几或者几十
• $P(\bar X \in A)$ 很小比如 $P(\sqrt{n}|\bar X -\mu| > 4\sigma)$

此时正态分布的近似就很差。而Edgeworth展开作为中心极限定理的自然延伸，从一定程度上改善了近似的效果

Tchebycheff–Hermite Polynomials

设 $\phi(z)$ 是标准正态分布 $N(0,1)$ 的密度函数，由归纳法容易证明，存在多项式 $ \{H(z) \}_k$ 成立
$${ d^k \over dz^k} \phi (z) = (-1)^k H_k(z) \phi (z) $$
其中 $\{ H_k(z): k \in \Bbb{N} \}$ 称为Tchebycheff–Hermite多项式。对上式求导，得到
$${d \over dx} [ H_k(z)\phi(z) ]=-H_{k+1}(z)\phi(z) $$
由此可以递推出 $H_k(z)$，前几项为
$$\begin{align}
H_0(z) &=1 & H_3(z)&= z^3 − 3z \\
H_1 (z) &= z & H_4 (z) &= z^4 − 6z^2 + 3 \\
H_2 (z) &= z^2 − 1 & H_5 (z) &= z^5 − 10z^3 + 15z
\end{align}$$

$\{ H_k(z): k \in \Bbb{N} \}$ 是正交多项式，在如下的加权內积意义下
$$<H_i(z), H_j(z)> = \int_{\Bbb{R} }H_i (z) H_j (z)\phi (z)dz=
\begin{cases}
0 & \text{if }i\neq j \\
i! & \text{if } i=j
\end{cases}$$

Inverse Fourier Transform

设 $X \sim g(x)$ 具有特征函数 $h(u)$，则 $h(u)$ 具有逆变换
$$g(x)={1 \over 2\pi } \int_{\Bbb{R}} e^{-iux} h(u)du$$
特别地，由Fourier变换的性质，成立
$$\begin{align}
{1 \over 2\pi } \int_{\Bbb{R}} e^{-iux} e^{-u^2/2} (iu)^k du
&= {(-1)^k \over 2 \pi} {d^k \over dx^k} \int_{\Bbb{R}} e^{-iux} e^{-u^2/2}du \\
&= (-1)^k {d^k \over dx^k} \phi(x) \\
&= H_k(x)\phi(x)
\end{align}$$
也就是说若 $X$ 具有特征函数 $e^{-u^2/2} (iu)^k$，则 $X$ 具有密度函数 $H_k(x)\phi(x)$。

Edgeworth Expansion 的推导

设 $X_i$ 有均值0方差1，有有限四阶矩，记 $\gamma=E X_j^3, \tau=E X_j^4$。现在考虑正规化均值 $S_n = \sqrt{n} \bar X$ 的特征函数。
$$h_{S_n} (u)=E \exp \{ ( iu/ \sqrt{n} ) \sum X_j \} =[ h_X(u/\sqrt{n}) ]^n$$
对 $\exp \{ iuX/\sqrt{n} \}$ 作泰勒展开得到
$$h_X \left( {u \over \sqrt{n}} \right) = E e^ { iu {X \over \sqrt{n}} }
= \left( 1- {u^2 \over 2n} \right) +
{(iu)^3 \gamma \over 6n \sqrt{n} } +
{(iu)^4 \tau \over 24n^2 } +
o \left( {1 \over n^2} \right)$$
两边取n次幂，再由二项展开式和
$$ \left( 1+ {a \over n} \right) ^{n-k} =
e^a \left( 1- {a(a+k) \over 2n} \right) +
o \left( {1 \over n} \right)$$
得到
$$h_{S_n}(u) = e^{-u^2/2} \left[ 1 +
{(iu)^3 \gamma \over 6 \sqrt{n} } +
{(iu)^4 ( \tau -3) \over 24n } +
{(iu)^6 \gamma^2 \over 72n} \right] +
o \left( {1 \over n} \right)$$
最后由上述关于特征函数和密度函数的逆变换，得到 $S_n$ 的密度函数为
$$\begin{align}
g_{S_n}(x) &= {1 \over 2\pi } \int_{\Bbb{R}} e^{-iux} h_{S_n}(u) du \\
&= \phi(x) \left( 1+
{1 \over \sqrt{n}} {\gamma H_3(x) \over 6} +
{1 \over n} \left[ {(\tau-3) H_4(x) \over 24} +
{\gamma^2 H_6(x) \over 72} \right] \right) +
o \left( {1 \over n} \right)
\end{align}$$
两边求积分，利用上面 $H_k(x) \phi(x)$ 的递推关系得到 $S_n$ 分布函数的近似
$$G_{S_n}(x)=
\Phi(x) - \phi(x) \left(
{1 \over \sqrt{n}} {\gamma H_2(x) \over 6} +
{1 \over n} \left[ {(\tau-3) H_3(x) \over 24} +
{\gamma^2 H_5(x) \over 72} \right] \right) +
o \left( {1 \over n} \right)
$$
包含上述 $1/\sqrt{n}$ 和 $ 1/n$ 的近似，分别叫做Edgeworth一阶/二阶展开。

此外值得一提的是

• 若 $X_i$ 具有对称分布，三阶矩 $\gamma=0$，通常的正态分布近似已经具有一阶Edgeworth近似的精度，
• $g_{S_n}(x)$ 在0点（非正规化的随机变量在均值 $\mu$ 点），由于奇数Hermite多项式在0点为0，即 $H_3(0)=0$，故一阶Edgeworth近似具有 $O(1/n)$ 的精度，
• Edgeworth展开并未进行概率意义上的规约，所以计算结果可能超出 $[0,1]$ 之外。
• 高阶的Edgeworth展开虽然具有理论上的高阶渐进精度，但也容易矫枉过正，使结果不可靠，因而仅具有理论意义。

指数分布的和

设 $X_i \sim \text{Exp}(1),i=1,…,5$ 则 $\mu = \sigma^2=1, E(X_i-\mu)^3/\sigma^3 = 2, E(X_i-\mu)^4/\sigma^4 = 6$。现在考虑 $X_i$ 的正规化和 $S_n = \sqrt{n} (\bar X -\mu)$ 的分布，由于 $\sum X_i \sim \Gamma(5,1)$，我们可以准确知道 $S_n$ 的分布。下面的代码和图1展示了正态分布和Edgeworth近似的效果。

n = 5
b1 = 2
b2 = 6
H3 = function(x) (x ^ 3 - 3 * x)
H4 = function(x) (x ^ 4 - 6 * x ^ 2 + 3)
H6 = function(x) (x ^ 6 - 15 * x ^ 4 + 45 * x ^ 2 - 15)

f1 = function(x) {
  dnorm(x) * (1 + b1 * H3(x) / 6 / sqrt(n))
}

f2 = function(x) {
  f1(x) + (b2 * H4(x) / 24 / n + b1 ^ 2 * H6(x) / 72 / n) * dnorm(x)
}
x = seq(-3, 3, 0.01)
fn_true = dgamma((x + sqrt(n)) * sqrt(n), n, 1) * sqrt(n)
fn_E0 = dnorm(x)
fn_E1 = f1(x)
fn_E2 = f2(x)
label = c('Exact',
          'Normal',
          'First-order Edgeworth',
          'Second-order Edgeworth')
data = data.frame(
  density = c(fn_true, fn_E0, fn_E1, fn_E2),
  x = rep(x, 4),
  estimate = rep(label, each = length(x))
)
ggplot(data, aes(x, density, color = estimate)) +
  scale_colour_manual(breaks = label,
                      values = c('#000000', '#1b9e77',  '#d95f02', '#7570b3')) +
  theme(legend.position = c(0.85, 0.82)) +
  geom_line()

由上图可以看出

• 正态分布的近似效果极差，由于样本量太小，无法纠正指数分布的偏度
• Edgeworth近似在分布的主体部分有不错的近似效果
• 在分布的尾部，Edgeworth的效果的仍不好，甚至在左尾出现了负值

下面的表格展示了一些区间概率的估计，印证了上面的论断

Laplace Approximation

统计中的计算问题大致可以分成两类，积分问题和最优化问题，积分问题常常和Bayes方法有关，最优化问题常常和似然方法相关。求积分的方法，除了基于模拟的Monte Carlo方法外，也有一些基于分析的方法在统计中常见，其中最古老也最有用的方法就是Laplace近似。

假设我们感兴趣的积分问题是
$$\int_A f(x|\theta) dx$$
其中 $f$ 是非负可积函数，$\theta$ 是固定的参数。

将被积函数写作 $f(x|\theta)=\exp \{ nh(x|\theta) \}$，其中n是样本大小或者某个可以趋近无穷大的参数，对 $h(x|\theta)$ 在 $x_0$ 点作Taylor展开，得到
$$h(x|\theta) =
h(x_0|\theta)+(x-x_0)h’(x_0|\theta)+{(x-x_0)^2 \over 2!} h’’(x_0|\theta)+
{(x-x_0)^3 \over 3!} h’’’(x_0|\theta) + R_n(x)$$
其中 $R_n(x) = o((x-x_0)^3), x \to x_0.$

现在选择 $x_0=\hat x_{\theta}$，满足 $h’(x_0|\theta)=0$ 且最大化 $h(x|\theta)$。在 $\hat x_{\theta}$ 的某个邻域成立
$$\int_A e^{nh(x|\theta)} dx \approx e^{nh(\hat x_{\theta}|\theta)} \int_A
e^{n {(x-x_0)^2 \over 2!} h’’(\hat x_{\theta}|\theta)}
e^{n {(x-x_0)^3 \over 3!} h’’’(\hat x_{\theta}|\theta)}$$
对于三阶导项，若考虑指数展开 $e^y \approx 1+y+y^2/{2!}$ ，得到近似
$$e^{n {(x-x_0)^3 \over 3!} h’’’(\hat x_{\theta}|\theta)} \approx
1+ n {(x-x_0)^3 \over 3!} h’’’(\hat x_{\theta}|\theta)+
n^2 {(x-x_0)^6 \over {2! (3!)^2 }} [ h’’’(\hat x_{\theta}|\theta) ]^2$$

带回去，我们得到
$$\int_A e^{nh(x|\theta)} dx \approx
e^{nh(\hat x_{\theta}|\theta)} \int_A
e^{n {(x-x_0)^2 \over 2!} h’’(\hat x_{\theta}|\theta)} \times
\left[ 1+ n {(x-x_0)^3 \over 3!} h’’’(\hat x_{\theta}|\theta)+
n^2 {(x-x_0)^6 \over {2! (3!)^2 }} [ h’’’(\hat x_{\theta}|\theta) ]^2 \right] dx$$
上式方括号内分别取前1/2/3项，分别叫做一阶/二阶/三阶Laplace近似。

特别地，对于一阶近似，注意到其和正态分布密度的关系，我们可以将其表示为
$$\int_A e^{nh(x|\theta)} dx \approx
e^{nh(\hat x_{\theta}|\theta)}
\sqrt{2 \pi \over {-n h’’(\hat x_{\theta}|\theta)}} \times
\left \{ \Phi [ \sqrt{-n h’’(\hat x_{\theta}|\theta)} (b- \hat x_{\theta}) ] -
\Phi [ \sqrt{-n h’’(\hat x_{\theta}|\theta)} (a- \hat x_{\theta}) ]\right \} $$
其中，取 $A=[ a,b ]$，$\Phi(\cdot)$ 为标准正态分布的cdf。

从上面的形式可以看出，Laplace近似的计算不存在本质困难，但二三阶近似的计算略显繁琐。

Gamma积分的近似

若以一例来说明Laplace近似的功效，我们考虑Gamma分布 $\mathcal{G}a(\alpha, 1/ \beta)$ 的积分
$$\int_a^b {x^{\alpha-1} \over {\Gamma (\alpha) \beta^\alpha}} e^{-x/ \beta} dx$$

略去分母的常数，首先求得 $h(x) = - {x \over \beta} +(\alpha-1) \log(x)$，一二三阶导数分别为
$$\begin{align}
h’(x) &= {\alpha -1 \over x}- {1 \over \beta} \\
h’’(x) &= - {\alpha-1 \over x^2} \\
h’’’(x) &= {2(\alpha-1) \over x^3}
\end{align}$$
对于 $\alpha=5, \beta=2$，选取 $\hat x_\theta = (\alpha-1)\beta=8$。对于不同的积分区间，计算Laplace近似，得到下表

由此可以看出，对于一般的包含的 $\hat x_\theta$ 的区间，Laplace近似效果不错，尤其的关于 $\hat x_\theta$ 对称的区间，由于二阶近似等于一阶近似，一阶近似的值已经相当精确。但是对于偏离 $\hat x_\theta$ 的区间，近似的偏差较大，对于尾概率的估计更是难以接受。相较于费时的Monte Carlo方法，快捷的Laplace近似可以看作最终解决方案的一个指导。

ARS (Adaptive Rejection Sampling)

ARS (Adaptive Rejection Sampling) 是一种生成具有对数凹密度 (log-concave density)的随机数的方法，ARS使用包围的接受拒绝 (Envelope Accept-Reject) 方法，而Envelope接受拒绝方法是一般接受拒绝方法的一个变形（改进）。

Envelope Accept-Reject

和普通的接受拒绝方法相比，EAR多了一个下界，即，若存在密度$g_m$，函数$g_l$和正常数$M$满足
$$g_l(x) \leq f(x) \leq Mg_m(x),$$
则下列算法产生随机变量服从分布$f$

1. 生成 $X \sim g_m(x), U \sim \mathcal{U}_{[0,1]}$
2. IF $U \leq g_l(X)/Mg_m(X)$
3. THEN 接受$X$
4. ELSEIF $U \leq f(X)/Mg_m(X)$
5. THEN 接受$X$

这个算法的成立是显然的，若有$U \leq g_l(X)/Mg_m(X)$，必有$U \leq f(X)/Mg_m(X)$，所以和普通的接受拒绝方法相比，EAR可能在第3行就结束，而不用计算$f(x)$的值，这个优点当$f(x)$很难计算时便很有用。

当然EAR并不能完全避免计算$f(x)$的值，具体来讲，节省的比例为
$$\begin{align}
P(U \leq g_l(X)/Mg_m(x)) &= E[P(U \leq g_l(X)/Mg_m(X) | X)] \\
&= E[g_l(X)/Mg_m(X)] \\
&= \frac{1}{M} \int{g_m(x)dx}
\end{align}$$

值得一提的是，上述的函数$g_l, f, g_m$，都可以之多差一个系数，但是这样的话，常数$M$不再具有上述的含义。以正态分布$N(0,1)$为例，略去常数，取 $f(x)=\exp (-x^2/2)$。由泰勒公式，
$$e^{-x^2/2} \geq 1-\frac{x^2}{2}$$
于是可取
$$g_l(x) = \left (1-\frac{x^2}{2} \right ) \lor 0$$
另一方面若取$g_m(x)$为Laplace(1)分布，即$g_m(x)=\lambda /2 \exp (-\lambda x), \lambda=1$，则容易计算接受率为76%，引入Envelope下界，可避免58%的$f(x)$的计算，对于复杂的目标密度，这个节省还是很可观。

Adaptive Rejection Sampling

自适应拒绝采样(ARS)的想法很简单，基于上述EAR的思路，对于具有对数凹密度的随机变量，由凹函数的几何性质，很容易构造出自适应的Envelope上下界，准确来讲是构造分段线性函数的上下界。

如上图所示，设 $h=\log f$ 凹，$S_n = \{x_i, i=0,1,…,n+1\}$ 是一列从小到大的 $f$ 支撑上的点。由 $h$ 的凹性，直线 $L_{i,i+1}$ 穿过点 $(x_i,h(x_i))$ 和点 $(x_{i+1},h(x_{i+1}))$，在区间 $[x_i,x_{i+1}]$ 上在 $h$ 的下方，在这个区间外在 $h$ 的上方，于是我们定义
$$ \overline{h}_n (x) = \min \{ L_{i-1,i}(x), L_{i+1,i+2}(x) \}$$

如上图红线，定义
$$\underline{h}_n(x) = L_{i,i+1}(x)$$

如上图蓝线所示。在集合 $[ x_0,x_{n+1} ]^c$ 上补充
$$\underline{h}_n(x) = - \infty, \quad \overline{h}_n (x) = \min \{ L_{0,1}(x), L_{n,n+1}(x)$$
这样，在 $f$ 的支撑上成立
$$\underline{h}_n(x) \leq h(x) \leq \overline{h}_n (x)$$
最后取指数，即 $\underline{f}_n(x) = \exp \underline{h}_n(x)$，$\overline{f}_n(x) = \exp \overline{h}_n(x)$，有
$$\underline{f}_n(x) \leq f(x) \leq \overline{f}_n (x)= \varpi g_n(x)$$
最右边的等号表示，将 $\overline{f}_n(x)$ 归一化为密度函数，$\varpi$是正规化常数。

这样我们得到了Envelope的上下界，带入到EAR算法中，我们得到如下的ARS算法

0. 初始化 $n$ 和 $S_n$
1. 生成 $X \sim g_n(x), U \sim \mathcal{U}_{[ 0,1 ]}$
2. IF $U \leq \underline{f}_n(X) / \varpi g_n(X)$
3. THEN 接受 $X$
4. ELSEIF $U \leq f(X) / \varpi g_n(X)$
5. THEN 接受 $X$
6. $\qquad$更新 $S_n \to S_{n+1}=S_n \cup \{ X \}$

ARS算法在生成随机数过程中会不断更新 $S_n$，使得Envelope上下界越来越接近 $f$，算法的效率越来越高，体现在接受率越来越接近1,和计算 $f(x)$ 的频率越来越少（趋近于0）。

另外要指出的是，算法中的常数 $\varpi$ 并不用计算，因为 $\varpi$ 总是和 $g_n(X)$ 一起，而 $\varpi g_n(X) = \overline{f}_n (x)$。

到此为止，ARS算法只差的步骤1的生成 $X \sim g_n(X)$，它被解决如下

Generate $X \sim g_n(x)$

$g_n(x)$ 是分段函数，我们将它重新写成
$$g_n(x) = \varpi^{-1} \left \{ \sum_{i=0}^{r_n} e^{\alpha_i x + \beta_i} \Bbb{I}_{[ x_i,x_{i+1} ]}(x) +
e^{\alpha_{-1} x + \beta_{-1}} \Bbb{I}_{[- \infty,x_0]}(x) +
e^{\alpha_{r_n+1} x + \beta_{r_n+1}} \Bbb{I}_{[x_n+1,+ \infty]}(x) \right \}$$

其中 $y=\alpha_i x + \beta_i$ 对应区间 $[ x_i,x_i+1 ]$，这其中包括相邻两个线段的交点，直线 $y=\alpha_i x + \beta_i$ 和 $y=\alpha_{i+1} x + \beta_{i+1}$ 的交点在 $x=-(\beta_{i+1}-\beta_i)/(\alpha_{i+1}-\alpha_i)$ 处。$y=\alpha_{-1} x + \beta_{-1}$ 和 $y=\alpha_{r_n+1} x + \beta_{r_n+1}$ 分别对应区间 $[- \infty, x_0]$ 和 $ [ x_{n+1}, + \infty ]$， $r_n$ 表示线段的数量

由 $\varpi g_n(X) = \overline{f}_n (x)$，两边积分可得
$$\begin{align} \varpi_n &= \int_{- \infty}^{x_0} e^{\alpha_{-1} x + \beta_{-1}} dx +
\sum_{i=0}^{r_n} \int_{x_i}^{x_{i+1}} e^{\alpha_{i} x + \beta_{i+1}} dx +
\int_{x_{n+1}}^{+ \infty} e^{\alpha_{r_n+1} x + \beta_{r_n+1}} dx \\
&= \frac{e^{\alpha_{-1} x_0 + \beta_{-1}}}{\alpha_{-1}} +
\sum_{i=0}^{r_n} e^{\beta_i} \frac{\alpha_{i} x_{i+1} - \alpha_{i} x_{i}}{\alpha_{i}} -
\frac{\alpha_{r_n+1} x_{n+1}}{\alpha_{r_n+1}}
\end{align}$$

由上面的计算知道，$g_n(x)$ 在区间 $[ x_i,x_{i+1} ]$ 上的概率
$$\omega_i = \begin{cases}
\frac{e^{\alpha_{-1} x_0 + \beta_{-1}}}{ \varpi \alpha_{-1}} &\text{if } i=-1, \\
e^{\beta_i} \frac{\alpha_{i} x_{i+1} - \alpha_{i} x_{i}}{\varpi \alpha_{i}} &\text{if } 0 \leq i \leq r_n, \\
-\frac{\alpha_{r_n+1} x_{n+1}}{\varpi \alpha_{r_n+1}} &\text{if } i = r_n+1.
\end{cases}$$

其中 $ x_{-1}=- \infty, x_{r_n+2}=+ \infty$。利用逆分布分布函数方法，注意到若 $ X \sim g_n(x)$，则 $X|x_i \leq X \leq x_{i+1} \sim h(x)=c^{-1} \exp \{ \alpha_i x +\beta_i \} $，其中 $ c=\alpha^{-1} ( e^{\alpha_i x_{i+1} +\beta_i} -e^{\alpha_i x_{i} +\beta_i} )$，我们得到如下的算法生成随机变量 $X \sim g_n(x)$

1. 依概率 $\omega_i$ 选择区间 $[ x_i,x_{i+1} ]$
2. 生成 $U \sim \mathcal{U}_{[ 0,1 ]}$，取
   $$ X = \alpha^{-1} \log [ e^{\alpha_i x_i} + U (e^{\alpha_i x_{i+1}} - e^{\alpha_i x_i}) ]$$

ARS算法的C语言实现

基于上述描述，我们用C语言实现ARS算法，由于算法中涉及许多动态动态内存的操作，考虑到结构化和性能的平衡，我们设计程序的结构如下。

首先是保存一个线段的结构 Segment，包含线段的始末横坐标，斜率，截距和取完指数后的积分。

1
2
3
4
5
6

typedef struct 
{
    float x1, x2;
    float alpha, beta;
    float area;
} Segment;

结构体Interval保存了区间 $[ x_i, x_{i+1} ]$ 上的envelope上下界，包含下界一段和上界（至多）两段线段的信息。

1
2
3
4
5
6
7

typedef struct
{
    double x1, x2;
    double alpha, beta;
    Position flag;
    Segment leftseg, rightseg;
} Interval;

最后我们使用二叉搜索树来组织这些Interval，节点的序关系由区间的前后定义。二叉树的实现使得对分段函数的查找和求值达到 $O(\log n)$ 的时间复杂度。

生成随机数的函数rars()接受一个指向求密度函数对数的函数指针(*log_density)(double) ，不同的函数指针可以包装成不同的随机数函数。

下面的显示了部分核心代码，完整的源代码参见 GitHub仓库。

最后值得一提的是ARS的算法我觉得唯一的优点就是泛用性，因为许多有用的分布都是对数凹密度。ARS算法看起来很理想，但是在实现上的缺点是它包含大量动态内存操作，这些内存操作的开销已经远超过计算本身，此外由于ARS算法构造分段的envelope上下界，随着集合 $S_n$ 不断增大，对分段函数进行求值得开销也越来越大。所以我觉得不如一开始就设定一个充分大的集合 $S_n$，比如包含几百个分割点，无需在生成随机数过程中更新 $S_n$。

An Optimization with SIMD Instructions

本文介绍利用SIMD指令（AVX2）优化运算的一点探索，以单精度float类型的加法和乘法为例。

计算模型：

• hardware: E5 2678 v3 2.5GHz; X99; DDR4 2133 MHz 16GB x1
• software: win10 17134; mingw-W64(gcc) 8.1.0

SIMD和AVX2简介

SIMD (single instrucion multiple data)单指令多数据是用来对多个数据执行单个指令（算术或逻辑）的向量化方法，是X86-64(AMD64)指令集的拓展。由最早的SSE（SSE有多个改进版本，最后的SSE版本是SSE4.2），到后来的AVX，AVX2，最新的是AVX512，由于AVX512在较新和高端的CPU才支持，相较而言，更多的CPU支持AVX2。

AVX2指令集配合16个256位ymm寄存器（低128位为xmm寄存器），可以存储8个4字节(float, int)或4个8字节(double, long)数据，并在一个指令内同时做8个或4个运算。

Basic Solution

现在来考虑一个任务，将长度均为N = 100000的float数组a, b, c做运算(a + b) * c，存放在数组r中。将这个操作重复NTIME = 50000次，记录消耗的时间。

由于这个任务作为编程问题是如此的简单和直接，我们很快写出c程序如code 0: simd0.c，并觉得毫无优化的余地，因为仅有三行代码没有任何多余的操作，只有17,19,20三行代码是在完成任务。其中的变量sum仅仅是用来防止编译器的自动优化会直接去掉重复50000次这个步骤，虽然这个优化很厉害，但是偏离了我们的原意。

使用命令gcc simd0.c -o simd0.exe编译，运行时间为13s 803ms 846μs，取5次计算的中位数，下同。如果我们满足于这个结果，觉得已经很快了，那生活也就没有了乐趣。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

#include <sys/time.h>
#include <stdio.h>

#define N 100000
#define NTIMES 50000

float a[N], b[N], c[N], r[N];

int main(void)
{
    int i, times, sum;
    struct timeval starttime, endtime;
    long s, ms, us, timeuse;
    gettimeofday(&starttime, 0);
    
    /* workhorse code  */
    for (times = 0; times < NTIMES; times++)
    {
        for (i = 0; i < N; ++i)
            r[i] = (a[i] + b[i]) * c[i];
        sum += sum;
    }

    /* elapsed time  */
    gettimeofday(&endtime, 0);
    timeuse = 1000000 * (endtime.tv_sec - starttime.tv_sec) + 
            endtime.tv_usec - starttime.tv_usec;
    s = timeuse / 1000000;
    ms = timeuse / 1000 % 1000;
    us = timeuse % 1000;
    printf("%ld s %ld ms %ld us\n", s, ms, us);

    return sum;
}

GCC自动优化

我们可以用命令gcc -S -o simd0.s simd0.c生成上述程序simd0.exe的汇编代码，看看哪里可以优化，其实这个汇编的蠢不在于没有用SIMD指令，主要问题是有太多的访问内存指令，代码没有贴出来。

如果我们使用-O1参数进行优化，即gcc -O1 -o simd0-O1.exe simd0.c，消耗的时间下降为3s 841ms 346μs。汇编代码如code 1: simd0-O1.s, 其中略去了计时，打印和某些伪指令，下同。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

main:
  movl $50000, %edx
  jmp .L2
.L6:
  addss %xmm1, %xmm1
  subl $1, %edx
  je .L4
.L2:
  movl $0, %eax
.L3:
  movss a(%rax), %xmm0
  addss b(%rax), %xmm0
  mulss c(%rax), %xmm0
  movss %xmm0, r(%rax)
  addq $4, %rax
  cmpq $400000, %rax
  jne .L3
  jmp .L6
.L4:
  cvttss2si %xmm1, %eax
  ret

simd0-O1.s的代码和不使用优化的代码相比，主要改进是从内存中取一个数比如a[i]，只需要一次访存，因为指标i存在寄存器%rax中，这一个改进大幅降低了访问指标i和数组元素a[i]和其他数组的时间。simd0-O1.s的代码非常符合我们思维，它使用movss从内存中取一个单精度数，使用addss和mulss执行单精度浮点数的加法和乘法。执行一遍计算后再做第二遍。

使用-O2优化和-O1对这个程序没有太大差异，我们略去。

SIMD的第一尝试

如果我们使用-O3级别优化，就会放发生某些不一样的事情，汇编代码见code 2: simd0-O3.s

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

main:
  movl $50000, %edx
.L2:
  xorl %eax, %eax
.L3:
  movaps b(%rax), %xmm0
  addps a(%rax), %xmm0
  addq $16, %rax
  mulps c-16(%rax), %xmm0
  movaps %xmm0, r-16(%rax)
  cmpq $400000, %rax
  jne .L3
  addss %xmm1, %xmm1
  subl $1, %edx
  jne .L2
  cvttss2si %xmm1, %eax
  ret

和code 1: simd0-O2.s相比，code 2: simd0-O3.s代码启用了SSE指令集的SIMD指令，它使用完整的xmm寄存器，而不是低32位，进行浮点数的加法和乘法。其中，movaps从内存中取4个单精度浮点数放在xmm寄存器中，或者相反方向，addps和mulps分别对xmm寄存器中的4个单精度数据做加法和乘法。使用SSE指令，运行时间又降为2s 612ms 38μs。

AVX2上场

虽然我们仍然可以使用参数让编译器自动向量化，比如gcc -O3 -maxv2 -o simd1-avx2-O3.exe simd0.c，但是这样的优化结果是不可预测的，因此在编程的时候使用适当的机制或约定引导编译器进行向量化，是更加稳健的方法。

内存对齐

将变量的存储地址对齐到8/16/32/64等整数的倍数，是某些指令的要求，同时可以使寻址更加有效。在GCC环境下，使用__attribute__((aligned(__BIGGEST_ALIGNMENT__)))修饰符引导编译器进行内存对齐。比如声明长度为N的float数组，对齐到64字节

1

float a[N] __attribute__((aligned(64)));

向量数据类型

在GCC环境下，可以声明基本数据类型char, int, float, double, long的数组作为向量，来引导编译器使用SIMD指令编译，比如下面的代码声明了一种8个float长的向量类型v8f，

1

typedef float v8f __attribute__((vector_size(32)));

向量数据类型可以使用运算符+ - * / & ^ | %等进行逐元素运算，他们会被编译成相应的SIMD指令。

现在我们可以利用上面的技巧，来重写程序，确保我们使用了SIMD指令，见code 3: simd1.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

#include <sys/time.h>
#include <stdio.h>

#define N 100000
#define NTIMES 50000
#define ALIGN 64
typedef float v8f __attribute__((vector_size(32)));
#define VLEN (sizeof(v8f) / sizeof(float))

float a[N] __attribute__((aligned(ALIGN)));
float b[N] __attribute__((aligned(ALIGN)));
float c[N] __attribute__((aligned(ALIGN)));
float r[N] __attribute__((aligned(ALIGN)));

int main(void)
{
    v8f *va, *vb, *vc, *vr;
    int i, times, sum;
    struct timeval starttime, endtime;
    long s, ms, us, timeuse;
    gettimeofday(&starttime, 0);

    /* workhorse code  */
    for (times = 0; times < NTIMES; times++)
    {
        va = (v8f *)a;
        vb = (v8f *)b;
        vc = (v8f *)c;
        vr = (v8f *)r;
        for (i = 0; i < N; i += VLEN){
            *vr = (*va + *vb) * (*vc);
            va++, vb++, vc++, vr++;
        }
        sum += sum;
    }
    
    /* elapsed time  */
    gettimeofday(&endtime, 0);
    timeuse = 1000000 * (endtime.tv_sec - starttime.tv_sec) + 
    endtime.tv_usec - starttime.tv_usec;
    s = timeuse / 1000000;
    ms = timeuse / 1000 % 1000;
    us = timeuse % 1000;
    printf("%ld s %ld ms %ld us\n", s, ms, us);

    return sum;
}

上面的程序使用指向向量数据类型v8f的指针来进行运算，这些操作可以很好的编译成SIMD指令。并且由于向量长度是8的整数倍，所以我们不用额外的操作。使用命令gcc -mavx2 -O3 -S -o simd1-O3.s simd1.c编译，得到和我们预想一样的汇编代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

main:
  xorl %eax, %eax
  movl $50000, %ecx
.L3:
  xorl %edx, %edx
.L2:
  vmovaps a(%rdx), %ymm1
  vaddps b(%rdx), %ymm1, %ymm0
  addq $32, %rdx
  vmulps c-32(%rdx), %ymm0, %ymm0
  vmovaps %ymm0, r-32(%rdx)
  cmpq $400000, %rdx
  jne .L2
  addl %eax, %eax
  subl $1, %ecx
  jne .L3
  vzeroupper
  ret

和code 2: simd0-O3.s不同的是，这里使用指令vmovaps vaddps vmulps对8个（组）float数据进行操作。可惜的是，虽然我们一顿操作猛如虎，但是效率提升并不明显，原因是因为这个例子中，最费时的是内存的读写而不是计算加法和乘法。

如果我们换一个顺序，对于某个i，我们重复NTIMES次r[i] = (a[i] + b[i]) * c[i]操作，然后在进行下一个i，由于编译器优化总是会去掉一层循环，所以我们直接修改汇编代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

main:
.L2:
    movl	$50000, %eax
.L3:
    vmovaps	a(%rdx), %ymm1
    vaddps	b(%rdx), %ymm1, %ymm0
    vmulps	c(%rdx), %ymm0, %ymm0
    .p2align 4,,10
    addl	%ebx, %ebx
    vmovaps	%ymm0, r(%rdx)
    subl	$1, %eax
    jne	.L3
    addq	$32, %rdx
    cmpq	$400000, %rdx
    jne	.L2
    xorl	%edx, %edx
    vzeroupper
    ret

用命令gcc -o simd2-0.exe simd2-0.s将汇编指令链接成二进制程序。虽然只是调换了两个循环的顺序，操作也没有少一个，但是由于L1/2/3 cache的存在，code 5显著降低了访问内存的次数，执行时间仅为0s 634ms 562μs，验证了我们上面关于效率瓶颈的论断。

如果我们显式地先读一次数据到寄存器中，做NTIMES次计算，最后再写到内存中，这样就避免了多次读写内存的操作。汇编代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

main:
.L2:
    movl	$50000, %eax
    vmovaps	a(%rdx), %ymm1
    vmovaps	b(%rdx), %ymm2
    vmovaps	c(%rdx), %ymm3
.L3:
    vaddps	%ymm2, %ymm1, %ymm0
    vmulps	%ymm3, %ymm0, %ymm0
    .p2align 4,,10
    addl	%ebx, %ebx
    vmovaps	%ymm0, r(%rdx)
    subl	$1, %eax
    jne	.L3
    addq	$32, %rdx
    cmpq	$400000, %rdx
    jne	.L2
    xorl	%edx, %edx
    vzeroupper
    ret

最后优化的程序执行时间为0s 324ms 504μs。和最初的13.8s相比，这个速度快了40多倍，虽然主要功劳并不归功于SIMD指令，但是合理的SIMD向量化有助于我们减少内存的访问，从而为程序加速。