伪随机数与准随机数的比较.pdf

32 计算机与信息技术 软件纵横 伪随机数与准随机数的比较 王水花 张煜东 吴乐南 东南大学 信息科学与工程学院 江苏 南京 210096 摘 要 传统的随机数生成法主要采用逆转法 首先生成 0 1 区间上的均匀分布U 然后令X F 1 U 则 X满足 分布F 然而 这样生成的均匀分布具有明显的差异性 在小样本或高维空间的情况下尤其严重 因此 引进一种新的随机 数生成方法 即准随机数生成器 实验验证了准随机数生成器得到的随机数的差异性优于传统方法 最后 提出一种基于准 随机数生成器的蒙特卡罗积分方法 结果优于传统的蒙特卡罗积分 关键词 伪随机数 准随机数 Kolmogorov Smirnov 假设检验 1 引言 随机数生成算法 1 是一类重要的算法 广泛应用于仿真技 术等场合 然而 目前的伪随机数生成器 Pseudo random number generator PRNG 2 存在一个重要缺陷 即样本分 布与真实分布不一致 这主要发生在以下两种情况 抽样 代价过高 样本数目较少 空间维数较高 3 因此 有必要寻找一类新的随机数发生器 准随机数发 生器 Quasi random number generator QRNG 4 能够生成 稳定 低差异性的 low discrepancy 样本 而与样本数目 或空间维数无关 5 故针对蒙特卡罗积分结果不稳定的情况 提出一种基于 QRNG 的蒙特卡罗积分 发现比传统方法性能 有所提升 2 伪随机数介绍 伪随机数是由确定的算法生成的 其分布函数与相关性 均能通过统计测试 与真实随机数的差别在于 它们是由算 法产生的 而不是一个真实的随机过程 一般地 伪随机数 的生成方法主要有以下 3 种 6 1 直接法 Direct Method 根据分布函数的物理意 义生成 缺点是仅适用于某些具有特殊分布的随机数 如二 项式分布 泊松分布 2 逆转法 Inversion Method 假设U服从 0 1 区间上的均匀分布 令X F 1 U 则X的累计分布函数 CDF 为F 该方法原理简单 编程方便 适用性广 3 接受拒绝法 Acceptance Rejection Method 假设 希望生成的随机数的概率密度函数 PDF 为f 则首先找 到一个 PDF 为g的随机数发生器与常数c 使得f x cg x 然后根据接收拒绝算法求解 由于算法平均运算c次 才能得到一个希望生成的随机数 因此c的取值必须尽可能 小 显然 该算法的缺点是较难确定g与c 因此 伪随机数生成器 PRNG 一般采用逆转法 其基 础是均匀分布 均匀分布 PRNG 的优劣决定了整个随机数体 系的优劣 7 下文研究均匀分布的 PRNG 3 伪随机数生成器的缺点 重复做N 10000 次试验 每次产生S 20 与S 100 个随 机分布的样本 同时采用 Kolmogorov Smirnov 假设检验 hypothesis test 来确定样本是否满足均匀分布 规定 0 假设 null hypothesis 为样本服从均匀分布 1 假设 alternative hypothesis 为样本不服从均匀分布 采用 P 值 0 1 衡量 P 值越趋近于 0 表示越 有理由拒绝 0 假设 即样本不服从均匀分布 P 值越趋近于 1 表示越有理由接受 0 假设 即样本服从均匀分布 如图 1 与图 2 所示 随着 P 值下降 样本也越来越不服 从均匀分布 实践中希望 P 值越大越好 然而统计学的结论 显示 P 值一定服从均匀分布 与N S大小无关 这表明由 于随机性 总会出现某次抽样得到的样本不服从 甚至远离 均匀分布 另外 样本大小的不同 造成检验标准的不同 直观上看 S 100 对应的均匀分布普遍比 S 20 对应的更均匀 因此 小样本情况下均匀分布 PRNG 的差异性尤为严重 软件纵横 计算机与信息技术 33 00 20 40 60 81 0 20 40 60 80 100 120 140 p values Number of Tests 00 20 40 60 81 0 5 10 15 样本值 频数 a KS 假设检验 P 值直方图 b P 0 95 对应结果 00 20 40 60 81 0 5 10 15 样 本 值 频数 00 20 40 60 81 0 5 10 15 20 样本值 频数 c 0 475 P 0 525 对应结果 d P0 95 对应结果 00 20 40 60 81 0 1 2 3 样本值 频数 00 20 40 60 8 0 1 2 3 4 5 样本值 频数 c 0 475 P 0 525 对应结果 d P 0 05 对应结果 图 2 10000 次 S 20 的均匀分布 PRNG 34 计算机与信息技术 软件纵横 1 0 500 51 1 0 5 0 0 5 1 x y a 二维 P 4 b 三维 P 6 图 3 维数灾难示意图 维数灾难 curse of dimensionality 是另一个造成差异性 的潜在原因 通常假设一个半径r维数为d的超球 被一个 边长为 2r维数为d的超立方体所包围 假设超立方体内存在 一个均匀分布的点 则由于超球的体积为 2r d d 2 d d 2 超立方体的体积为 2r d 因此该点同时也落在超球内的概 率P为 2 1 2 2 d d P dd 1 令维数d由 2 逐步增长到 20 则对应的概率P如图 4 所 示 显然 当d 20 时 P仅为 2 46 10 8 因此 若在 2 维 空间中 1 个样本在半径r的意义下能逼近一个正方形 则在 20 维空间内 则需要 1 2 46 10 8 4 06 107个样本才能在半 径r的意义下逼近超立方体 因此 在高维空间中的大样本 甚至不如低维空间中的小样本 2468101214161820 10 6 10 4 10 2 d维数 P概率 图 4 概率 P 与维数 d 的关系 4 准随机数发生器 上节讨论了造成差异性的两个情况 小样本与高维空间 本节讨论如何构建一类新的随机数发生器 使其具有较低的 差异性 PRNG 缺陷的根源在于 随机性 与 均匀性 的矛盾 因此 不要求新的发生器模拟真实的均匀分布 而力求任意 大小的样本 尤其是小样本 都能满足低差异性 换言之 以牺牲随机性为代价 换来均匀性的提高 称其为准随机数 发生器 QRNG 均匀分布 QRNG 的优势在于 其生成的样 本更趋于均匀分布 在其基础上构建的各类分布 包括高斯 分布 的 QRNG 其生成的样本也更趋于服从对应的分布 目前有 3 种准随机序列 Quasi random sequency 可用 来辅助生成均匀分布随机数 分别是 Halton 序列 Sobol 序 列 Latin 超立方体序列 以 Halton 序列为例 8 简单介绍其 原理 假定存在互质基底b1 b2 bs 则整数n可以 展开为 0 i ij i naj n b 2 若记第j个基底的逆函数为 1 j i bij i naj n b 3 则定义序列 12 s nbbb xnnn 4 为 Halton 序列 这样生成的序列有强相关性 为了消除 相关性 一般采用 忽略 skip 跳跃 leap 与 置 乱 scramble 3 个步骤 忽略与跳跃的物理意义见图 5 另 外常采用 RR2 算法对 Halton 序列进行置乱 图 5 忽略与跳跃的物理意义 5 实验 实验在主频 3 2GHz 内存 2GB 的 IBM P4 电脑上运行 实验平台采用 Matlab2009b 并借助 Statistical Toolbox 软件纵横 计算机与信息技术 35 5 1 P 值比较 对QRNG重新进行第0节的试验 S 20 令Skip为1000 Leap 为 100 图 6 显示了基于 QRNG 的 P 值直方图 此时 P 值几乎全部落在点 1 上 放大可见 P 值分布在 0 994 1 区 间内 对比图 6 a 与图 2 a 可见 QRNG 效果优于 PRNG 其差异性更小 任取其中某次抽样结果 绘成直方图 示于图 6 b 可见较图 2 b c d 更为均衡 0 9950 9960 9970 9980 9991 0 200 400 600 800 1000 1200 1400 1600 1800 p values Number of Tests 00 20 40 60 81 0 0 5 1 1 5 2 2 5 3 样 本 值 频数 a KS 假设检验 P 值直方图 b 某次抽样结果直方图 图 6 10000 次 S 20 的均匀分布 QRNG 5 2 抱团现象 图 7 a 给出了另一个直观例子 在二维空间中绘出 500 个随机数的散点图 可见基于 PRNG 的均匀分布有一部分点 互相重叠 反之另一部分区域为空 这种 抱团现象 正是 由于随机性引起 图 7 b 给出了基于 QRNG 的散点图 可 见此时 抱团现象 完全消除 在高维空间中抱团现象能导 致 子超立方体 为空 使得样本分布与均匀分布相差极大 a PRNG b QRNG 图 7 500 个均匀分布的样本点 6 应用 QRNG 可以应用在对随机数的低差异性要求极高的场 合 下面给出一个例子 传统的积分计算通过将自变量划分 为一个个小块 计算对应的函数值 再累加求得 这样对一 维函数或二维函数比较适用 但是对高维函数则无能为力 例 如 对 10 维的函数求积分 假设每维划分 100 个区间 则一 共需要计算大约 100 10 1020 在通常的计算机上根本无法承受 因此 往往采用蒙特卡罗方法 9 求解 与 插针法 类似 以二维为例 考虑平面上一个边长为 1 的正方形及其内部的 一个形状不规则图形 为求出该图形的面积 蒙特卡罗方法 向该正方形随机地投掷N个点 统计落于图形内的M个点 则该图形的面积近似为M N 然而 传统蒙特卡罗方法中的随机数均采用 PRNG 生成 差异性较高 因此提出一种基于 QRNG 的蒙特卡罗方法 即 算法中的随机数采用 QRNG 生成 假设求解 10 维空间中一 个半径为 1 的超球的体积 理论结果是 2 5502 分别运行两 种算法 规定每运行 10 万次记录当前结果 结果如图 8 所示 很明显 传统方法得到的结果偏小 且与理想结果偏差较大 而本文方法的结果更趋于分布在理想结果两侧 且偏差更小 小 明显优于传统方法 0 20 40 60 811 21 41 61 82 x 106 2 48 2 49 2 5 2 51 2 52 2 53 2 54 2 55 2 56 计 算 次 数 积分结果 蒙 特 卡 罗 方 法 理 想 结 果 0 20 40 60 811 21 41 61 82 x 106 2 5 2 51 2 52 2 53 2 54 2 55 2 56 2 57 2 58 2 59 计 算 次 数 积分结果 蒙 特 卡 罗 方 法 理 想 结 果 a 传统方法 b 本文方法 图 8 两种蒙特卡罗方法比较 36 计算机与信息技术 软件纵横 结论 为了弥补 PRNG 的缺陷 介绍了 QRNG 的原理与算法 通过实验验证了 QRNG 优于 PRNG 并提出一种基于 QRNG 的蒙特卡罗积分方法 能够较好地实现高维函数的积分 Maaranen 提出将 QRNG 用于计算遗传算法的初始值 10 发现其性能优于传统的遗传算法 这也指导了我们今后进一 步的研究方向 即将 QRNG 用于其他优化算法 参考文献 1 Dwyer G P Williams K B Portable random number generators J Journal of Economic Dynamics and Control 2003 27 4 645 650 2 Wu P C Huang K C Parallel use of multiplicative congruential random number generators J Computer Physics Communications 2006 175 1 25 29 3 Marchi A Liverani A Giudice A Polynomial pseudo random number generator via cyclic phase J Mathematics and Computers in Simulation 2009 79 11 3328 3338 4 Lemieux C Handbooks in Operations Research and Management Science M 2006 North Holland 351 379 5 Ogawa S Lecot C A quasi random walk method for one dimensional reaction diffusion equations J Mathematics and Computers in Simulation 2003 62 3 487 494 6 Diaz N C Gil A V Vargas M J Assessment of the suitability of different random number generators for Monte Carlo simulations in gamma ray spectrometry J Applied Radiation and Isotopes 2010 68 3 469 473 7 Sun X P Chen Z Z Spherical basi