Alias Sampling

1 简介

alias sampling 是一种离散概率分布采样算法。所谓采样就是通过数学方法按照指定概率分布得到一组或多组样本。离散概率分布采样就是采样得到的样本服从离散概率分布。

一般的离散概率分布采样方法是基于某个一维随机变量服从均匀分布的思想(后面会讲到),然后通过顺序查找或者二分查找的方式确定样本值。

均匀分布是其他基本所有分布的基础。如果一个随机变量服从均匀分布,那么就可以通过它生成与这个变量相关的随机变量,使其服从其他分布。

alias sampling 算法将随机变量扩展到二维,且巧妙地使用“割补”的方式在保持概率不变的情况下将每个维度的大小保持一致,从而降低时间复杂度。

2 基石

2.1 如何生成随机数

目前基本所有的随机数生成算法都是基于一个种子进行的,比如线性同余法(最简单且最快)、梅森旋转算法(号称能够生成最好的伪随机数)。如果每次生成随机数序列所使用的种子相同,那么得到的序列完全相同。也因此,以这种方式生成的随机数被称为“伪随机数”。

理论上,只要保证每次产生随机数前的种子不一样,就可以生成完全不同的随机数序列。C/C++ 中一些做法是通过 time 函数来生成种子,该函数会返回当前距离 UTC 时间的秒数。但如果两次产生随机数的操作时在 time 函数返回的时间内完成的,也即两次产生随机数时用的种子一样,就会产生相同的随机数。

C++11 开始,提供了 random_device(非确定性随机数生成设备),在类 Unix 下会读取 /dev/urandom 设备获取种子。类 Unix 下会搜集环境信息,比如时钟中断、WiFi 信号强弱、鼠标位置变化、键盘敲击速度等,并维护在一个 熵池中。详见:知乎 - 电脑取随机数是什么原理,是真正的随机数吗?

/dev/urandom 产生的随机数可通过如下命令查看: cat /dev/urandom | od -X | head -n 20

因此,在某种程度上来说,可以通过这种方式产生 真随机数

2.2 根据随机生成数生成器构造另一个生成器

给定一个生成 [a, b] 范围内随机数(服从均匀分布)的函数 f,如何生成 [c, d] 范围内随机数(服从均匀分布)的函数 g(其中,a <= c < d <= b)?

方法是:通过 f 生成一个数,然后判定其是否在 [c, d] 之间。如果是,则返回这个数;否则,继续用 f 生成一个数。直到生成的数在 [c, d] 之间为止。伪代码如下:

1
2
3
4
5
6
7
int g() {
int x;
do {
x = f();
} while (x <= d && x >= c);
return x;
}

那么这样得到的数是否在 [c, d] 间等概率呢?证明如下:

假设,a = c = 1,b = 7,d = 5。对于 [1, 5] 间的任意一个数 m,判定生成它的概率。

  • 如果是第一次调用 f 生成了这个数,则概率为 1/7;
  • 如果是第二次调用 f 生成了这个数,那说明第一次调用 f 生成的数没在 [1, 5] 之间(否则第一次就返回了),则概率为 2/7 * 1/7;
  • 如果是第三次调用 f 生成了这个数,那说明第二次调用 f 生成的数没在 [1, 5] 之间(否则第二次就返回了),则概率为 (2/7)^2 * 1/7
  • ......

那么将所有情况下的概率求和,得到:

\[ \begin{align*} p(X = m) & = \frac{1}{7} + \frac{2}{7} \times \frac{1}{7} + (\frac{2}{7})^{2} \times \frac{1}{7} + ... \\ & = \frac{1}{7} \times \sum_{n = 0}^{\infty} (\frac{2}{7})^{n} \\ & = \frac{1}{7} \times \lim_{n \to \infty}{\frac{1 - (\frac{2}{7})^{n}}{1 - \frac{2}{7}}} \\ & = \frac{1}{5} \end{align*} \]

因此,可以知道通过这种方式生成的所有 [1, 5] 范围的数的概率均为 1/5。

详见:Hawstein's Blog - CTCI 解析

3 一般的离散概率分布采样算法

3.1 问题转换

假设现在要按照 {1/2, 1/2} 的概率分布来进行采样。

最为通常的,可以通过公平地掷硬币的方式来得到样本,硬币正面或者背面分别代表两种样本值。

更为普遍的,可以通过数学的方式进行。假设有一个 (0, 1) 的开区间,可通过计算机随机生成这个范围内的数(服从均匀分布)。因此,数落在 (0, 1/2) 这个范围的概率为 1/2,落在 [1/2, 1) 这个范围的概率为 1/2。这两个范围分别代表两种样本值。

这种方法还可以扩展到随机变量取值有 2 种但不相同,或者两种以上的情况。

假设一个概率分布为 {1/2, 1/3, 1/12, 1/12},那么可以在 (0, 1) 上构造多个区间,时间复杂度为 \(O(n)\),如下所示:

Screen Shot 2018-07-25 at 8.10.29 PM.png-17.4kB

只要能够随机地生成 [0, 1) 的数(服从均匀分布),那么就一定能够保证:

  • 数落在 [0, 1/2)(对应一种样本值)的概率为 1/2
  • 数落在 [1/2, 5/6)(对应一种样本值)的概率为 1/3
  • 数落在 [5/6, 11/12)(对应一种样本值)的概率为 1/12
  • 数落在 [11/12, 1)(对应一种样本值)的概率为 1/12

(1/2, 5/6, 11/12) 这个序列是有序的,因此不仅可以顺序查找,还可以采用二分查找。假设采用数组的方式进行存储,则空间复杂度为 \(O(n)\),顺序查找的时间复杂度为 \(O(n)\),二分查找的时间复杂度为 \(O(logn)\)

3.2 代码

二分查找代码如下:

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
//
// Created by Jason on 2018/10/5.
//

#include <iostream>
#include <vector>
#include <random>

using namespace std;

/**
* 普通离散概率分布采样算法 (二分查找)
*/

/**
* 二分查找, 找到 seq 中从左至右第一个比 target 小的数(比 target 大最小数)的下标 - 1 (即为样本值).
* @param seq 构造的有序序列
* @param target 随机生成的 [0, 1) 之间的实数
* @return 样本值 (原概率分布中的数值的下标)
*/
int binary_search(const vector<double> &seq, double target);

/**
* 根据 probability_vec 构造 seq
* @param probability_vec 给出的离散分布
* @param seq 累积值. 0, p1, p1 + p2, ..., 1
*/
void init(const vector<double> &probability_vec, vector<double> &seq);


int binary_search(const vector<double> &seq, double target) {

if (seq.empty() || target > seq.back()) {
return -1;
}

int begin = 0, end = static_cast<int>(seq.size() - 1);
while (begin <= end) {
int mid = begin + (end - begin) / 2;
if (seq[mid] > target) {
end = mid - 1;
} else {
begin = mid + 1;
}
}

return begin - 1;
}

void init(const vector<double> &probability_vec, vector<double> &seq) {
// 构造有序序列
// 0, p1, p1 + p2, ..., 1
for (const auto &num : probability_vec) {
seq.push_back(seq.back() + num);
}
}


int main() {

vector<double> probability_vec = {1.0 / 2, 1.0 / 3, 1.0 / 12, 1.0 / 12};
vector<double> seq = {0};
init(probability_vec, seq);

random_device rd;
default_random_engine random(rd());
uniform_real_distribution<double> distribution(0.0, 1.0);
for (int i = 0; i < 1000000; i++) {
// 随机生成一个在 [0, 1) 间服从均匀分布的的实数
double rand_num = distribution(random);
// 二分查找. 找到 seq 中第一个比 target 大的数的下标 - 1 (即为样本值).
cout << binary_search(seq, rand_num) << endl;
}

return 0;
}

4 Alias Sampling 算法

4.1 问题转换

一般的离散采样算法是将分段放在一维空间上进行处理和计算的,时间复杂度依赖于离散概率分布中点的个数。而 alias sampling 算法将其扩展到二维空间,且时间复杂度为 \(O(1)\)

下面仍然以概率分布 {1/2, 1/3, 1/12, 1/12} 为例。将其扩展到二维空间,如下图所示:

Screen Shot 2018-10-05 at 2.16.23 PM.png-15.2kB

每个小矩形的宽度相同,高度为其对应概率值。根据均匀分布概率计算规则,只要保证 (X, Y) 在所有的小矩形构造的空间 A 中服从均匀分布,然后以此概率分布来生成随机数,那么就能够推导出样本值落在各个小矩形框的概率与它们的面积成正比。

由于 X 和 Y 相互独立,因此可以各自单独生成随机数。但是各个矩形的高度不一,因此这样很难方便地生成在 A 范围内的数。一个想法就是只要能够构造一个包围它们的大矩形 B,且某个二维联合随机变量在 B 范围内服从均匀分布,然后根据 2.2 中的方法,就一定能够得到 (X, Y) 在 A 服从均匀分布。最简单的构造如下:

Screen Shot 2018-10-05 at 2.33.02 PM.png-16.2kB

但是这样仍然没达到目的。每次随机生成一个 B 内的数,然后判定是否在各个小矩形范围内。而且这个过程还是循环进行的。所以看上去似乎比一般的离散采样算法更复杂了。如何解决这两个问题呢?

令 B 中除了 A 的部分为 C。由于随机得到的数落在各个小矩形范围内的概率只与其面积成正比,而与位置无关。因此,一个自然的想法就是对 A 中的部分进行拆分,然后填充到 C 中。这样不能直接办到,因为 B 就是由 A 和 C 构造得到的。因此需要减小 B 的范围,然后将 A 在 B 的部分填充到 B 内,使得 A 和 B 重合,A 的面积正好等于 B 的面积。这样,就不需要循环过程了。

为了更为直观,一般会先将每个矩形的高度扩展为原来的 k 倍(k 为矩形个数),使得 A 的面积为原来的 k 倍。虽然各个小矩形的面积放大了,但由于总面积也跟着放大了,因此生成的数落在每个小矩形的概率仍然是不变的。然后再构造出同样面积的 B。这样构造完成后,B 的高度正好是 1

但是这种构造方式有很多,选择哪种方式最好呢?一个原则是,构造后的每列中最多只有两个小矩形构成,它们的高度之和为 1(这样在数值上看上去就更为直观了)。这样就只需要根据这两个小矩形在高度上的分界点,一次性地判定生成的点是落在哪个小矩形内了。

对于本例,k = 4,则构造后的 B 如下所示:

ret.jpg-29.7kB

解决了这两个问题,那么就可以在 \(O(1)\) 时间复杂度内知道生成的点是属于哪个颜色的小矩形内的,进而可以采样出对应颜色的值。

这样,就可以将原问题转换为如何构造 B 了。其实也就是如何构造两个数组:

  • 一个数组设为 prob,存储每列中原来的小矩形被选择的概率,也即第二层的值。比如这里就是 [1, 2/3, 1/3, 1/3];
  • 另一个数组设为 alias,如果某列需要由其他列补充, 则存储补充过来的列的下标。比如这里就是 [-1, 0, 0, 1],其中 -1 表示不需要补充

可以验证下这种构造是否和原问题等价。可以计算出随机生成的数落在蓝色区域的概率: \[\frac{1}{4} \times 1 + \frac{1}{4} \times \frac{1}{3} + \frac{1}{4} \times \frac{2}{3} = \frac{1}{2}\]

其他类似。

4.2 判定过程

假设现在已经构造出了这两个数组,具体的 判定过程 如下:

  • 先随机生一个 [0, k - 1](k 为数组长度,即概率分布中值的个数)范围内的数用于选出第 i 列(横坐标),概率为 1/k,此时也就可以得到两个数组中对应列的值 prob[i] 和 alias[i];
  • 然后再随机生成一个指定高度范围内的数
    • 如果生成的数小于 prob[i],则选择第 i 列对应的颜色;
    • 否则,选择 alias[i] 对应的颜色。

时间复杂度为 \(O(1)\),空间复杂度为 \(O(1)\)

4.3 如何构造两个数组

伪代码如下:

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
// 每个元素扩大 p.length 倍。同时构造两个栈,一个存储 p 中元素小于 1 的下标,另一个存储 p 中元素大于等于 1 的下标
stack<int> small, large;
for (int i = 0; i < p.length; i++) {
p[i] *= p.length;
if (p[i] < 1) {
small.push(i);
} else {
large.push(i);
}
}

// 构造两个长为 p.length 的数组。prob 存储每列中原来的小矩形被选择的概率,alias 存储每个矩形中另外拆分过来的小矩形的颜色
double[] prob = new double[p.length];
double[] alias = new double[p.length];

while (small 和 large 均不为空) {
分别取出 small 和 large 的栈顶元素 l 和 g
更新 prob[l] = p[l];
更新 alias[l] = g;
p[g] = p[g] − 1 + p[l];
if (p[g] < 1) {
small.push(g);
} else {
large.push(g);
}
}

while (large 非空) {
取出 large 的栈顶元素 g;
更新 prob[g] = 1;
// 表示不需要其他颜色填充
更新 alias[g] = -1;
}

while (small 非空) {
取出 small 的栈顶元素 l;
更新 prob[l] = 1;
// 表示不需要其他颜色填充
更新 alias[l] = -1;
}

时间复杂度为 \(O(n)\),空间复杂度为 \(O(n)\)

以概率分布 {1/2, 1/3, 1/12, 1/12} 为例:

初始化:

init.jpg-53.5kB

从 small 和 large 中分别取出栈顶元素,对应 p 中的下标为 3 和 1,p[3] 为 1/3,p[1] 为 4/3。更新 prob[3] = 1/3,alias[3] = 1。更新 p[1] = p[1] - 1 + p[3] = 2/3。由于 2/3 < 1,因此将 1(粉色)插入到 small 栈顶。处理完成后如下图所示:

1.jpg-57.1kB

从 small 和 large 中分别取出栈顶元素,对应 p 中的下标为 1 和 0,p[1] 为 2/3,p[0] 为 2。更新 prob[1] = 2/3,alias[1] = 0。更新 p[0] = p[0] - 1 + p[1] = 5/3。由于 5/3 >= 1,因此将 0(蓝色)插入到 large 栈顶。处理完成后如下图所示:

2.jpg-60.5kB

从 small 和 large 中分别取出栈顶元素,对应 p 中的下标为 2 和 0,p[2] 为 1/3,p[0] 为 5/3。更新 prob[2] = 1/3,alias[2] = 0。更新 p[0] = p[0] - 1 + p[2] = 1。由于 1 >= 1,因此将 0(蓝色)插入到 large 栈顶。处理完成后如下图所示:

3.jpg-58kB

由于 small 已经空了,因此只处理 large。取出 large 栈顶元素,对应 p 中的下标为 0,p[0] 为 1。更新 prob[0] = 1,alias[] = -1(表示不需要用其他颜色补充)。处理完成后如下图所示:

4.jpg-41.4kB

由于 large 也空了,因此整个过程处理完成。

注:构造出的结果并不只一种。算法中只要初始化时的遍历顺序稍微改下,结果就不一样。但无论用哪种方法,都不影响采样结果的正确性。

4.4 代码

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
//
// Created by Jason on 2018/10/6.
//

#include <iostream>
#include <vector>
#include <stack>
#include <random>

using namespace std;

/**
* alias sampling
*/


/**
* 初始化, 扩大 p 中每个元素为原来的 p.size() 倍, 同时构造 small, large, prob, alias
* @param p 概率分布数值
* @param small p 中小于 1 的元素的下标
* @param large p 中大于等于 1 的元素的下标
* @param prob 存储每列中原来的小矩形被选择的概率
* @param alias 如果某列需要由其他列补充, 则存储补充过来的列的下标
*/
void init(vector<double> &p, stack<int> &small, stack<int> &large, vector<double> &prob, vector<int> &alias);

/**
* 根据构造出的 prob 和 alias, 以及随机得到列标和实数, 来抽样得到样本值 (原概率分布中的数值的下标)
* @param col 用于选择列
* @param p [0, 1) 之间的实数
* @param prob 存储每列中原来的小矩形被选择的概率
* @param alias 如果某列需要由其他列补充, 则存储补充过来的列的下标
* @return 样本值 (原概率分布中的数值的下标)
*/
int search(int col, double p, const vector<double> &prob, const vector<int> &alias);


void init(vector<double> &p, stack<int> &small, stack<int> &large, vector<double> &prob, vector<int> &alias) {

for (int i = 0; i < p.size(); i++) {
p[i] *= static_cast<int>(p.size());
if (p[i] < 1) {
small.push(i);
} else {
large.push(i);
}
}

while (!small.empty() && !large.empty()) {
int small_idx = small.top();
small.pop();
int large_idx = large.top();
large.pop();
prob[small_idx] = p[small_idx];
alias[small_idx] = large_idx;

p[large_idx] -= (1 - p[small_idx]);
if (p[large_idx] < 1) {
small.push(large_idx);
} else {
large.push(large_idx);
}
}
while (!small.empty()) {
int idx = small.top();
small.pop();
prob[idx] = 1;
alias[idx] = -1;
}
while (!large.empty()) {
int idx = large.top();
large.pop();
prob[idx] = 1;
alias[idx] = -1;
}
}

int search(int col, double p, const vector<double> &prob, const vector<int> &alias) {

if (col >= prob.size() || col < 0) {
return -1;
}

return p < prob[col] ? col : alias[col];
}


int main() {

vector<double> p = {1.0 / 2, 1.0 / 3, 1.0 / 12, 1.0 / 12};
stack<int> small;
stack<int> large;
vector<double> prob(p.size());
vector<int> alias(p.size());
init(p, small, large, prob, alias);

random_device rd;
default_random_engine col_random(rd());
default_random_engine random(rd());
uniform_int_distribution<> col_distribution(0, static_cast<int>(p.size() - 1));
uniform_real_distribution<double> distribution(0.0, 1.0);

for (int i = 0; i < 1000000; i++) {
// 随机生成一个在 [0, p.size() - 1] 间服从均匀分布的整数
int col = col_distribution(col_random);
// 随机生成一个在 [0, 1) 间服从均匀分布的的实数
double p_temp = distribution(random);
// 采样
cout << search(col, p_temp, prob, alias) << endl;
}
}

Reference