libc隨機數發生器有瑕疵？

Question

考慮一種算法來測試在特定次數的嘗試後從一組N個唯一數字中挑選出某個數字的概率（例如，在N = 2的情況下，輪盤賭中的概率是多少（無0）試圖讓黑方獲勝？）。

正確的分佈是pow（1-1/N，X-1）*（1/N）。

但是，當我使用下面的代碼測試它時，總是在X = 31處有一個深溝，與N無關，並且獨立於種子。

這是一個內在的缺陷，由於使用PRNG的實現細節無法防止，這是一個真正的bug，還是我忽略了一些明顯的東西？

// 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); 
    } 
} 

測試上了最新的Xubuntu 12.10與libc6的軟件包2.15-0ubuntu20和Intel Core i5-2500 SandyBridge的，但我在幾年前就已經發現了這個舊的Ubuntu的機器上。

我也在Windows 7上使用Unity3D/Mono測試了這個（不知道哪個單聲道版本），這裏溝渠在使用System.Random時發生在X = 55，而Unity的內置Unity.Random沒有可見的溝渠（至少不適用於X < 100）。

分佈：

的差異：

Answer 1

這是由於不是glibc的的random()功能足夠隨機的。據this page，由random()返回的隨機數，我們有：

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

或：

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

現在取xi = oi % 36，並假設上面的第一個公式是使用的（這發生在每個數字有50％的機會）。現在如果xi-31=0和xi-3!=0，那麼xi=0小於1/36的機會。這是因爲時間oi-31 + oi-3 50％將小於2^31，並且當發生這種情況，

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

這是非零。這會導致你在0個樣本後看到31個樣本。

Answer 2

在這個實驗中測量的是伯努利實驗的成功試驗之間的間隔，其中一些k（OP中的36）的成功定義爲random() mod k == 0。不幸的是，random()的實施意味着伯努利試驗不是統計獨立性的事實，這一事實令人沮喪。

我們將編寫rndi爲`隨機（）」的ith輸出，我們注意到：

rndi = rndi-31 + rndi-3    的概率爲0.75

rndi = rndi-31 + rndi-3 + 1的概率爲0。25

（請參閱下面的證明輪廓。）

讓我們假設rndi-31 mod k == 0，我們目前正在觀察rndi。那麼一定是rndi-3 mod k ≠ 0的情況，因爲否則的話，我們會把這個週期算作長度k-3。 （大部分時間）(mod k): rndi = rndi-31 + rndi-3 = rndi-3 ≠ 0。因此，目前的試驗在統計學上並不依賴於以前的試驗，成功後的試驗成功的可能性要比無偏差的伯努利試驗系列試驗成功的可能性小得多。

使用線性同餘發生器的常用建議，其實際上並不適用於random()算法，因爲高位是「更隨機的」，所以使用高位而不是低位。 「（即與連續值較少相關）。但是在這種情況下，這也不起作用，因爲對於功能，功能mod k == low log k bits的上述身份同樣適用。實際上，我們可能會希望線性同餘發生器能夠更好地工作，特別是如果我們使用輸出的高階位，因爲雖然LCG在蒙特卡洛模擬中不是特別好，但它不會受到線性反饋random()。

---

random算法，默認情況下：

讓state是無符號多頭的載體。使用種子初始化state0...state30，一些固定值和混合算法。爲了簡單起見，我們可以認爲狀態向量是無限的，儘管只有最後31個值被使用，所以它實際上被實現爲環形緩衝區。

爲了產生rndi: (Note: ⊕是加法模2 32 ）

statei = statei-31 ⊕ statei-3

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

現在，請注意：

(i + j) mod 2 = i mod 2 + j mod 2   如果i mod 2 == 0或j mod 2 == 0

(i + j) mod 2 = i mod 2 + j mod 2 - 2如果i mod 2 == 1和j mod 2 == 1

i如果和j均勻地分佈，第一殼體將發生的75％的時間，和第二殼體25％。

所以，由生成公式中替換：

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或

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

這兩種情況可以進一步簡化爲：

rndi = rndi-31 ⊕ rndi-3

RND 我 = RND I-31 ⊕ RND I-3 + 1

如上所述，第一種情況發生的時間的75％，假定RND _I-31和rnd _i-3獨立地從均勻分佈（它們不是，但它是合理的第一近似）獨立地繪製。

Answer 3

正如其他人指出的，random()不是隨機的。

在這種情況下，使用較高位而不是較低位不起作用。根據手冊（man 3 rand），舊的rand()的實現在較低的位中有問題。這就是推薦使用random()的原因。儘管rand()的當前實現使用與random()相同的生成器。

我試過老rand()推薦正確使用：

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

...並得到了相同的深溝裏，在X = 31

Interstingly，如果我混rand()的數字與另一個序列，我擺脫了溝：

unsigned x=0; 
//... 

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

我使用舊的Linear Congruential Generator。我從素數表中隨機選擇了79,179和997。這應該會生成一個長度爲997的重複序列。

這就是說，這個技巧可能會引入一些非隨機性，一些足跡......由此產生的混合序列肯定會失敗其他統計測試。 x在連續迭代中從不採用相同的值。事實上，重複每個值需要997次迭代。

''[..]隨機數字不應隨隨機選擇的方法生成。有些理論應該被使用。」（DEKnuth，‘計算機程序設計藝術’，下冊）

爲了仿真，如果你想可以肯定，使用Mersenne Twister