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$。