李小兵组，Halton序列

                       蒙特卡洛算法中随机点的产生

一  伪随机数介绍

   随机数生成算法是一类重要的算法，广泛应用于仿真技术等场合。

伪随机数是由确定的算法生成的，其分布函数与相关性均能通过统计测试。与真实随机数的差别在于，它们是由算法产生的，而不是一个真实的随机过程。一般地，伪随机数的生成方法很多，有线性同余法，直接法，逆转法等

C/C++中的随机数产生方法：

C/C++语言中伪随机数生成算法实际上是采用了"线性同余法。

具体的计算如下：

 Xi = (Xi-1 * A + C ) mod M 

其中A,C,M都是常数（一般会取质数）。当C=0时，叫做乘同余法。引出一个概念叫seed，它会被作为X0被代入上式中，然后每次调用rand()函数都会用上一次产生的随机值来生成新的随机值。可以看出实际上用rand()函数生成的是一个递推的序列，一切值都来源于最初的 seed。所以当初始的seed取一样的时候，得到的序列都相同。C语言里面有RAND_MAX这样一个宏，定义了rand()所能得到的随机值的范围。在C里可以看到RAND_MAX被定义成0x7fff，也就是32767。rand()函数里递推式中M的值就是32767。线性同余法生成的是伪随机数，粗略符合均匀分布。

函数原型
1、C++标准函数库提供的随机数生成器rand，返回0－RAND_MAX之间均匀分布的伪随机整数。 RAND_MAX值为32767。rand()函数不接受参数，默认以1为种子（即初值）。随机数生成器总是以相同的种子开始，所以形成的伪随机数列也相同，失去了随机意义。（但这样便于程序调试）

2、C++中另一函数srand（），可以指定不同的数（无符号整数变元）为种子。但是如果种子相同，伪随机数列也相同。一个办法是让用户输入种子，但是仍然不理想。

3、 比较理想的是用变化的数，比如时间来作为随机数生成器的种子。time的值每时每刻都不同。所以种子不同，所以，产生的随机数也不同。引入time.h头文件，然后使用srand(time(0))来使用当前时间使随机数发生器随机化，这样就可以保证每两次运行时可以得到不同的随机数序列。

#include

srand(time(0))；

 

二  伪随机算法都存在差异性，不均匀性。因此，不要求新的发生器模拟真实的均匀分布，而力求任意大小的样本（尤其是小样本）都能满足低差异性。换言之，以牺牲随机性为代价，换来均匀性的提高，称其为准随机数发生器。

目前有3种准随机序列可用来辅助生成均匀分布随机数，分别是Halton序列、Sobol序列、Latin超立方体序列。以Halton序列为例^[8]，简单介绍其原理。

假设是以质数2为基底的Halton序列，范围在0～1之间。则从产生的序列为

1/2  1/4  3/4  1/8  5/8  3/8  7/8  1/16  9/16……..

怎么产生的呢？

3=1+2

5=1+4

7=3+4

9=1+8

………..

假设以质数3为基底，Halton序列如下

1/2  2/3  1/9  4/9  7/9  2/9  5/9  8/9  1/27  10/27  19/27……..

 

在二维中，产生的点的序列就是（1/2，1/3）（1/4，2/3）（3/4，1/9）（1/8，4/9）（5/8，7/9）（3/8，2/9）（7/8，5/9）（1/16，8/9）（9/16，1/27）（，）（，）….

三  简单比较伪随机算法和准随机产生

 

  伪随机序列产生的点                        准随机序列产生的点

在二维空间中绘出个随机数的散点图。可见基于伪随机算法的均匀分布有一部分点互相重叠，反之另一部分区域为空，这种“抱团现象”正是由于随机性引起。而基于准随机算法的散点图，可见此时“抱团现象”完全消除。在高维空间中抱团现象能导致“子超立方体”为空，使得样本分布与均匀分布相差极大。

四  应用举例

以二维为例，考虑平面上一个边长为1的正方形及其内部的一个形状不规则图形，为求出该图形的面积，蒙特卡罗方法向该正方形随机地投掷N个点，统计落于图形内的M个点，则该图形的面积近似为M/N。

    然而，传统蒙特卡罗方法中的随机数均采用PRNG生成，差异性较高。因此提出一种基于准随机的蒙特卡罗方法，即算法中的随机数采用准随机生成。假设求解10维空间中一个半径为1的超球的体积，理论结果是2.5502。分别运行两种算法，规定每运行10万次记录当前结果，结果如图8所示。很明显，传统方法得到的结果偏小，且与理想结果偏差较大；而本文方法的结果更趋于分布在理想结果两侧，且偏差更小小，明显优于传统方法。

             

（a）传统方法

（b）本文方法