蒙特卡罗方法完整教程(WORD文档内附有源码)
蒙特卡罗方法 (Monte Carlo simulation)教材
2019/4/15
Monte Carlo模拟
9
2.Monte Carlo方法简史
Stanislaw Ulam (1909-1984)
S. Ulam is credited as the inventor of Monte Carlo method in 1940s, which solves mathematical problems using statistical sampling.
1
Monte Carlo模拟
第一章 引言 （Introduction)
1. 2. 3. 4. Monte Monte Monte Monte Carlo方法 Carlo方法简史 Carlo模拟的应用 Carlo算法的主要组成部分
2019/4/15
Monte Carlo模拟
2
1.Monte Carlo方法
2019/4/15
Monte Carlo模拟
5
2.Monte Carlo方法简史 Buffon投针实验
1768年，法国数学家Comte de Buffon利用投针实验估 计的值
L
d
p
2019/4/15
d
2L
Monte Carlo模拟
6
2.Monte Carlo方法简史
Problem of Buffon’s needle: If a needle of length l is dropped at random on the middle of a horizontal surface ruled with parallel lines a distance d > l apart, what is the probability that the needle will cross one of the lines?
蒙特拉罗方法
用该方法计算π的基本思路是： 1 根据圆面积的公式： s=πR^2 ,当R=1时， S=π。 由于圆的方程是：x^2+y^2=1(x^2为x的平方 的意思),因此1/4圆面积为x轴、y轴和上述方 程所包围的部分。 如果在1*1的正方形中均匀地落入随机点，则 落入1/4圆中的点的概率就是1/4圆的面积。其 4倍，就是圆面积。 由于半径为1，该面积的值为π的值。
蒙特卡罗(Monte Carlo)方法 蒙特卡罗 方法
蒙特卡洛法是什么？ 蒙特卡洛法是什么？
• 蒙特卡洛 蒙特卡洛(Monte Carlo)方法，或称计算机随机 方法， 方法 模拟方法，是一种基于“随机数”的计算方法。 模拟方法，是一种基于“随机数”的计算方法。 这一方法源于美国在第二次世界大战中研制原 子弹的“曼哈顿计划” 该计划的主持人之一、 子弹的“曼哈顿计划”。该计划的主持人之一、 数学家冯·诺伊曼用驰名世界的赌城 诺伊曼用驰名世界的赌城—摩纳哥的 数学家冯 诺伊曼用驰名世界的赌城 摩纳哥的 Monte Carlo—来命名这种方法，为它蒙上了 来命名这种方法， 来命名这种方法 一层神秘色彩。 一层神秘色彩。
m G内，则随机点落入G内的概率 I ≈ n
一道积分题
• 一道证明积分不等式的题：
π
1 1 −(x2 +y2 ) 23 (1− ) < ∫ e dxdy < 4 e 0 30
• 中间的积分值可以用蒙特卡洛法求得 • 因为它是一个二重积分，其几何直观为一 个立体的体积，很巧的是它可以完全包含 于一个棱长为1的正方体中，我们在其中产 生随机点，其中落于所求体积的点数与正 方体中产生的点数之比即为所求的积分值。
1
需要计算的积分为I = ∫ f ( x)dxΒιβλιοθήκη ,积分I等于图中的面积G。0
蒙特卡罗方法
)
其中 c0 2.515517, c1 0.802853,c2 0.010328;
d1 1.432788; d2 0.189269; d3 0.001308
§2 随机数的产生和随机变量的抽样
随机变量的抽样
连续型随机变量的抽样：
3.正态随机变量的抽样方法
(2)基于中心极限定理的方法
N i 1
Xi )2
§1 蒙特卡洛方法该概述---减小误差
减小方差的各种技巧:
显然当给定置信度α（λα）后，误差ε由σ和N决定。要减小ε：
（1）增大试验次数N。在σ固定的情况下，要把精 度提高一个数量级，试验次数N需增加两个数量级。 因此，单纯增大N不是一个有效的办法。
N
（2）减小估计的均方差σ，比如降低一半，那误差就减小一半，这 相当于N增大四倍的效果。
理论依据： 大数定理：均匀分布的算术平均收敛于真值 中心极限定理：置信水平下的统计误差
两个例子： Buffen投针实验求π 射击问题（打靶游戏）
§1 蒙特卡罗方法概述---基本思想
Buffon 投针问题:平面上画很多平行线，间距为 a．向此平面投掷长为l (l < a) 的
针,此针与任一平行线相交的概率 p。
数r,计算满足条件：p(i-1) r p(i)的i值，所
对应的i即为离散型随机变量的一个样本值。
§2 随机数的产生和随机变量的抽样
随机变量的抽样
离散型随机变量的抽样：
2.泊松分布的抽样方法
若N是服从泊松分布的离散型随机变量，其取值为n的概
率为
P(N n) e n (n 0,1, 2,...)
3.计算机方法
根据数论方法，通过数学递推公式运算来实现。 这种方法得到的随机数其实是一种伪随机数。
（完整word版）蒙特卡洛方法及其在风险评估中的应用
（完整word版）蒙特卡洛⽅法及其在风险评估中的应⽤蒙特卡洛⽅法及其应⽤1风险评估及蒙特卡洛⽅法概述１.１蒙特卡洛⽅法。
蒙特卡洛⽅法，⼜称随机模拟⽅法或统计模拟⽅法，是在20世纪40年代随着电⼦计算机的发明⽽提出的。
它是以统计抽样理论为基础，利⽤随机数，经过对随机变量已有数据的统计进⾏抽样实验或随机模拟，以求得统计量的某个数字特征并将其作为待解决问题的数值解。
蒙特卡洛模拟⽅法的基本原理是：假定随机变量X1、X2、X3……X n、Y，其中X1、X2、X3……X n 的概率分布已知，且X1、X2、X3……X n、Y有函数关系：Y=F（X1、X2、X3……X n），希望求得随机变量Y的近似分布情况及数字特征。
通过抽取符合其概率分布的随机数列X1、X2、X3……X n带⼊其函数关系式计算获得Y的值。
当模拟的次数⾜够多的时候，我们就可以得到与实际情况相近的函数Y的概率分布和数字特征。
蒙特卡洛法的特点是预测结果给出了预测值的最⼤值,最⼩值和最可能值，给出了预测值的区间范围及分布规律。
１.２风险评估概述。
风险表现为损损益的不确定性，说明风险产⽣的结果可能带来损失、获利或是⽆损失也⽆获利，属于⼴义风险。
正是因为未来的不确定性使得每⼀个项⽬都存在风险。
对于⼀个公司⽽⾔，各种投资项⽬通常会具有不同程度的风险，这些风险对于⼀个公司的影响不可⼩视，⼩到⼀个项⽬投资资本的按时回收，⼤到公司的总风险、公司正常运营。
因此，对于风险的测量以及控制是⾮常重要的⼀个环节。
风险评估就是量化测评某⼀事件或事物带来的影响的可能程度。
根据“经济⼈”假设，收益最⼤化是投资者的主要追求⽬标，⾯对不可避免的风险时，降低风险，防⽌或减少损失，以实现预期最佳是投资的⽬标。
当评价风险⼤⼩时，常有两种评价⽅式：定性分析与定量分析法。
定性分析⼀般是根据风险度或风险⼤⼩等指标对风险因素进⾏优先级排序，为进⼀步分析或处理风险提供参考。
这种⽅法适⽤于对⽐不同项⽬的风险程度，但这种⽅法最⼤的缺陷是在于，在多个项⽬中风险最⼩者也有可能亏损。
Monte Carlo(蒙特卡洛方法)
P(n 1) U P(n)
则令 X取值
xn.
例1:
离散型随机变量X有如下分布律: X 0 1 2 P(x) 0.3 0.3 0.4 设 U1 ,U 2 ,,U 是 (0,1)上均匀分布的随机数,令 N
0, 0 U i 0.3 xi 1, 0.3 U i 0.6 2, 0.6 U i
ˆ f n ( A) 。 在 n 次中出现的频率。假如我们取 fn ( A) 作为 p P( A) 的估计，即 p
ˆ 然后取 2l a.s. ˆ fn ( A) 作为 的估计。根据大数定律，当 n 时， p p. af n ( A) 2l P 。这样可以用随机试验的方法求得 的估计。历史上 af成器的周期 长度是 10，而后两个生成器的周期长度只有 它的一半。我们自然希望生成器的周期越长 越好，这样我们得到的分布就更接近于真实 的均匀分布。
在给定 m 的情况下，生成器的周期与 a 和 初值 x0 (种子)选择有关。
线性同余生成器(混合同余法) (Linear Congruential Generator )
证明： 由 F 1 (U ) 的定义和均匀分布的分布函数可得： P ( X x) P ( F 1 (U ) x) P (U F ( x )) F ( x )
由定理 1 ，要产生来自 F ( x) 的随机数，只要先 产生来自U (0,1) 随机数 u ，然后计算 F 1 (u ) 即 可。具体步骤如下：
一般形式： xi 1 (axi c) mod m ui 1 xi 1 / m
1. c是非负整数.通过适当选取参数c可以改善 随机数的统计性质(独立性,均匀性).
2. 线性同余器可以达到的最长周期为 m 1 ，我们 可以通过适当的选择 m 和 a ，使无论选取怎样的 初值 x0 都可以达到最大周期(一般选取 m 为质数)
蒙特卡罗法
统计学模拟法之一
01 概念
03 优缺点 05 应用举例
目录
02 基本思路 04 步骤
蒙特卡罗法也称统计模拟法、统计试验法。是把概率现象作为研究对象的数值模拟方法。是按抽样调查法求 取统计值来推定未知特性量的计算方法。蒙特卡罗是摩纳哥的著名赌城，该法为表明其随机抽样的本质而命名。 故适用于对离散系统进行计算仿真试验。在计算仿真中，通过构造一个和系统性能相近似的概率模型，并在数字 计算机上进行随机试验，可以模拟系统的随机特性。
解：希望能用某种方法把我方将要对敌人实施的20次打击结果显示出来，确定有效射击的比率及毁伤敌方火 炮的平均值。这是一个概率问题，可以通过理论计算得到相应的概率和期望值。但这样只能给出作战行动的最终 静态结果，而显示不出作战行动的动态过程。
为了显示我方20次射击的过程，必须用某种方式模拟出以下两件事：一是观察所对目标的指示正确或不正确； 二是当指示正确时，我方火力单位的射击结果。对第一件事进行模拟试验时有两种结果，每一种结果出现的概率 都是1/2。因此，可用投掷1枚硬币的方式予以确定。当硬币出现正面时为指示正确，反之为不正确。对第二件事 进行模拟试验时有3种结果，毁伤1门火炮的可能为1/3，毁伤2门火炮的可能为1/6，没能毁伤敌火炮的可能为1/2。 这时，可用投掷骰子的办法来确定，如果出现的是1、2、3三个点则认为没能击中敌人，如果出现的是4、5点则 认为毁伤敌1门火炮，如果出现6点则认为毁伤敌2门火炮。
应用举例
在我方某前沿防守地域，敌人以1个炮兵排(含两门火炮)为单位对我方进行干扰和破坏。为躲避我方打击， 敌方对其指挥所进行了伪装并经常变换射击地点。经过长期观察发现，我方指挥所对敌方目标的指示有50%是准 确的，而我方火力单位在指示正确时，有1/3的射击效果能毁伤敌人1门火炮，有1/6的射击效果能全击的过程动态地显现出来。
Monte-Carlo模拟教程
举例
例1 在我方某前沿防守地域，敌人以一个炮排（含两门火炮） 为单位对我方进行干扰和破坏．为躲避我方打击，敌方对其阵地 进行了伪装并经常变换射击地点．
经过长期观察发现，我方指挥所对敌方目标的指示有50％是准 确的，而我方火力单位，在指示正确时，有1/3的射击效果能毁 伤敌人一门火炮，有1/6的射击效果能全部毁伤敌人火炮．
蒙特卡罗方法的关键步骤在于随机数的产生， 计算机产生的随机数都不是真正的随机数(由算 法确定的缘故)，如果伪随机数能够通过一系列 统计检验，我们也可以将其当作真正的随机数 使用。
rand('seed',0.1);
rand(1) ％每次运ra行nd程('s序tat产e',s生um的(1值00*是clo相ck同)*r的and);
E = P(A0) = P(j=0)P(A0∣j=0) + P(j=1)P(A0∣j=1)
= 1 0 1 1 0.25 2 22
P(A1) = P(j=0)P(A1∣j=0) + P(j=1)P(A1∣j=1)
= 10 11 1 2 23 6
P(A2) = P(j=0)P(A2∣j=0) + P(j=1)P(A2∣j=1)
非常见分布的随机数的产生
• 逆变换方法
由定理 1 ，要产生来自 F(x) 的随机数，只要先 产生来自U (0,1) 随机数 u ，然后计算 F 1(u) 即 可。具体步骤如下：
(1) 生成 (0,1)上 均匀分布的随机数U。
(2) 计算 X F -1(U ) ,则 X 为来自 F(x) 分布的随机数.
蒙特卡罗方法的基本思想很早以前就被人们所发现和 利用。早在17世纪，人们就知道用事件发生的“频率” 来决定事件的“概率”。19世纪人们用蒲丰投针的方法 来计算圆周率π，上世纪40年代电子计算机的出现，特别 是近年来高速电子计算机的出现，使得用数学方法在计算 机上大量、快速地模拟这样的试验成为可能。
第六章 M onte-Carlo 方法
10
1、离散型分布随机变量的直接抽样 对一个可以取两个值的随机变量x，如果它以几率p1取值x1, 而以几率p2取值x2。则：p2=(1-p1)。如果取(0,1)间一个随机数， 若满足： x < p 1 , 则取： x = x 1
第六章 Monte-Carlo 方法 第一节 Monte-carlo 方法概述 Monte-Carlo(蒙特卡罗)是摩纳哥闻名的赌城的名字，其本意具 有“随机”、“机遇”之意，从而Monte-Carlo方法又称为随机 抽样技巧或统计模拟方法(statistical simulation method )。 是利用随机数进行数值模拟的方法。 是由Metropolis在二次世界大战期间提出的,Nouman命名。在 Manhattan计划中，研究与原子弹有关的中子输运过程。 Monte Carlo方法是现代计算技术的最为杰出的成果之一。
由于试验次数不能太少，进行大量模拟就有很大的运算量， 从而只有在计算机出现和发展后，该方法才得到有效应用，所 以说，Monte-Carlo方法是和计算机紧密联系在一起的。
5
三. Monte-Carlo 方法的适用范围非常广泛
由于空间维数的多少对于Monte-Carlo方法的影响不大，且受问 题 的条件限制小，另外用该方法解决问题所编写的程序结构简 单，所以该方法已广泛应用在许多领域。 它可以解决一些典型的数学问题， 如多重积分的计算、线性代 数方程组、线性积分方程求解、齐次线性积分方程本征值的计 算、微分方程边值的计算等； 另外生物、 物理、材料、化学、经济、通讯等 科学方面许多 复杂问题用该方法来解决相对来说比较简单。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Monte Carlo 方法法§1 概述Monte Carlo 法不同于确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。
Monte Carlo 方法（MCM ），也称为统计试验方法，是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态。
它是用一系列随机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。
运用该近似方法所获得的问题的解in spirit 更接近于物理实验结果，而不是经典数值计算结果。
普遍认为我们当前所应用的MC 技术，其发展约可追溯至1944年，尽管在早些时候仍有许多未解决的实例。
MCM 的发展归功于核武器早期工作期间Los Alamos （美国国家实验室中子散射研究中心）的一批科学家。
Los Alamos 小组的基础工作刺激了一次巨大的学科文化的迸发，并鼓励了MCM 在各种问题中的应用[2]-[4]。
“Monte Carlo ”的名称取自于Monaco （摩纳哥）内以赌博娱乐而闻名的一座城市。
Monte Carlo 方法的应用有两种途径：仿真和取样。
仿真是指提供实际随机现象的数学上的模仿的方法。
一个典型的例子就是对中子进入反应堆屏障的运动进行仿真，用随机游动来模仿中子的锯齿形路径。
取样是指通过研究少量的随机的子集来演绎大量元素的特性的方法。
例如，)(x f 在b x a <<上的平均值可以通过间歇性随机选取的有限个数的点的平均值来进行估计。
这就是数值积分的Monte Carlo 方法。
MCM 已被成功地用于求解微分方程和积分方程，求解本征值，矩阵转置，以及尤其用于计算多重积分。
任何本质上属随机组员的过程或系统的仿真都需要一种产生或获得随机数的方法。
这种仿真的例子在中子随机碰撞，数值统计，队列模型，战略游戏，以及其它竞赛活动中都会出现。
Monte Carlo 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。
§2 随机数和随机变量的产生[5]－[10]全面的论述了产生随机数的各类方法。
其中较为普遍应用的产生随机数的方法是选取一个函数)(x g ，使其将整数变换为随机数。
以某种方法选取0x ，并按照)(1k k x g x =+产生下一个随机数。
最一般的方程)(x g 具有如下形式：m c ax x g mod )()(+= （1）其中=0x 初始值或种子（00>x ）=a 乘法器（0≥a ）=c 增值（0≥c ）=m 模数对于t 数位的二进制整数，其模数通常为t2。
例如，对于31位的计算机m 即可取1312-。
这里a x ,0和c 都是整数，且具有相同的取值范围0,,x m c m a m >>>。
所需的随机数序{}n x 便可由下式得m c ax x n n m od )(1+=+ （2） 该序列称为线性同余序列。
例如，若70===c a x 且10=m ，则该序列为7，6，9，0，7，6，9，0…… （3） 可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。
（3）式中序列周期长度为4。
当然，一个有用的序列必是具有相对较长周期的序列。
许多作者都用术语乘同余法和混合同余法分别指代0=c 和0≠c 时的线性同余法。
选取c a x ,,0和m 的法则可参见[6,10]。
这里我们只关心在区间)1,0(内服从均匀分布的随机数的产生。
用字符U 来表示这些数字，则由式（2）可得mx U n 1-=（4） 这样U 仅在数组{}m m m m /)1(,......,/2,/1,0-中取值。
（对于区间（0,1）内的随机数，一种快速检测其随机性的方法是看其均值是否为0.5。
其它检测方法可参见[3,6]。
）产生区间),(b a 内均匀分布的随机数X ，可用下式U a b a X )(-+= （5） 用计算机编码产生的随机数（利用式（2）和（4））并不是完全随机的；事实上，给定序列种子，序列的所有数字U 都是完全可预测的。
一些作者为强调这一点，将这种计算机产生的序列称为伪随机数。
但如果适当选取c a ,和m ，序列U 的随机性便足以通过一系列的统计检测。
它们相对于真随机数具有可快速产生、需要时可再生的优点，尤其对于程序调试。
Monte Carlo 程序中通常需要产生服从给定概率分布)(x F 的随机变量X 。
该步可用[6]，[13]－[15]中的几种方法加以实现，其中包括直接法和舍去法。
直接法（也称反演法或变换法），需要转换与随机变量X 相关的累积概率函数)()(x X prob x F ≤=（即：)(x F 为x X ≤的概率）。
1)(0≤≤x F 显然表明，通过产生(0,1)内均匀分布随机数U ，经转换我们可得服从)(x F 分布的随机样本X 。
为了得到这样的具有概率分布)(x F 的随机数X ，不妨设)(x F U =，即可得)(1U FX -= （6）其中X 具有分布函数)(x F 。
例如，若X 是均值为μ呈指数分布的随机变量，且 ∞<<-=-x ex F x 0,1)(/μ（7）在)(x F U =中解出X 可得)1ln(U X --=μ （8） 由于)1(U -本身就是区间(0,1)内的随机数，故可简写为U X ln μ-= （9）有时(6)式所需的反函数)(1x F-不存在或很难获得。
这种情况可用舍去法来处理。
令dxx dF x f )()(=为随机变量X 的概率密度函数。
令b x a >>时的0)(=x f ，且)(x f 上界为M （即：M x f ≤)(），如图1所示。
我们产生区间(0,1)内的两个随机数),(21U U ，则 11)(U a b a X -+= M U f 21= （10） 分别为在(a,b)和(0,M)内均匀分布的随机数。
若)(11X f f ≤ （11） 则1X 为X 的可选值，否则被舍去，然后再试新的一组),(21U U 。
如此运用舍去法，所有位于)(x f 以上的点都被舍去，而位于)(x f 上或以下的点都由11)(U a b a X -+=来产生1X 。
图1 舍去法产生概率密度函数为)(x f 的随机变量例1 设计一子程序使之产生0,1之间呈均匀分布的随机数U 。
用该程序产生随机变Θ，其概率分布由下式给定πθθθ<<-=0),cos 1(21)(T 解：生成U 的子程序如图2所示。
该子程序中，,0,21474836471221==-=c m 且1680775==a 。
应用种子数（如1234），主程序中每调用一次子程序，就会生成一个随机数U 。
种子数可取1到m 间的任一整数。
0001 C********************************************************** 0002 C PROGRAM FOR GENERA TING RANDOM V ARIABLES WITH 0003 C A GIVEN PROBABILITY DISTRIBUTION0004 C**********************************************************00050006 DOUBLE PRECISION ISEED 00070008 ISEED=1234.D0 0009 DO 10 I=1,1000010 CALL RANSOM(ISEED,R) 0011 THETA=ACOSD(1.0-2.0*R) 0012 WRITE(6,*)I,THETA 0013 10 CONTINUE 0014 STOP 0015 END0001 C**********************************************************0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN 0003 C THE INTERV AL (0,1)0004 C********************************************************** 00050006 SUBROUTINE RANDOM (ISEED,R) 0007 DOUBLE PRECISION ISDDE,DEL,A 0008 DATA DEL,A/2147483647.D0,16807.D0/ 00090010 ISDDE=DMOD(A*ISDDE,DEL) 0011 R=ISDDE/DEL 0012 RETURN 0013 END图2 例1的随机数生成器图2的子程序只是为了说明本章所介绍的一些概念。
大多数计算机都有生成随机数的子程序。
为了生成随机变量Θ，令 )cos 1(21)(Θ-=Θ=T U 则有 )21(cos )(11U U T-==Θ--据此，一系列具有给定分布的随机变量Θ便可由图2所示主程序中生成。
§3 误差计算Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。
因此，该解中包含了平均值附近的浮动量，而且不可能达到100％的置信度。
要计算Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。
我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计[13,16]。
设X 是随机变量。
则X 的期望值或均值x 定义为 ⎰∞∞-=dx x xf x )( （12）这里)(x f 是X 的概率密度分布函数。
如果从)(x f 中取些独立的随机样本N x x x ,...,,21，那么的x 估计值就表现为N 个样本值的均值。
∑==Nn nxNx11ˆ （13）x 是X 的真正的平均值，而xˆ只是x 的有着准确期望值的无偏估计。
虽然x ˆ的期望值等于x ，但x x≠ˆ。
因此，我们还需要x ˆ的值在x 附近的分布测度。
为了估计X 以及x ˆ在x 附近的的值的分布，我们需要引入X 的方差，其定义为X 与x差的平方的期望值，即 ⎰∞∞--=-==dx x f x x x x x Var )()()()(222σ（14）由2222)(x x x x x x +-=-，故有 ⎰⎰⎰∞∞-∞∞-∞∞-+-=dx x f x dx x xf x dx x f x x )()(2)()(222σ （15）或者222)(x x x -=σ （16）方差的平方根称为标准差，即2/122)()(x x x -=σ （17）标准差给出了x 在均值x 附近的分布测度，并由此给出了误差幅度的阶数。