[VB6] 亂數真的夠亂嗎?

[meto000 (陽貨欲見孔子)]

大家應該都知道
用 X=Rnd 指令就可以得到一個介於0~1之間的均勻分布亂數
也就是 X~Uniform(0,1)
然而 VB 只能產生均勻分布的亂數
因此小弟最近寫了個各種分布的亂數產生程式來用
但發現生出來的東西似乎不是很精確
因此進一步對 Rnd 這個指令加以檢驗
才發現大事不妙

理論上
由於X~Uniform(0,1)
故E(X)=1/2, Var(X)=1/12
因此
當我們反覆取n個亂數後
再計算其平均和變方
隨著 n->∞
所得的平均和變方應該會趨近1/2和1/12
這也是所謂蒙地卡羅法(Monte Carlo)和拔靴法(Bootstrap)等的理論基礎
在各種科學領域的應用多不勝數
最簡單的蒙地卡羅法應用
大概就屬在單位平面中逢機取點來估算圓周率了
這個我在高中時就玩過了
應該不需要在這邊班門弄斧才對

離題了
話說現在要驗證Rnd是否真的符合均勻分布
於是我寫了這樣的小程式加以驗證：

Open "rnd.txt" For Output As #1
  X1 = 0: X2 = 0: Randomize
  For I = 1 To 10000000
      X = Rnd
      X1 = X1 + X
      X2 = X2 + X ^ 2
      If I Mod 100000 = 0 Then
            EX = X1 / I
            VarX = (X2 - EX ^ 2 * I) / I
            'Print #1, I, X
            Print #1, I, EX, VarX * 12
     End If
  Next I
Close #1

全部的變數都宣告為倍精度
這個程式基本上就是反覆取1000萬個亂數
並同時累積計算 韈X 和 韈X^2
然後每隔10萬個就輸出一次平均EX和變方VarX*12
結果如下：
 100000        0.500341235609055           0.998641767848771
 200000        0.499809045362473           1.00078372215905
 300000        0.500107454074224           0.998437901086435
 400000        0.500276461744309           0.998385742258827
 500000        0.500290068372726           0.999079838731831
 600000        0.500419607292811           1.00021925961703
 700000        0.50026993564742            1.00019207389547
 800000        0.50015473200798            1.000139160073
 900000        0.500144928914176           1.00052572698579
 1000000       0.500105085889816           1.00056944920825
  略
 9000000       0.500131491163254           0.999874573474332
 9100000       0.500137886596617           0.999861888329618
 9200000       0.500145418414655           0.999857654740129
 9300000       0.500141381108992           0.999838666710201
 9400000       0.50014871621416            0.999868309864149
 9500000       0.50014508351587            0.99988266186096
 9600000       0.500142511142095           0.999940878307397
 9700000       0.50014840750419            0.999973340103397
 9800000       0.500148654139577           0.999944556581512
 9900000       0.500138692927601           0.999910212634631
 10000000      0.500136948072433           0.999931865426466

按理
所得的平均和12倍變方應該會趨近於0.5和1
然而
一開始做10萬次模擬發現沒有收斂(在0.49-0.51間亂跳)
做到1000萬次也只收斂到小數第四位
而且是有偏歪的 0.50013-0.50014
也難怪轉換出來的各種分布亂數是有問題的了
這麼一來
用VB來做模擬研究
不就有問題了嗎?

那是 Rnd 指令本身有問題
還是我的程式哪邊搞錯了呢?
可以請大家幫忙看看嗎?

請問各位先進
有人知道 VB 的 Rnd 亂數產生指令是怎麼運作的嗎?
它是經由怎樣的公式運算產生的呢？
不知哪邊可以查得到?
懇請賜教，感激不盡

--

迷時師渡  悟時自渡

                      ～ 六祖惠能


--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 220.140.18.156
※ meto000:轉錄至看板 Statistics                                   06/25 00:08
→ MOONRAKER:簡單虛擬亂數產生器大多是 XOR + Shift
→ MOONRAKER:或者最多是 R(n) = ( R(n-1) * K ) mod M 這種形式
→ MOONRAKER:你得出rnd()的bias可能是因為浮點進位的誤差所致
→ MOONRAKER:但我認為你可以自己寫一個輸出整數的亂數產生器
→ MOONRAKER:簡單的可以找超有名的numerical recipes系列來看