转载—--Alias Method离散分布随机取样

本文转载自天空的城

在图的随机游走中,有一块需要随机取样, 比如当前到达v节点,那么下一次随机会到达哪个节点。这种问题其实就是离散分布的随机变量的取样问题。 查了一些资料, 发现Alias Method是一种很高效的方式。在初始好之后,每次取样的复杂度为O(1)。 下面介绍下这种方法的使用,具体原理暂时没有深究,有兴趣转Darts,Dice, and coids

Alias-Method 原理

大部分资料都会以这个例子讲解:游戏中经常遇到按一定的掉落概率来随机掉落指定物品的情况,例如掉落 银币25%,金币20%, 钻石10%,装备5%,饰品40%。现在要求下一个掉落的物品类型,或者说一个掉落的随机序列,要求符合上述概率。

其实方法很简单,最直接的方法就是直接取样,先构建一个累加概率分布列表:[0.25,0.45,0.55, 0.60, 1.0],之后产生一个随机数(0-1),假设为0.7,那么在列表中找到第一个大于0.7的数为1.0,对应的类别是第5类。 这种方法很简单直接,但是每次取样复杂度为$O(K)$,使用二分搜索之后,降为$O(\log K)$。

论文中经常使用的是另一种很巧妙的方法: Alias Method,它在初始化之后每次随机取样的复杂度为O(1)。 下面以Darts,Dice, and coids文中例子用图示说明整个步骤,原理太繁杂,不做介绍,参考博客即可:

假设概率分布为: $\frac{1}{2},\frac{1}{3},\frac{1}{12}, \frac{1}{12}$。

  • 初始概率分布: 类别数目K=4K=4,以颜色表示不同的类别

  • 每个类别概率乘以K=4,使得总和为4. 这样分为两类,大于1:第一列与第二列; 小于1: 第三列与第四列。

  • 下面通过拼凑,使得每一列的和都为1,但是每一列中,最多只能是两种类型的拼凑,就是每一列最多两种颜色存在

    • 将第一列拿出$\frac{2}{3}$给最后一列,使其变为1,如下:(棕色表示空缺)

    • 将第一列拿出$\frac{2}{3}$给第三列,使之变为1,如下:

    • 最后一次,第二列给第一列$\frac{1}{3}$, 最后每一列都是1,且每一列最多两种类型,其中下面一层表示原类的概率,上面层表示另外一种类型的概率,如只有一种比如第二列,那么第二层就是NULL:

    用图表示如下:

​ 到此为止,得到Prob, Alias表示初始化完成。

  • 采样过程: 随机取某一列k(即[1,4]的随机整数),再随机产生一个[0-1]的小数c,如果Prob[k] 大于 c,那么采样结果就是k,反之则为Alias[k]。

需要说明的是,该过程一定可以结束(原文有证明)。此外在初始化完成之后每次采样的复杂度为O(1),因此应用很广。

举个例子说明采样过程,比如随机取得第1列, 随机产生的小数为$0.5<\farc{2}{3}$,那么采样的结果就是第一类。 如果随机产生的小数为$0.8>\frac{2}{3}$,那么采样结果就是第一列的第二层的类别,也就是Alias[1]=2(紫色对应的类别: 第二列)。

再简单验证下,看看该方法的采样是不是满足原始的概率分布$\frac{1}{2},\frac{1}{3},\frac{1}{12.}, \frac{1}{12}$。

对于采到第一种的概率,采到第一种有三部分组成,第一列,第三,四列分别求概率求和: $\frac{2}{3}\frac{1}{4}+(1-\frac{1}{3})\frac{1}{4}+(1-\frac{2}{3})*\frac{1}{4}=\frac{1}{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
import numpy as np
import numpy.random as npr

def alias_setup(probs):
'''
probs: 某个概率分布
返回: Alias数组与Prob数组
'''
K = len(probs)
q = np.zeros(K) # 对应Prob数组
J = np.zeros(K, dtype=np.int) # 对应Alias数组
# Sort the data into the outcomes with probabilities
# that are larger and smaller than 1/K.
smaller = [] # 存储比1小的列
larger = [] # 存储比1大的列
for kk, prob in enumerate(probs):
q[kk] = K*prob # 概率
if q[kk] < 1.0:
smaller.append(kk)
else:
larger.append(kk)

# Loop though and create little binary mixtures that
# appropriately allocate the larger outcomes over the
# overall uniform mixture.

# 通过拼凑,将各个类别都凑为1
while len(smaller) > 0 and len(larger) > 0:
small = smaller.pop()
large = larger.pop()

J[small] = large # 填充Alias数组
q[large] = q[large] - (1.0 - q[small]) # 将大的分到小的上

if q[large] < 1.0:
smaller.append(large)
else:
larger.append(large)

return J, q

def alias_draw(J, q):
'''
输入: Prob数组和Alias数组
输出: 一次采样结果
'''
K = len(J)
# Draw from the overall uniform mixture.
kk = int(np.floor(npr.rand()*K)) # 随机取一列

# Draw from the binary mixture, either keeping the
# small one, or choosing the associated larger one.
if npr.rand() < q[kk]: # 比较
return kk
else:
return J[kk]

K = 5
N = 100

# Get a random probability vector.
probs = npr.dirichlet(np.ones(K), 1).ravel()

# Construct the table.
J, q = alias_setup(probs)

# Generate variates.
X = np.zeros(N)
for nn in xrange(N):
X[nn] = alias_draw(J, q)
Reference

http://shomy.top/2017/05/09/alias-method-sampling/

写的还不错?那就来个红包吧!
0%