(1条消息) 伪蒙特卡洛(Quasi
分享一道由群员“Melbourne”，外号 “Paper Machine”，有数学小王子之称的小伙伴分享的题目！
特别说明：本文非原创，经投稿者同意后发表。
01
PART
算法介绍
期望：在概率论和统计学中，数学期望(mean)（或均值，亦简称期望）是试验中每次可能结果的概率乘以其结果的总和，是最基本的数学特征之一。它反映随机变量平均取值的大小。
题目：在1*1的正方形中随机撒三个点，两两点都可构成长方形的一组对顶点，这样一共有三个长方形，需要求面积第二大的长方形的面积的期望。
算法：每次随机三个点，计算第二大面积，最后统计期望。
02
PART
蒙特卡洛
蒙特卡罗法也称统计模拟法、统计试验法。是把概率现象作为研究对象的数值模拟方法。是按抽样调查法求取统计值来推定未知特性量的计算方法。蒙特卡罗是摩纳哥的著名赌城，该法为表明其随机抽样的本质而命名。故适用于对离散系统进行计算仿真试验。在计算仿真中，通过构造一个和系统性能相近似的概率模型，并在数字计算机上进行随机试验，可以模拟系统的随机特性。
蒙特卡洛方法(Monte Carlo Method) 指的是一类使用随机变量解决概率问题的方法。比较常见的是计算积分、计算概率、计算期望等问题。
常见的蒙特卡洛方法依赖于随机变量的“随机性”，即未发生的事件无法根据已有信息进行预测，比如抛硬币、掷骰子等。在计算机中，常见的随机数是由一系列确定性算法进行生成的，通常称之为伪随机数(pseudo random number)。由于计算精度有限，且这些随机数在统计意义上“不够随机”，会出现可预测的重复序列，这些数在统计意义上收敛精度有限。
与常见的蒙特卡洛方法不同的是，伪蒙特卡洛使用了低差异序列(low discrepancy sequence，常见的有halton序列、sobol序列等)，不使用常见的（伪）随机数，其收敛速率更快（记 N 为样本数量，伪蒙特卡洛收敛速率可达
，而普通蒙特卡洛方法收敛速率仅为
。另一个最重要的性质是伪蒙特卡洛使用的低差异序列是可复现的(replicable)，即不会随环境改变而改变，没有随机种子；而普通蒙特卡洛使用的伪随机数会因随机种子不同而导致结果不同，收敛效果也不尽相同。
03
PART
题目分析
本算法利用伪蒙特卡洛完成。
(CPP代码如下)
1#include <cmath>2#include <cstdio>3#include <vector>4#include <cassert>5#include <omp.h>6const int UP=100;7bool sieve[UP+100];8int primes[UP],top=0;9void init()10{11 for (int i=2;i<=UP;++i)12 if (!sieve[i])13 {14 primes[top++]=i;15 for (int j=i;j<=UP/i;++j)16 sieve[i*j]=true;17 }18}19std::vector<double> halton(long long i,const int &dim)20{21 assert(dim<=top);22 std::vector<double> prime_inv(dim,0),r(dim,0);23 std::vector<long long> t(dim,i);24 for (int j=0;j<dim;++j)25 prime_inv[j]=1.0/primes[j];26 auto f=[](const std::vector<long long> &t)->long long {27 long long ret=0;28 for (const auto &e:t)29 ret+=e;30 return ret;31 };32 for (;f(t)>0;)33 for (int j=0;j<dim;++j)34 {35 long long d=t[j]%primes[j];36 r[j]+=d*prime_inv[j];37 prime_inv[j]/=primes[j];38 t[j]/=primes[j];39 }40 return r;41}42double experiment(long long idx)43{44 std::vector<double> li=halton(idx,6);45 double area1=fabs((li.at(0)-li.at(2))*(li.at(1)-li.at(3)));46 double area2=fabs((li.at(0)-li.at(4))*(li.at(1)-li.at(5)));47 double area3=fabs((li.at(2)-li.at(4))*(li.at(3)-li.at(5)));48 double w=area1+area2+area3-std::max(std::max(area1,area2),area3)-std::min(std::min(area1,area2),area3);49 return w;50 }51const int BATCH=100000;52const int THREADS=40;53int main()54{55 init();56 double total=0;57 for (long long trial=0;;)58 {59 std::vector<double> li(THREADS,0);60 omp_set_dynamic(0);61 omp_set_num_threads(THREADS);62 #pragma omp parallel for63 for (long long thread=0;thread<THREADS;++thread)64 {65 for (long long i=0;i<BATCH;++i)66 li.at(thread)+=experiment(trial+thread*BATCH+i);67 }68 for (const auto &d:li)69 total+=d;70 trial+=THREADS*BATCH;71 printf("%lld: %.10f\n",trial,total/trial),fflush(stdout);72 }73 return 0;74}
分析：使用了并行计算，批量跑随机实验，速度大大提升。其中halton函数会生成halton低差异序列，其值域为[0，1]，参数i表示第i个抽样，dim表示生成数据的维度（本例中每次实验需要6个点，使用6维数据点即可），不同样本之间互不影响，故可使用并行计算提速。
#表示随机试验次数×10^7，Avg表示第二大面积的平均值，Err表示与真实值的绝对误差×10^（-10）。
| # | Avg | Err | # | Avg | Err |
|---|---|---|---|---|---|
| 1 | 0.1017786804 | 55 | 2 | 0.1017786707 | 152 |
| 3 | 0.1017786905 | 46 | 4 | 0.1017786889 | 30 |
| 5 | 0.1017786809 | 50 | 6 | 0.1017786836 | 23 |
| 7 | 0.1017786849 | 10 | 8 | 0.1017786868 | 9 |
| 9 | 0.1017786799 | 60 | 10 | 0.1017786837 | 22 |
| 11 | 0.1017786845 | 14 | 12 | 0.1017786839 | 20 |
| 13 | 0.1017786874 | 15 | 14 | 0.1017786839 | 20 |
| 15 | 0.1017786848 | 11 | 16 | 0.1017786868 | 9 |
| 17 | 0.1017786851 | 8 | 18 | 0.1017786863 | 4 |
| 19 | 0.1017786854 | 5 | 20 | 0.1017786887 | 28 |
| 21 | 0.1017786858 | 1 | 22 | 0.1017786844 | 15 |
| 23 | 0.1017786841 | 18 | 24 | 0.1017786852 | 7 |
| 25 | 0.1017786849 | 10 | 26 | 0.101778684 | 19 |
| 27 | 0.1017786838 | 21 | 28 | 0.1017786852 | 7 |
| 29 | 0.1017786838 | 21 | 30 | 0.1017786846 | 13 |
| 31 | 0.1017786859 | 0 | 32 | 0.1017786862 | 3 |
| 33 | 0.1017786859 | 0 | 34 | 0.1017786853 | 6 |
| 35 | 0.1017786854 | 5 | 36 | 0.1017786859 | 0 |
| 37 | 0.101778685 | 9 | 38 | 0.1017786854 | 5 |
| 39 | 0.1017786853 | 6 | 40 | 0.1017786858 | 1 |
| 41 | 0.1017786848 | 11 | 42 | 0.1017786851 | 8 |
| 43 | 0.1017786847 | 12 | 44 | 0.1017786841 | 18 |
| 45 | 0.101778685 | 9 | 46 | 0.1017786842 | 17 |
| 47 | 0.1017786852 | 7 | 48 | 0.1017786848 | 11 |
| 49 | 0.1017786854 | 5 | 50 | 0.1017786851 | 8 |
| 51 | 0.1017786842 | 17 | 52 | 0.1017786844 | 15 |
可以看到，在实验 次之后，收敛精度可达9位小数，非常精确。由于使用的随机数“不够随机”，普通的蒙特卡洛在同样的实验次数下仅能收敛至五位小数的精度。
上述方法可扩展至其他随机问题中，非常实用且高效，欢迎大家讨论。