伪蒙特卡洛(Quasi-Monte Carlo, QMC)随机
分享一道由群员“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 for
63 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位小数,非常精确。由于使用的随机数“不够随机”,普通的蒙特卡洛在同样的实验次数下仅能收敛至五位小数的精度。
上述方法可扩展至其他随机问题中,非常实用且高效,欢迎大家讨论。
推荐一些其他的算法文章:
漫画:位运算系列篇(只出现一次的数字 - 进阶版)
漫画:位运算系列篇(只出现一次的数字)、
漫画:百度从Google学来的面试题,想进大厂必备!
漫画:美团面试题(整数拆分)
所以,今天的问题你学会了吗?评论区留下你的想法!
小浩算法,每日一题
关注领取《图解算法》高清版
进群的小伙伴请加右侧私人微信(备注:进群)
伪蒙特卡洛(Quasi-Monte Carlo, QMC)随机相关推荐
- Variance Reduction Methods: a Quick Introduction to Quasi Monte Carlo——完结
https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods- ...
- 蒙特卡洛积分(Monte Carlo Integration)应用:利用蒙特卡洛积分生成 McBeth表
蒙特卡洛积分(Monte Carlo Integration)应用 蒙特卡洛积分 通常函数f(x)的积分: 可以解释为计算函数曲线下方的面积: 而我们的蒙特卡洛积分则是通过近似的方式来获取一个函数的积 ...
- 蒙特卡洛(Monte Carlo)法求定积分
蒙特卡洛(Monte Carlo)法是一类随机算法的统称.随着二十世纪电子计算机的出现,蒙特卡洛法已经在诸多领域展现出了超强的能力.在机器学习和自然语言处理技术中,常常被用到的MCMC也是由此发展而来 ...
- Stata: 蒙特卡洛模拟(Monte Carlo Simulation)没那么神秘
作者:侯新烁 湘潭大学 [编译] (知乎 | 简书 | 码云) Stata连享会 精彩推文1 || 精彩推文2 资料参考来源: The Stata Blog » Monte Carlo simulat ...
- MATLAB蒙特卡洛(Monte Carlo)方法求椭圆面积
MATLAB蒙特卡洛方法求椭圆面积 代码 代码 在某个规定的范围内随机打点,找到满足条件的点,并数一下这些点的数量与总的随机点数量的比,就OK了.关键是设置条件. 代码 clear;clc; n=10 ...
- 【路径追踪】数学工具--蒙特卡洛方法(Monte Carlo)
Intro 蒙特卡洛方法是一类通过随机采样来求解问题的算法, 要求解的问题是某随机事件的概率或某随机变量的期望. 现在认为最早记载的一个蒙特卡洛计算示例是由蒲丰在 1777 年完成的投针试验. 在实验 ...
- java计算椭圆的面积_java算法3_蒙特卡洛方法(Monte Carlo method)求PI和椭圆面积
蒙特卡洛方法,是一种以概率统计理论为指导的一类非常重要的数值计算方法.是指使用随机数来解决很多计算问题的方法.蒙特卡洛方法的名字来源于摩纳哥的一个城市蒙特卡洛,该城市以×××业闻名,而蒙特卡洛方法正是 ...
- 心得复述知识体系:《强化学习》中的蒙特卡洛方法 Monte Carlo Methods in Reinforcement Learning
前言: 刚刚读完 Sutton 的<强化学习(第二版)>第5章:蒙特卡洛方法.为了巩固本章收获,笔者将在本文中用尽量简单直白的语言复述本章的思想,各个知识点之间的关系.同时,这方便笔者日后 ...
- 强化学习(四) - 蒙特卡洛方法(Monte Carlo Methods)及实例
强化学习(四) - 蒙特卡洛方法(Monte Carlo Methods)及实例 4. 蒙特卡洛方法 4.1 蒙特卡洛预测 例4.1:Blackjack(21点) 4.2 动作价值的蒙特卡洛估计 4. ...
最新文章
- K8s Ingress Provider 为什么选择 MSE 云原生网关?
- Oracle入门(十)之概要文件
- Mvc NuGet 数据迁移
- 计量经济学计算机输出结果,计量经济学作业答案A..doc
- 提交注册信息到数据库中
- JAVA连接数据库 遍历集合数组!!!
- 常见排序算法的时间复杂度汇总
- fluidsim元件库下载_FluidSIM
- DC-DC电压基准芯片和REF芯片
- 怎么关闭惠普暗影精灵OMEN 8的主机灯
- sonarqube如何导入规则_sonar如何添加自定义JAVA规则
- win10插入耳机没声音,如何设置声音
- C# 制作贪吃蛇小游戏,最简单的实现
- 中路由怎么配置_【国土语录 第64期】茶队点金大鱼后发生了什么;否定Fy关于茶队中路看法...
- Spark教程——(10)Spark SQL读取Phoenix数据本地执行计算
- 用 Elasticsearch 统计做了几次核酸检测?怎么破?
- 易语言编写“文本文档”
- 便利蜂智能制作策略平台的探索与实践
- 动态规划—1.3 九宫格最短路径
- 怎么解决Tekla节点内部材质改变却外部无法改变