作者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