作者meto000 (陽貨欲見孔子)
看板Statistics
標題[VB6 ] 亂數真的夠亂嗎?
時間Thu Jun 25 00:08:28 2009
這個問題原先是貼在 VB 版的
不過我想它跟統計應該有點關係
或許這邊的先進會有不同的看法也說不定
所以也貼來這邊請教一下
ps:
我用EXCEL驗證過它的RAND指令
每次計算一行65535個亂數的平均和12倍變方
五次的結果如下
0.499098904 0.999599624
0.499093231 0.997547013
0.501684552 1.003165544
0.496511423 0.995069681
0.500108072 0.998571983
看來EXCEL的亂數好像也只能收斂到兩位小數的樣子
※ [本文轉錄自 Visual_Basic 看板]
作者: meto000 (陽貨欲見孔子) 看板: Visual_Basic
標題: [VB6] 亂數真的夠亂嗎?
時間: Thu Jun 25 00:01:45 2009
大家應該都知道
用 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 來自: 220.140.18.156 (06/25 00:13)
1F：→ yhliu:Var(Xbar) = Var(X)/n = (1/12)/65535 = 0.000001272 06/25 00:37
2F：→ yhliu:sd(Xbar) = 0.00113. 就 "平均數" 結果來看沒有問題. 06/25 00:38
3F：→ yhliu:不過, "是否夠亂" 有許多需要 test 的. 06/25 00:39
4F：→ yhliu:另, 在你的 VB 程式中, 計算的捨入誤差當 n 大時可能重要. 06/25 00:44
5F：→ meto000:感謝回答，我也想過是不是四捨五入的問題，但倍精度有15位 06/25 00:45
6F：→ yhliu:以 X 的累加來說,因小數點對準,後面累加的數字會損失精確度. 06/25 00:46
7F：→ meto000:小數，而平均從80萬之後就已經收斂在0.5001xx了，好怪 06/25 00:46
8F：→ yhliu:而 X^2 的累加精確度損失更嚴重. 06/25 00:48
9F：→ yhliu:VB 結果 mean 幾乎都在 0.5 之上可能與 VB 計算特性有關, 但 06/25 00:50
10F：→ yhliu:因我沒學過 VB, 因此不知其故. 06/25 00:50
11F：→ yhliu:哦...你的 VB 計算是累計的, 難怪 mean 幾乎都在 0.5 之上. 06/25 00:52
12F：→ meto000:感謝,所以我們無法期望由大量取樣而能收斂到小數點很多位 06/25 00:53
13F：→ meto000:由SD(Xbar)來看頂多兩三位而已, 是不是這樣呢? 06/25 00:54
14F：→ yhliu:n=10000000, Var(Xbar)=0.8333E-8, SD(Xbar)=0.0000913 06/25 00:54
15F：→ meto000:所以我取樣1000萬次,就只能有四位的精準度, 是這樣嗎? 06/25 00:55
16F：→ yhliu:Xbar=0.500137 在合理範圍 (0.5(+/-)2*SD) 06/25 00:56
17F：→ yhliu:是的, MonteCarlo 計算並不能取得 "很精確" 的結果, 但可以 06/25 00:57
18F：→ yhliu:"快速" 取得 "尚可" 的結果. 特別適合多變量計算如 R^k 之區 06/25 00:59
19F：→ yhliu:域上的定積分. 06/25 00:59
20F：→ meto000:感謝指教,多變量定積分的模擬我也有作過,只是對精確度有點 06/25 01:03
21F：→ meto000:疑惑，原來不是拿時間跟他拼就能提升多少的，頂多只有 06/25 01:03
22F：→ meto000:mean±2sd範圍內的精確度而已,我還以為VB的亂數歪掉了說 06/25 01:04
23F：→ yhliu:要得精確積分值還是數值積分較實在. 但數值積分當 dimension 06/25 01:06
24F：→ yhliu:提高時所需時間很恐怖! 如一維定積分需分 n 段, k維要分 n^k 06/25 01:08
25F：→ yhliu:個區塊, 而誤差仍不及一維時. 但 MonteCarlo 方法則與維度 06/25 01:09
26F：→ yhliu:相關不大. 換言之,一維積分需取 n 點, k維取個 k*n 可能嫌多 06/25 01:11
27F：→ meto000:了解,montecarlo只是提供難纏問題的另一種思考方向而已 06/25 01:17
28F：→ clickhere:可以好慮Variance reduction的方法可多個2-3位精準度 06/25 09:01
29F：→ clickhere: 考 06/25 09:02
30F：→ meto000:可以麻煩clickhere大大介紹一下作法嗎？ 06/25 23:10
※ meto000:轉錄至看板 Visual_Basic 06/25 23:11
31F：→ clickhere:ask google 06/26 11:05
33F：推 yuekun:不只EXCEL連JAVA內建的亂數都不夠亂吧 網路上有些專門產生 06/28 00:05
34F：→ yuekun:亂數的Library拿來用(例如SSJ) 這個問題討論很久了 06/28 00:05
35F：→ meto000:這裡討論的其實是VB，Excel只是順便試試看而已... 06/28 00:53
36F：→ meto000:關於Variance reduction的部份之前有看過，只是熊熊沒想到 06/28 00:54
37F：→ meto000:就是這個東西... 感謝幾位先進指教 06/28 00:55