史上最全采样方法详细解读与代码实现
史上最全采样方法详细解读与代码实现
1.什么是采样
在信号系统、数字信号处理中，采样是每隔一定的时间测量一次声音信号的幅值，把时间连续的，模拟信号转换成时间离散、幅值连续的采样信号。如果采样的时间间隔相等，这种采样称为均匀采样。
在计算机系统中，有一个重要的问题就是给定一个概率分布p(x) ， 我们如何在计算机中生成它的样本。平时我们接触比较多的场景是，给定一堆样本数据，求出这堆样本的概率分布p(x)。而采样刚好是个逆命题：给定一个概率分布p(x)，如何生成满足条件的样本？
2.均匀分布采样
均匀分布式是一种最简单的分布。在计算机中生成[0,1]之间的伪随机数序列，就可以看做是一种均匀分布。随机数生成的方法有很多，比较简单的一种方式比如：
xn+1=(axn+c) mod mx_{n+1} = (ax_n + c) \ mod \ mxn+1=(axn+c) mod m
当然计算机中产生的随机数一般都是伪随机数，不过在绝大部分场景下也够用了。
3.离散分布采样
离散分布是比均匀分布稍微复杂一点的情况。例如令p(x)=[0.1,0.2,0.3,0.4]p(x) = [0.1, 0.2, 0.3, 0.4]p(x)=[0.1,0.2,0.3,0.4]，那么我们可以把概率分布的向量看做一个区间段，然后从x∼U(0,1)x \sim U_{(0,1)}x∼U(0,1)分布中采样，判断x落在哪个区间内。区间长度与概率成正比，这样当采样的次数越多，采样就越符合原来的分布。
举个例子，假设p(x)=[0.1,0.2,0.3,0.4]p(x) = [0.1, 0.2, 0.3, 0.4]p(x)=[0.1,0.2,0.3,0.4]，分别为"hello", “java”, “python”, "scala"四个词出现的概率。我们按照上面的算法实现看看最终的结果
import numpy as np
from collections import defaultdict
dic = defaultdict(int)
def sample():
u = np.random.rand()
if u <= 0.1:
dic["hello"] += 1
elif u <= 0.3:
dic["java"] += 1
elif u <= 0.6:
dic["python"] += 1
else:
dic["scala"] += 1
def sampleNtimes():
for i in range(10000):
sample()
for k,v in dic.items():
print k,v
sampleNtimes()
- 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
代码输出结果：
python 3028
java 2006
hello 971
scala 3995
- 1
- 2
- 3
- 4
上面的采样方法，基本满足我们的要求。
4. Box-Muller算法
如果概率密度分布函数p(x)p(x)p(x)是连续分布，如果这个分布可以计算累积分布函数(Cumulative Distribution Function, CDF)，可以通过计算CDF的反函数，获得采样。
如果随机变量U1U_1U1,U2U_2U2是IID独立同分布，且U1,U2∼Uniform[0,1]U_1, U_2 \sim Uniform[0, 1]U1,U2∼Uniform[0,1]，那么有:
Z0=−2lnU1−−−−−−√cos(2πU2)Z_0 = \sqrt{-2ln U_1} cos(2 \pi U_2)Z0=−2lnU1