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

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

1. 看似高大上、神秘兮兮的算法，都是paper tiger；
2. 计算机的计算方式（动辄几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;}