从伪随机数的产生到高大上的蒙特卡洛算法（C语言实现）

一 准备
1 生成任意区间任意大小的伪随机数

2 什么是蒙特卡洛算法

二蒙特卡洛算法的实现
1 pi的蒙特卡洛计算方式

2 特殊图形的蒙特卡洛计算方式

通过这篇短文想说明两个道理：

看似高大上、神秘兮兮的算法，都是paper tiger；

计算机的计算方式（动辄几Ghz的主频）简直就是为蒙特卡洛度身定做；

一、 准备

1.1 生成任意区间任意大小的伪随机数

C语言中的

rand()

更深远的意义在于其对应于数学（概率论）中的均匀分布（uniformed distributed）。

C语言生成伪随机数的函数：

int rand(void);

该函数随机生成

0~RAND_MAX

之间内的整数：

#define RAND_MAX 0x7fff  // 0x7fff == 32767

产生随机数需要设置种子：

void srand(unsigned int _Seed);

这两个函数所在的头文件是

stdlib.h

或者

cstdlib

，后者又被包含在

iostream

头文件中。

有了

rand()

这个可以生成

0-RAND_MAX

随机数（整数）的函数，经过一定的四则运算和取模运算，便可很容易地得到任意区间的随机数。

以生成

(2, 5)

之间的随机数（可整可小）为例：

double x = 3*(double(rand())/RAND_MAX)+2;

先通过

double(rand())/RAND_MAX

使随机数区间转换为

(0, 1)

，再通过一定的伸缩平移实现对任意区间的仿真，这里的

double

类型转换不可省略，否则整数之间的除法运算得到的结果仍是整数。

double

vs

float

两者的区别在于对浮点数表达的精度不同，

double

是双精度，

float

为单精度。

sizeof(double) == 8;
sizeof(float) == 4;

c语言中的浮点数默认是

double

类型的，除非显式声明为

l

（或者

L

）

float x = 1l;

long

vs

int

16位系统：long是4字节，int是2字节

32位系统：long是4字节，int是4字节

64位系统：long是8字节，int是4字节

更详细的讨论见

long

vs

int

。

更详细的内容参见之前的一篇博客C++伪（pseudo）随机数生成及简单应用。

1.2 什么是蒙特卡洛算法？

这部分内容会比较枯燥，如果读不下去，可先看后边的实验，再读这部分内容会很容易理解和接受。

蒙特卡罗方法（Monte Carlo method），也称统计模拟方法，是一种数值计算方法。是二十世纪四十年代中期由于科学技术的发展和电子

计算机

的发明，而被提出的一种以

概率统计理论

为指导的一类非常重要的

数值计算

方法。是指使用

随机数

（或更常见的

伪随机数

）来解决很多计算问题的方法。

通常蒙特卡罗方法：因所求解的问题本身具有内在的随机性，借助计算机的运算能力可以直接模拟这种随机的过程。

更详尽的解释，参见Monte Carlo method, 蒙特卡洛算法在机器学习中的应用可参见增强学习-蒙特卡洛方法。

二、蒙特卡洛算法的实现

2.1 π的蒙特卡洛计算方式

Fig. 1. 求π的近似值

如图所示的正方形其面积A=1；还有14的圆，其面积是π4。这里如何利用几何图形的概率特性，即蒙特卡洛算法，来近似计算圆周率π的值呢。

想象这是一张纸，其中的圆弧线，将纸划分为两部分，在下雨时将这张纸放置室外，经过一段时间，雨点落在14圆的个数为C，落在整张纸上的雨点个数为D。则有

CD=AB=π41=π4

可得π=4CD

可通过对大量重复随机实验来仿真或者近似计算CD的真实值。让计算机产生随机数x,y,x≤1,y≤1, 模拟雨点的分布情况。这里的关键问题是如何表示或者判断雨点落在扇形区域，即：

x2+y2−−−−−−√≤1

易写出如下的程序：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
int main(int, char**)
{
long c = 0, d = 0, N = 10000;
double x = 0, y = 0, pi = 0;
srand(unsigned(time(NULL)));
for (long i = 0; i < N; ++i)
{
d += 1;
x = double(rand()) / RAND_MAX;
y = double(rand()) / RAND_MAX;
if (sqrt(x*x + y*y) <= 1)    // x^2表示异或；
c += 1;
}
printf("π= %f\n", 4.*c / d);
return 0;
}

这里有一份迭代出来的近似值：

iteration	π
100	2.9600
1000	3.116364
10000	3.150270
100000	3.138326
1000000	3.139696
10000000	3.141699
100000000	3.141511
1000000000	3.141521

可见随着计算迭代次数的增加，估算的精度越来越高。

2.2 特殊图形的蒙特卡洛计算方式

Fig. 2. 计算区域B的面积

继续沿用计算π的思路，模拟雨点落在阴影区域B的概率。此时的阴影区域应满足：

对圆心在(0,0)的扇形而言，

x2+y2−−−−−−√>1

对圆心在(1,0)的扇形而言，

(x−1)2+y2−−−−−−−−−−−√>1

转化为程序语言即是：

if (sqrt(x*x+y*y)>1 && sqrt((x-1)*(x-1)+y*y)>1)
c+=1;

真实的区域面积应当等于：

1−2∗π12−3√4≈0.0434

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
int main(int, char**)
{
long c = 0, d = 0, N = 100000;
double x = 0, y = 0, pi = 0;
srand(unsigned(time(NULL)));

for (unsigned i = 0; i < N; ++i)
{
d += 1;
x = double(rand()) / RAND_MAX;
y = double(rand()) / RAND_MAX;
if (sqrt(x*x + y*y) > 1 && sqrt((x - 1)*(x - 1) + y*y) > 1)
c += 1;
}
printf("s = %f\n", double(c) / d);
return 0;
}