Máy phát điện số ngẫu nhiên libc thiếu sót?

Question

Xem xét một thuật toán để kiểm tra xác suất mà một số nhất định được chọn từ một tập hợp N số duy nhất sau một số lần thử cụ thể (ví dụ, với N = 2, xác suất trong Roulette (không 0) là bao nhiêu cố gắng cho Black để giành chiến thắng?).

Phân bố chính xác cho điều này là pow (1-1/N, X-1) * (1/N).

Tuy nhiên, khi tôi kiểm tra điều này bằng cách sử dụng mã sau, luôn có một rãnh sâu tại X = 31, độc lập với N và độc lập với hạt giống.

Đây có phải là một lỗ hổng nội tại không thể ngăn chặn do các chi tiết cụ thể của PRNG đang sử dụng, đây có phải là lỗi thực sự hay tôi nhìn thấy điều gì đó hiển nhiên?

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

Thử nghiệm trên một up-to-date Xubuntu 12.10 với libc6 gói 2.15-0ubuntu20 và Intel Core i5-2500 SandyBridge, nhưng tôi phát hiện này đã là một vài năm trước đây trên một máy tính Ubuntu cũ.

Tôi cũng đã thử nghiệm điều này trên Windows 7 bằng Unity3D/Mono (không chắc chắn phiên bản Mono nào), và ở đây mương xảy ra ở X = 55 khi sử dụng System.Random, trong khi Unity's Unity.Random không có mương nhìn thấy được (ít nhất là không cho X < 100).

Sự phân bố:

Sự khác biệt:

Answer 1

Điều này là do chức năng glibc của random() không đủ ngẫu nhiên. Theo this page, cho những con số ngẫu nhiên được trả về bởi random(), ta có:

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

hay:

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

Bây giờ hãy xi = oi % 36 và giả sử phương trình đầu tiên ở trên là phương trình được sử dụng (điều này xảy ra với 50% cơ hội cho mỗi số). Bây giờ nếu xi-31=0 và xi-3!=0, thì cơ hội xi=0 nhỏ hơn 1/36. Điều này là do 50% thời gian oi-31 + oi-3 sẽ ít hơn 2^31, và khi điều đó xảy ra,

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

đó là khác không. Điều này gây ra mương bạn thấy 31 mẫu sau khi lấy mẫu 0.

Answer 2

Điều được đo trong thử nghiệm này là khoảng thời gian giữa các thử nghiệm thành công của thử nghiệm Bernoulli, nơi thành công được xác định là random() mod k == 0 đối với một số k (36 trong OP). Thật không may, nó bị hủy hoại bởi thực tế là việc thực hiện random() có nghĩa là các thử nghiệm Bernoulli không độc lập về mặt thống kê.

Chúng tôi sẽ viết rndi cho ith đầu ra của `ngẫu nhiên()' và chúng tôi lưu ý rằng:

rndi = rndi-31 + rndi-3     với xác suất 0,75

rndi = rndi-31 + rndi-3 + 1 với xác suất 0.25

(Xem dưới đây để biết một phác thảo bằng chứng.)

Giả sử rndi-31 mod k == 0 và hiện tại chúng tôi đang tìm kiếm tại rndi. Sau đó, nó phải là trường hợp rndi-3 mod k ≠ 0, bởi vì nếu không chúng tôi sẽ tính chu kỳ là chiều dài k-3.

Nhưng (hầu hết thời gian) (mod k): rndi = rndi-31 + rndi-3 = rndi-3 ≠ 0.

Vì vậy, bản dùng thử hiện tại không độc lập về mặt thống kê so với các thử nghiệm trước đó, và thử nghiệm 31 ^st sau khi thành công ít có khả năng thành công hơn so với thử nghiệm Bernoulli không thiên vị. Lời khuyên thông thường khi sử dụng các máy tạo đồng đẳng tuyến tính, không thực sự áp dụng cho thuật toán random(), là sử dụng các bit bậc cao thay vì các bit bậc thấp, bởi vì các bit bậc cao là "ngẫu nhiên hơn" "(có nghĩa là, ít tương quan với các giá trị liên tiếp). Nhưng điều đó cũng không có tác dụng trong trường hợp này, bởi vì các thông tin nhận dạng trên giữ tốt như nhau đối với hàm high log k bits đối với hàm mod k == low log k bits. Trên thực tế, chúng ta có thể mong đợi một máy phát điện tuyến tính để làm việc tốt hơn, đặc biệt nếu chúng ta sử dụng các bit thứ tự cao của đầu ra, bởi vì mặc dù LCG không đặc biệt tốt ở mô phỏng Monte Carlo, nó không bị ảnh hưởng phản hồi tuyến tính của random().

---

random thuật toán, đối với trường hợp mặc định:

Hãy state là một vector của chờ đợi unsigned. Khởi tạo state0...state30 sử dụng hạt giống, một số giá trị cố định và thuật toán trộn. Để đơn giản, chúng ta có thể xem xét vector trạng thái là vô hạn, mặc dù chỉ có 31 giá trị cuối cùng được sử dụng để nó thực sự được thực hiện như một bộ đệm vòng.

Để tạo rndi: (Note: ⊕ là bổ sung mod 2)

statei = statei-31 ⊕ statei-3

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

Bây giờ, lưu ý rằng:.

(i + j) mod 2 = i mod 2 + j mod 2    nếu i mod 2 == 0 hoặc j mod 2 == 0

(i + j) mod 2 = i mod 2 + j mod 2 - 2 nếu i mod 2 == 1 và j mod 2 == 1

Nếu i và j được phân bố đều, trường hợp đầu tiên sẽ diễn ra 75% thời gian, và trường hợp thứ hai là 25%.

Vì vậy, bằng cách thay thế trong công thức thế hệ:

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 hoặc

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

Hai trường hợp có thể giảm hơn nữa để:

rndi = rndi-31 ⊕ rndi-3

.210

rnd i = rnd i-31 ⊕ rnd i-3 + 1

Như trên, trường hợp đầu tiên xảy ra 75% thời gian, giả sử rnd rằng _i-31 và rnd _i-3 được độc lập rút ra từ một bản phân phối đồng đều (mà chúng không có, nhưng đó là một xấp xỉ đầu tiên hợp lý).

Answer 3

Như những người khác đã chỉ ra, random() không đủ ngẫu nhiên.

Sử dụng các bit cao hơn thay vì các bit thấp hơn không giúp ích trong trường hợp này. Theo hướng dẫn sử dụng (man 3 rand), cũ việc triển khai rand() gặp sự cố trong các bit thấp hơn. Đó là lý do tại sao random() được đề xuất thay thế. Mặc dù, việc triển khai hiện tại rand() sử dụng cùng một trình phát như random().

tôi đã cố gắng đề nghị đúng sử dụng của cũ rand():

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

... và có những rãnh sâu cùng một lúc X = 31

Interstingly, nếu tôi trộn số rand() 's với một chuỗi khác, tôi thoát khỏi mương:

unsigned x=0; 
//... 

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

Tôi đang sử dụng cũ Linear Congruential Generator. Tôi đã chọn ngẫu nhiên 79, 179 và 997 từ một bảng số nguyên tố. Điều này sẽ tạo ra một chuỗi lặp lại có chiều dài 997.

Điều đó nói rằng, thủ thuật này có thể giới thiệu một số không ngẫu nhiên, một số dấu chân ... x không bao giờ có cùng giá trị trong các lần lặp liên tiếp. Thật vậy, phải mất chính xác 997 lần lặp để lặp lại mọi giá trị.

'' [..] Không được tạo số ngẫu nhiên bằng phương pháp được chọn ngẫu nhiên. Một số lý thuyết nên được sử dụng."(DEKnuth, "The Art of Computer Programming", vol.2)

Đối với mô phỏng, nếu bạn muốn chắc chắn, sử dụng Mersenne Twister