随机数的随机性(randomness of random number)

云风有两篇关于随机性的文章,我觉得挺有意思的:

coolshell里面提到的 `fisher_yates` 算法,让我想起了 `reservior sampling` 算法,都是一遍下来和交换N-1次,我隐隐地感觉两者之间是有联系的。

先说下第二篇文章里面的一个问题,在N次实验中,连续出现5次正面的概率多大,云风在最后面给了解法:

然后是关于文章里面可视化随机性的程序,可以在wiki上参考 PBM文件格式, 它同时支持ASCII和二进制,还支持位图/灰度/彩色(RGB)三种格式。

python下面有个 fmder/ghalton: Quasi Random Number Generator 准随机数生成器,我对比了一下这个生成器和numpy.random的随机性差异。

这里面有个很小的问题,是如何将(0,1)的浮点数映射成为整数。我没有特别好的办法,采用最简单的概率映射:

这个查找可以做成二分,但是如果范围比较小的话,就直接搜索吧。

def gen_halton_random(low, high, size):
    unit = 1 / (high - low)
    pvts = [i * unit for i in range(1, high - low + 1)]
    pvts[-1] = 1

    def to_int(p):
        for i in range(len(pvts)):
            if p <= pvts[i]:
                return i

    xs = gahlton_gen.get(size)
    ys = [to_int(p[0]) for p in xs]
    return ys


第一个实现是可视化这两个随机数生成器的随机性 代码

第一张图是numpy产生的,第二张图是halton产生的,差异还是蛮明显的

numpy_random.jpg halton_random.jpg

为了更好地产生随机数,ghalton可以指定dimension, 因为我们要产生二维图像所以inDim=2.

ghalton_gen = ghalton.Halton(inDim=2)


def to_int(ps, n):
    unit = 1 / n
    pvts = [i * unit for i in range(1, n)]
    pvts[-1] = 1

    def fx(p):
        for i in range(len(pvts)):
            if p <= pvts[i]:
                return i

    return [fx(p) for p in ps]


def gen_halton_random(n, m, size):
    zs = ghalton_gen.get(size)
    xs = to_int([p[0] for p in zs], n)
    ys = to_int([p[1] for p in zs], m)
    return zip(xs, ys)

第二个实现是对两个随机数分布做卡方检验 代码

下面函数是计算卡方值,deg是自由度

def compute_x2(fn, deg, size):
    n = deg + 1
    xs = fn(0, n, size)
    counter = Counter(xs)
    res = 0
    for k in range(n):
        v = counter[k]
        res += v * v
    return res * n / size - size

两个随机数生成器得到的卡方值差别也蛮大的,这个准随机数生成器效果就是好。

degree = 5, size = 4000
===== numpy random =====
[2.378, 2.891, 3.587, 4.142, 4.97, 6.086, 6.743, 8.852, 11.216, 12.473]
===== halton random =====
[0.002, 0.002, 0.002, 0.002, 0.002, 0.002, 0.002, 0.008, 0.008, 0.008]

卡方分布表可以在wiki里面找到

自由度k \ p-value 0.95 0.90 0.80 0.70 0.50 0.30 0.20 0.10 0.05 0.01 0.001
1 0.004 0.02 0.06 0.15 0.46 1.07 1.64 2.71 3.84 6.64 10.83
2 0.10 0.21 0.45 0.71 1.39 2.41 3.22 4.60 5.99 9.21 13.82
3 0.35 0.58 1.01 1.42 2.37 3.66 4.64 6.25 7.82 11.34 16.27
4 0.71 1.06 1.65 2.20 3.36 4.88 5.99 7.78 9.49 13.28 18.47
5 1.14 1.61 2.34 3.00 4.35 6.06 7.29 9.24 11.07 15.09 20.52
6 1.63 2.20 3.07 3.83 5.35 7.23 8.56 10.64 12.59 16.81 22.46
7 2.17 2.83 3.82 4.67 6.35 8.38 9.80 12.02 14.07 18.48 24.32
8 2.73 3.49 4.59 5.53 7.34 9.52 11.03 13.36 15.51 20.09 26.12
9 3.32 4.17 5.38 6.39 8.34 10.66 12.24 14.68 16.92 21.67 27.88
10 3.94 4.86 6.18 7.27 9.34 11.78 13.44 15.99 18.31 23.21 29.59

对于自由度k=5来说,只有值不大于11.07才能在p<0.05置信区间中(这个置信区间应该是单向的)。也就是说,多次试验当中, 如果卡方值如果出现大于11.07的话,那么说明每个变量并不是相互独立的,反过来说明随机数并不随机。