20

考虑一种算法来测试在特定次数的尝试后从一组 N 个唯一数字中选择某个数字的概率（例如，在 N=2 的情况下，轮盘赌（没有 0）中需要 X 次尝试的概率是多少黑赢？）。

正确的分布是 pow(1-1/N,X-1)*(1/N)。

但是，当我使用以下代码对此进行测试时，在 X=31 处总是有一个深沟，与 N 无关，也与种子无关。

这是由于使用中的 PRNG 的实现细节而无法避免的内在缺陷，这是一个真正的错误，还是我忽略了一些明显的东西？

// C

#include <sys/times.h>
#include <math.h>
#include <stdio.h>

int array[101];
void main(){

    int nsamples=10000000;
    double breakVal,diffVal;
    int i,cnt;

    // seed, but doesn't change anything
    struct tms time;
    srandom(times(&time));

    // sample
    for(i=0;i<nsamples;i++){
        cnt=1;
        do{
            if((random()%36)==0) // break if 0 is chosen
                break;
            cnt++;
        }while(cnt<100);
        array[cnt]++;
    }

    // show distribution
    for(i=1;i<100;i++){
        breakVal=array[i]/(double)nsamples; // normalize
        diffVal=breakVal-pow(1-1/36.,i-1)*1/36.; // difference to expected value
        printf("%d %.12g %.12g\n",i,breakVal,diffVal);
    }
}

在带有 libc6 包 2.15-0ubuntu20 和 Intel Core i5-2500 SandyBridge 的最新 Xubuntu 12.10 上进行了测试，但几年前我已经在一台较旧的 Ubuntu 机器上发现了这一点。

我还在 Windows 7 上使用 Unity3D/Mono 进行了测试（但不确定是哪个 Mono 版本），这里使用 System.Random 时，沟渠发生在 X=55，而 Unity 的内置 Unity.Random 没有可见沟渠（至少没有对于 X<100)。

分布：在此处输入图像描述

区别：在此处输入图像描述

4

3 回答 3

11

这是因为 glibc 的random()函数不够随机。根据这个页面，对于返回的随机数random()，我们有：

oi = (oi-3 + oi-31) % 2^31

或者：

oi = (oi-3 + oi-31 + 1) % 2^31.

现在取，并假设上面的第一个等式是使用的（每个数字都有 50% 的机会发生这种情况）。现在如果和，那么小于 1/36 的机会。这是因为 50% 的时间会小于 2^31，而当这种情况发生时，xi = oi % 36xi-31=0xi-3!=0xi=0oi-31 + oi-3

xi = oi % 36 = (oi-3 + oi-31) % 36 = oi-3 % 36 = xi-3,

这是非零。这会导致您在 0 个样本之后看到 31 个样本。

于 2013-02-04T01:10:24.457 回答
7

在这个实验中测量的是伯努利实验的成功试验之间的间隔，其中成功被定义random() mod k == 0为一些k（OP 中的 36 个）。random()不幸的是，执行意味着伯努利试验在统计上不独立，这一事实破坏了这一事实。

我们将编写`random()'的输出，我们注意到：rndiith

rndi = rndi-31 + rndi-3    概率为 0.75

rndi = rndi-31 + rndi-3 + 1概率为 0.25

（有关证明大纲，请参见下文。）

让我们假设，我们目前正在研究. 那么它一定是这样的，因为否则我们会把循环算作长度。rndi-31 mod k == 0rndirndi-3 mod k ≠ 0k-3

但是（大多数时候）。(mod k): rndi = rndi-31 + rndi-3 = rndi-3 ≠ 0

因此，当前的试验在统计上并不独立于之前的试验，成功后的第 31 次试验成功的可能性远低于伯努利试验的无偏系列。

使用实际上并不适用于random()算法的线性同余生成器的通常建议是使用高位而不是低位，因为高位“更随机”（即，与连续值的相关性较小）。但这在这种情况下也不起作用，因为上述恒等式同样适用于 functionhigh log k bits和 function mod k == low log k bits。

事实上，我们可能期望线性同余生成器工作得更好，特别是如果我们使用输出的高阶位，因为虽然 LCG 在蒙特卡罗模拟方面不是特别好，但它不会受到线性反馈的影响random().

random算法，对于默认情况：

让state是一个无符号长的向量。使用种子、一些固定值和混合算法进行初始化。为简单起见，我们可以认为状态向量是无限的，尽管只使用了最后 31 个值，因此它实际上是作为环形缓冲区实现的。state0...state30

生成rndi: (Note: ⊕ is addition mod 232.)

statei = statei-31 ⊕ statei-3

rndi = (statei - (statei mod 2)) / 2

Now, note that:

(i + j) mod 2 = i mod 2 + j mod 2    if i mod 2 == 0 or j mod 2 == 0

(i + j) mod 2 = i mod 2 + j mod 2 - 2 if i mod 2 == 1 and j mod 2 == 1

If i and j are uniformly distributed, the first case will occur 75% of the time, and the second case 25%.

So, by substitution in the generation formula:

rndi = (statei-31 ⊕ statei-3 - ((statei-31 + statei-3) mod 2)) / 2

     = ((statei-31 - (statei-31 mod 2)) ⊕ (statei-3 - (statei-3 mod 2))) / 2 or

     = ((statei-31 - (statei-31 mod 2)) ⊕ (statei-3 - (statei-3 mod 2)) + 2) / 2

The two cases can be further reduced to:

rndi = rndi-31 ⊕ rndi-3

rndi = rndi-31 ⊕ rndi-3 + 1

如上所述，第一种情况发生率为 75%，假设 rnd i-31和 rnd i-3是从均匀分布中独立得出的（它们不是，但它是一个合理的第一近似值）。

于 2013-02-04T02:31:03.613 回答
1

正如其他人指出的那样，random()不够随机。

在这种情况下，使用较高位而不是较低位没有帮助。根据手册 ( man 3 rand)，旧的实现rand()在低位有问题。这就是为什么random()建议改为。虽然，当前的实现rand()使用与random().

我尝试了建议的正确使用旧的rand()：

if ((int)(rand()/(RAND_MAX+1.0)*36)==0)

...并在 X=31 处获得相同的深沟

有趣的是，如果我将rand()的数字与另一个序列混合，我会摆脱困境：

unsigned x=0;
//...

        x = (179*x + 79) % 997;
        if(((rand()+x)%36)==0)

我正在使用旧的Linear Congruential Generator。我从素数表中随机选择了 79、179 和 997。这应该生成长度为 997 的重复序列。

也就是说，这个技巧可能引入了一些非随机性，一些足迹......由此产生的混合序列肯定会通过其他统计测试。x在连续迭代中永远不会采用相同的值。实际上，重复每个值需要 997 次迭代。

''[..] 随机数不应使用随机选择的方法生成。应该使用一些理论。”（DEKnuth，“计算机编程的艺术”，第 2 卷）

对于模拟，如果您想确定，请使用Mersenne Twister

于 2013-02-04T11:52:31.407 回答