Monte Carlo 法
§8.1 概述
Monte Carlo 法不同于前面几章所介绍的确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。Monte Carlo 方法（MCM ），也称为统计试验方法，是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态[1]。它是用一系列随机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。运用该近似方法所获得的问题的解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 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。
§8.2 随机数和随机变量的产生
[5]－[10]全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数)(x g ，使其将整数变换为随机数。以某种方法选取
0x ，并按照)(1k k x g x =+产生下一个随机数。最一般的方程)(x g 具有如下形式：
m c ax x g mod )()(+=
（8.1）
其中
=0x 初始值或种子（00>x ）
=a 乘法器（0≥a ）
=c 增值（0≥c ）
=m 模数
对于t 数位的二进制整数，其模数通常为t 2。例如，对于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+=+ （8.2）
该序列称为线性同余序列。例如，若70===c a x 且10=m ，则该序列为
7，6，9，0，7，6，9，0…… （8.3）
可以证明，同余序列总会进入一个循环套；也就是说，最终总会出现一个无休止重复的数字的循环。（8.3）式中序列周期长度为4。当然，一个有用的序列必是具有相对较长周期的序列。许多作者都用术语乘同余法和混合同余法分别指代
0=c 和0≠c 时的线性同余法。选取c a x ,,0和m 的法则可参见[6,10]。
这里我们只关心在区间)1,0(内服从均匀分布的随机数的产生。用字符U 来表示这些数字，则由式（8.2）可得
m
x U n 1
-= （8.4）
这样U 仅在数组{}m m m m /)1(,......,/2,/1,0-中取值。（对于区间（0,1）内的随机数，一种快速检测其随机性的方法是看其均值是否为0.5。其它检测方法可参见[3,6]。）产生区间),(b a 内均匀分布的随机数X ，可用下式
U a b a X )(-+= （8.5）
用计算机编码产生的随机数（利用式（8.2）和（8.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 F X -= （8.6）
其中X 具有分布函数)(x F 。例如，若X 是均值为μ呈指数分布的随机变量，且 ∞<<-=-x e x F x 0,1)(/μ （8.7）
在)(x F U =中解出X 可得
)1ln(U X --=μ （8.8）
由于)1(U -本身就是区间(0,1)内的随机数，故可简写为
U X ln μ-= （8.9）
有时(8.6)式所需的反函数)(1x F -不存在或很难获得。这种情况可用舍去法来
处理。令dx
x dF x f )
()(=
为随机变量X 的概率密度函数。令b x a >>时的0)(=x f ，且)(x f 上界为M （即：M x f ≤)(），如图8.1所示。我们产生区间(0,1)内的两个随机数),(21U U ，则
1
1)(U a b a X -+=
M U f 21=
（8.10）
分别为在(a,b)和(0,M)内均匀分布的随机数。若
)(11X f f ≤ （8.11）
则1X 为X 的可选值，否则被舍去，然后再试新的一组),(21U U 。如此运用舍去法，所有位于)(x f 以上的点都被舍去，而位于)(x f 上或以下的点都由
11)(U a b a X -+=来产生1X 。
图8.1 舍去法产生概率密度函数为)(x f 的随机变量
例8.1 设计一子程序使之产生0,1之间呈均匀分布的随机数U 。用该程序产生随机变Θ，其概率分布由下式给定
πθθθ<<-=0),cos 1(21
)(T
解：生成U 的子程序如图8.2所示。该子程序中，,0,21474836471221==-=c m 且1680775==a 。应用种子数（如1234），主程序中每调用一次子程序，就会生成一个随机数U 。种子数可取1到m 间的任一整数。
0001
C**********************************************************
0002 C PROGRAM FOR GENERATING RANDOM V ARIABLES
WITH
0003 C A GIVEN PROBABILITY DISTRIBUTION 0004
C********************************************************** 0005
0006 DOUBLE PRECISION ISEED 0007
0008 ISEED=1234.D0 0009 DO 10 I=1,100
0010 CALL RANSOM(ISEED,R) 0011 THETA=ACOSD(1.0-2.0*R) 0012 WRITE(6,*)I,THETA 0013 10 CONTINUE 0014 STOP 0015 END
0001
C**********************************************************
0002 C SUBROUTINE FOR GENERATING RANDOM NUMBERS IN
0003 C THE INTERV AL (0,1) 0004
C********************************************************** 0005
0006 SUBROUTINE RANDOM (ISEED,R) 0007 DOUBLE PRECISION ISDDE,DEL,A 0008 DATA DEL,A/2147483647.D0,16807.D0/ 0009
0010 ISDDE=DMOD(A*ISDDE,DEL) 0011 R=ISDDE/DEL 0012 RETURN 0013 END
图8.2 例8.1的随机数生成器
图8.2的子程序只是为了说明本章所介绍的一些概念。大多数计算机都有生成随机数的子程序。
为了生成随机变量Θ，令
)cos 1(21
)(Θ-=Θ=T U
则有 )21(cos )(11U U T -==Θ--
据此，一系列具有给定分布的随机变量Θ便可由图8.2所示主程序中生成。
§8.3 误差计算
Monte Carlo 程序给出的解按大量的检测统计都达到了平均值。因此，该解中包含了平均值附近的浮动量，而且不可能达到100％的置信度。要计算Monte Carlo 算法的统计偏差，就必须采用与统计变量相关的各种统计方法。我们只简要介绍期望值和方差的概念，并利用中心极限定理来获得误差估计[13,16]。
设X 是随机变量。则X 的期望值或均值x 定义为
?∞
∞-=dx x xf x )(
（8.12）
这里)(x f 是X 的概率密度分布函数。如果从)(x f 中取些独立的随机样本
N x x x ,...,,21，那么的x 估计值就表现为N 个样本值的均值。
∑==N
n n
x
N
x
1
1?
（8.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 )()()()(22
2σ
（8.14）
由2222)(x x x x x x +-=-，故有
???∞
∞
-∞
∞
-∞
∞
-+-=dx x f x dx x xf x dx x f x x )()(2)()(222σ
（8.15）
或者 222)(x x x -=σ （8.16）
方差的平方根称为标准差，即 2/122)()(x x x -=σ （8.17）
标准差给出了x 在均值x 附近的分布测度，并由此给出了误差幅度的阶数。x
?的标准差与x 的标准差的关系表示为
N
x x
)
()?(σσ=
（8.18）
该式表明，如果用根据(8.13)式由N 个n x 值构成的x
?来估计x ，那么结果中x ?在x 附近的扩散范围便与)(x σ成比例，且随着样本数N 的增加而降低。
为了估计x
?的扩散范围，我们定义样本方差为 ∑=--=N
n n x x N S 1
22
)?(11 （8.19）
由此式还可看出2S 的期望值等于)(2x σ。因此样本方差是)(2x σ的无偏估计。将（8.19）式中平方项乘出来，便可得样本标准差为
2
/11222
/1?11??
????-?
??
??-=∑=N n n x x N N N S
（8.20）
当N 较大时，系数)1/(-N N 可设为1。
作为获得中心极限定理的一种方法，概率论的一个基本解可考虑二次项函数 M N M q p M N M N M B --=)!
(!!
)(
（8.21）
该式表明N 次独立随机试验中有M 次成功的概率。(8.21)中，p 是一次试验中成
功的概率，且p q -=1。当M 和M N -都较大时，便可用Stirling 公式
n e n n n n π2~!- （8.22）
因此(8.21)式可近似为正态分布[17]
??
?
???--=≈)?(2)?(exp 2)?(1
)?()(2
2x x x x x f M B σπσ (8.23)
其中p N x =且Npq x =)?(σ。因此当∞→N 时，中心极限定理表明，描述由N 点Monte Carlo 算法获得的x
?的分布的概率密度函数是（8.23）式所示的正态分布函数)(x f 。也就是说，大量随机变量的集合趋于呈正态分布。将（8.18）式代入（8.23）式可得
??
????--=)(2)?(exp )(1
2)?(2
2x x x N x N x
f σσπ （8.24）
正态（高斯）分布在工程，物理以及统计学的各类问题中都非常有用。高斯模型的显著特性源于中心极限定理。因此，高斯模型经常用于其影响程度取决于由许多不规则的、浮动的元素构成的集合的情况。例8.2中我们给出了根据中心极限定理产生高斯随机变量的算法。
由于样本数N 是有限的，所以Monte Carlo 算法不可能达到绝对的确定性。
故而在x 附近估计出某一范围或区间以确保我们估计的x
?落入该范围内。假设要得到x
?位于δ-x 和δ+x 之间的概率。由定义 {}?+-=+<<-δδ
δδx x x d x
f x x
x ob ?)?(?Pr （8.25） 令()()
x N x x
σλ/2?-=
可得
{}(
)
()
?
-=+<<-σδλλπ
δδ/2/0
2
2
?Pr N d e x x
x ob
()????
??=x N erf σδ2/ （8.26a ） 或
ασσαα-=???
?
??+≤≤-1?Pr 2/2/N z x x N z x ob
（8.26b ）
其中)(x erf 是误差函数，且2/αz 是标准正态差的上2/α分位点。对（8.26）式可做如下解释：如果用来呈现独立随机观测值并构建相关随机区间δ±x 的Monte
Carlo 程序以较大的N 值反复运行，则这些随机区间的()????
??x N erf σδ2分位点将
近似等于x ?。随机区间δ±x 称为置信区间，()???
? ??x N erf σδ2称为置信度。大多数
的Monte Carlo 算法都使用误差()N x /σδ=，这表明x
?是在实际均值x 的标准差范围内的。由（8.26）式可得，样本均值x
?位于区间()N x x /?σ±内的概率是0.6826或68.3％。若要求更高的置信度，可用两个或三个标准差的值。如
???
??=????
??+<<-,
997.0,954.0,
6826.0)(?)(Pr N x M x x N x M x ob σσ
3
21
===M M M （8.27）
其中M 是标准差的个数。
在（8.26）和（8.27）式中均假设总体标准差σ为已知。由于这种情况很少出现，故必须用由（8.20）式算得的样本标准差S 来估计σ的值，从而使学生氏t －分布取代正态分布。我们知道当N 很大，比如30>N 时，t －分布近似趋于正态分布。（8.26）式等价于
ααα-=??????
+≤≤---1?Pr 1;2/1;2/N St x x N St x ob N N （8.28）
其中1;2/-N t α为自由度是1-N 的学生氏t －分布的上2/α分位点；其值在任何标准统计学课本中都有表列可查。这样置信区间的上、下限便可由下式给出 上限＝N
St x N 1
;2/-+
α
（8.29）
下
限
＝
N
St x N 1
;2/--
α
（8.30）
Monte Carlo 算法中关于误差估计的进一步讨论，可参见[18,19]。
例8.2 利用中心极限定理产生一高斯（正态）分布的随机变量X 。根据中心极限定理，大量均值附近的独立随机变量的总和，无论其个体变量的分布如何，
总趋向于一种高斯分布。也就是说，对于任意随机数N i Y i ,....2,1,=，均值为Y ，方差为)(Y Var ，
()
Y Var N Y
N Y
Z N
i i
-=∑=1
（8.31）
渐渐与N 合并为零均值、单位标准差的正态分布。若i Y 是均匀分布的随机变量（即i i U Y =），则12/1)(,2/1==Y Var Y ，故而
12
/2
/1
N N U
Z N
i i
-=∑=
（8.32） 且变量
μσ+=Z X （8.33）
近似等于均值为μ、方差为2σ的正态变量。N 小于3时近似为我们所熟知的钟形高斯分布。为简化计算，通常实际中设12=N ，因为这样可消除（8.32）式中的平方根项。然而N 的这种取值截掉了σ6±边界处的分布，且无法产生超过σ3的值。对于分布曲线的两端比较重要的仿真，就必须用其它方案来产生高斯分布（参见[20]－[22]）。 这样，要产生一个均值为μ、标准差为σ的高斯随机变量X ，就要遵循以下步骤：
（1） 生成12个均匀分布的随机数1221,....,,U U U （2） 求得612
1-=∑=i i U Z
（3） 令μσ+=Z X
§8.4 数值积分
对于一维积分，现已有一些求积公式，如3.10节中所述。而对多重积分的公式则相对较少。Monte Carlo 技术对此类多重积分的重要性至少存在两方面的原因。多重积分的求积公式非常复杂，而多重积分的MCM 几乎保持不变。Monte Carlo 积分的收敛性与维数无关，但对求积公式并非如此。人们已经发现，积分
的统计方法是计算天线问题尤其是那些与较大结构相关的问题中的二维或三维积分的很有效的方法[23]。这里将讨论两类Monte Carlo 积分方法，即简易MCM 和具有对立变量的MCM 。其它类型，如成功－失败法和控制变量法，可参见[24]－[26]。本节还将简要介绍MCM 在不规范积分中的应用。
§8.4.1 简易Monte Carlo 积分
设要计算积分 ?=R
f I
（8.34）
其中R 是n 维空间。令()n X X X X ,....,21=是在R 内均匀分布的随机变量。则
)(X f 也是随机变量，其均值由下式给出[27,28]
R
I
f R
X f R
=
=?
1)( (8.35) 方差为
()()2
2
11???
? ?
?-=?
?
R
R
f R
f R
X f Var （8.36）
其中 ?=R
dX R
（8.37）
如果取X 的N 个独立样本，即N X X X ,....,,21，它们与X 具有相同的分布，且构成平均项
()()()()∑==+++N
i i
N X f N
N X f X f X f 1
211
....
（8.38）
我们便会期望该平均项接近于)(X f 的均值。这样，由（8.35）和（8.38）式，即可得
()∑==
N
i i
X f N
R
I 1
（8.39）
Monte Carlo 公式可用于有限区域R 上的任何积分。为了举例说明，现将（8.39）式应用于一维和二维积分中。 对于一维积分，设
?=b
a
dx x f I )(
（8.40）
由（8.39）式可得
()∑=-=N
i i X f N a b I 1
（8.41）
其中i X 是区间()b a ,内随机变量，即
()10,<<-+=U U a b a X i （8.42）
对于二维积分，有
()
??=b a
d c
dX dX X X
f I 2121
,
（8.43）
则相应的Monte Carlo 公式为
()()
∑=--=N
i i i X X f N c d a b I 1
21,)( （8.44）
其中
1
0,)(10,)(2
2
2111<<-+=<<-+=U U c d c X U U a b a X i
i
（8.45） （8.39）式中无偏估计值I 收敛的很慢，这是由于估计值方差的级数是N /1。一种改良的方法，即控制变量法，可通过减小估计值方差来提高其准确性和收敛性。
§8.4.2 具有对立变量的Monte Carlo 积分方法 术语对立变量[29,30]用来描述任一系列估计值，它们能互相抵消掉彼此的方差。方便起见，我们假设积分范围为区间(0,1)。设要得到下面单重积分的估计值 ?=1
0)(dU U g I
（8.46） 我们期望表达式
[])1()(2
1
U g U g -+与)(U g 相比具有更小的方差。如果)(U g 太小，那反过来)1(U g -就很可能太大。因此，我们定义估计值
()()[]∑=-+=
N
i i i U g U g N
I 112
1
1
其中i U 是0和1之间的随机数。此估计值方差的级数是4
1
N ，这是在（8.39）式基础上的一个极大的进步。对于两维积分 (
)
??=101
2121,dU dU U U g I
（8.48）
则相应的估计值为
()()()()[]
∑=--+-+-+=
N
i i i i i i i i i U U g U U g U U g U U g N
I 1212121211,1,11,,4
1
1
（8.49）
根据相似性，可将该方法延拓至更高阶的积分。对于不是(0,1)的区间，可应用诸如（8.41）式到（8.45）式的转换。例如，
()()()??-=1
dU U g a b dx x f b a
()()[]∑=-+-≈N i i i U g U g N a b 112
1
（8.50）
其中)()(X f U g =，且()U a b a X -+=。由（8.47）和（8.49）式可得，当积分维数增加时，每一维用来在简易Monte Carlo 方法上提高效率的对立变量的最小值也会增加。这样使得简易Monte Carlo 方法在多维运算中更具优越性。
§8.4.3 不规范积分 积分式
?∞
=0)(dx x g I
（8.51）
可用Monte Carlo 仿真法进行估计[31]。对于概率密度为)(x f 的随机变量X ，其中)(x f 在区间(0,∞)上的积分等于1，则有
()()
??
∞∞
=00
)(dx x g dx x f x g （8.52）
因此，为计算（8.51）式中的I 值，首先要得到N 个服从同一在区间(0,∞)上的积分等于1的概率密度函数)(x f 的独立随机变量。其样本均值
()()
()
∑
==
N
i i i x f x g N
x g 1
1
便是I 的估计值。
例8.3 用Monte Carlo 方法计算积分 ?
?
=1020
cos πφαρφρρd d e I j
解：该积分式表示振幅和相位分别服从某一分布的圆形孔径天线的辐射状况。之所以选择该式，是因为它至少是辐射积分的一部分。且其解的闭合形式亦为可得，由此便可得Monte Carlo 解的精确性。其闭合解为 ()()
α
απα12J I =
其中()α1J 是一阶Bessel 函数。
图8.3给出由（8.44）和（8.45）式计算该积分的一个简单的程序，其中
0,1,0===c b a 以及π2=d 。该程序在Vax 11/780 中调用子程序RANDU 来产
生随机数1U 和2U 。对于不同的N 值，简易和对立变量Monte Carlo 方法都可用于计算辐射积分。表8.1中对5=α时的结果进行了精确的比较。在应用（8.49）式时，用到了下面的对应式：
()()111221111,,U a b X b U X U X U --=-≡-≡≡ ()()22211U c d X d U --=-≡-
表8.1 例8.3辐射积分的MC 方法积分结果
N 简易MCM 对立变量MCM 500 -0.2892-j0.0742
-0.2887-j0.0585
1000 -0.5737+j0.0808
-0.4982-j0.0080
2000 -0.4922-j0.0040
-0.4682-j0.0082
4000 -0.3999-j0.0345
-0.4216-j0.0323
6000 -0.3608-j0.0270
-0.3787-j0.0440
8000 -0.4327-j0.0378
-0.4139-j0.0241
10,000 -0.4229-j0.0237
-0.4121-j0.0240 精确解：-0.4116+j0
§8.5 位势问题的解
位势理论与布朗运动（或随机游动）的关系是由Kakutani在1944年首次提出的[32]。自此，所谓的概率位势理论便在诸如热传导[33]－[38]、静电学[39]－[46]以及电力工程等许多学科领域得到广泛应用。不同方程式的概率解或MC解的一个基本概念就是随机游动。不同类型的随机游动应用不同的Monte Carlo方法。最常见的类型是固定随机游动和自动定位随机游动。其它不常见类型有迁离法、缩减边界法、内切形法以及表面密度法。
0001
C***************************************************************** 0002 C INTEGRATION USING CRUDE MONTE CARLO
0003 C AND ANTITHETIC METHODS
0004 C ONLY FEW LINES NEED BE CHANGED TO USE THIS PROGRAM
0005 C FOR ANY MULTI—DIMENSIONAL INTEGRATION
0006
C**************************************************************** 0007
0008 DATA IS1,IS2,IS3,IS4/1234,5678,9012,3456/
0009 DATA A,B,C/0.0,1.0,0.0/
0010
0011 C SPECIFY THE INTEGRAND
0012 F(RHO,PHI)=RHO*CEXP(J*ALPHA*RHO*COS(PHI))
0013
0014 J=(0.0,1.0)
0015 ALPHA=5.0
0016 PIE=3.1415927
0017 D=2.0*PIE
0018 DO 30 NRUN=500,10000,500
0019 SUM1=(0.0,0.0)
0020 SUM2=(0.0,0.0)
0021 DO 10 I=1,NRUN
0022 CALL RANDU(IS1,IS2,U1)
0023 CALL RANDU(IS3,IS4,U2)
0024 X1=A+(B-A)*U1
0025 X2=C+(D-C)*U2
0026 X3=(B-A)*(1.0-U1)
0027 X4=(D-C)*(1.0-U2)
0028 SUM1=SUM1+F(X1,X2)
0029 SUM2=SUM2+F(X1,X2)+F(X1,X4)+F(X3,X2)+F(X3,X4)
0030 10 CONTINUE
0031 AREA1=(B-A)*(D-C)*SUM1/FLOAT(NRUN)
0032 AREA2=(B-A)*(D-C)*SUM2/(4.0*FLOAT(NRUN)) 0033 PRINT *,NRUN,AREA1,AREA2 0034 WRITE(6,*) NRUN,AREA1,AREA2 0035 WRITE(6,20) NRUN,AREA1,AREA2 0036
20
FORMAT(2X,’NRUN=’,I5,3X,’AREA1=’,F12.6,3X,F12.6,’AREA2=’, 0037 1 F12.6,3X,F12.6,/) 0038 30 CONTINUE 0039 STOP 0040 END
图8.3 例8.3中用Monte Carlo 方法计算二维积分的程序 §8.5.1 固定随机游动
为具体起见，假设用固定随机游动的MCM 解拉普拉斯方程
02=?V （在区域R 内） （8.54a ）
满足Dirichlet 边界条件
P V V =
（在边界
B 上）
（8.54b ）
首先将R 分成网状结构，并用其有限差当量替代2?。（8.54a ）在二维R 内的有限差表达式可由（3.26）式给出，即
()()()()()?-+?++?-+?+=-+-+y x V p y x V p y x V p y x V p y x V y y x x ,,,,, （8.55a ） 其中
4
1
=
===-+-+y y x x p p p p （8.55b ）
（8.55）式中，假设网络的一个方格尺寸是?，如图8.4所示。下面给出该方程的概率解释。若随机游动粒子在某一瞬时处于()y x ,位置处，它从此点向),(),,(y x y x ?-?+),(?+y x 及()?-y x ,移动的概率分别是+-+y x x p p p ,,和-y p 。决定粒子移动方式的方法是产生一个随机数10,<
()()
()()
()()()()
?-→?+→?-→?+→y x y x y x y x y x y x y x y x ,,,,,,,,
1
75.075.05.05.025.025
.00<<<<<<<
（8.56）
如果不用方格而用矩形格，则有-+=x x p p 且-+=y y p p ，但y x p p ≠。在用立方体晶格表示的三维问题中，6
1
======-+-+-+z z y y x x p p p p p p 。两种情况中都依据概率对区间10<
图8.4 固定随机游动示意图
为了计算()y x ,处的位势，让一随机游动粒子在该点出发开始游动。粒子便开始从一点到另一点在网格中蜿蜒前行直到到达边界。此时，粒子终止游动，记下边界点的指定位势P V 。令第一个粒子游动结束时P V 的值记为()1P V ，如图8.4所示。然后将第二个粒子从()y x ,释放，使其游动直到到达边界点，游动终止且该点P V 的值相应的记为()2P V 。依次第三个，第四个，K K ，第N 个粒子从()y x ,释放重复此过程，并记下相应的指定位势()()4,3P P V V ，K K ()N V P 。根据Kakutani [32]，()()2,1P P V V ，K K ()N V P 的期望值是Dirichlet 问题在()y x ,的解，即
()()∑==
N
i P i V N
y x V 1
1,
（8.57）
其中N 较大，为游动总次数。其收敛速度随N 而改变，所以需要保证许多随机游动都有精确的结果。 若要解Poisson 方程
()y x g V ,2-=?
在
R 内
（8.58a ） 且V 满足条件
P V V =
在
B 上
（8.58b ）
则式（8.55）中的有限差分表示为
()()()y x V p y x V p y x V x x ,,,?-+?+=-+
()()4
,,2g
y x V p y x V p y y ?+?-+?++-+
（8.59）
其中概率仍为式（8.55b ）所示。（8.59）式的概率解释也近似于（8.55）式。但随机游动的每一步都要记录下（8.59）式中的4/2g ?项。如果第i 次随机游动从
()y x ,出发至到达边界需要i m 步，则记录下
()()∑=?+
i
m j j
j
P y x g i V 1
2,4
（8.60）
由此可得，()y x V ,的Monte Carlo 解为
()()()∑∑∑=-==???????+=N i m j j j N i P i y x g N i V N y x V 11
1
21,41,
（8.61）
对刚才所描述MCM 的一个有趣的分析就是游走醉鬼问题[15,35]。我们将随机游动粒子当作一个“醉鬼”，将网格方块当作“城市街区”，结点当作“十字路口”，边界B 当作“城市边界”，而且把边界B 上的结点当作“警察”。尽管醉鬼想步行走回家，但他却醉的一塌糊涂以至于在整个城市中随机游走。警察的工作是在城市边界上一看见醉鬼便将其抓获，并勒令醉鬼交付罚金P V 。那么醉鬼要支付的预期罚金将会是多少呢？其答案便是（8.57）式。 在电介质边界上，指定边界条件为n n D D 21=。
()y x ,
x
?x ?x ?x ?x ?x ?x ?x ?x ? x x x x x x x x x x x x x x
X X X X X X X X X X X X X X X X X X X X
[1].R.Hersch and R.J.Griego,“Brownian motion and potential theory,”
Sci.Amer.,Mar.1969,pp.67-74.
当科学家们使用计算机来试图预测复杂的趋势和事件时, 他们通常应用一类需要长串的随机数的复杂计算。设计这种用来预测复杂趋势和事件的数字模型越来越依赖于一种称为蒙特卡罗模似的统计手段, 而这种模拟进一步又要取决于可靠的无穷尽的随机数目来源。 蒙特卡罗模拟因摩纳哥著名的赌场而得名。它能够帮助人们从数学上表述物理、化学、工程、经济学以及环境动力学中一些非常复杂的相互作用。数学家们称这种表述为“模式”, 而当一种模式足够精确时, 他能产生与实际操作中对同一条件相同的反应。但蒙特卡罗模拟有一个危险的缺陷: 如果必须输入一个模式中的随机数并不像设想的那样是随机数, 而却构成一些微妙的非随机模式, 那么整个的模拟(及其预测结果)都可能是错的。 最近, 由美国佐治亚大学的费伦博格博士作出的一分报告证明了最普遍用以产生随机数串 的计算机程序中有5个在用于一个简单的模拟磁性晶体中原子行为的数学模型时出现错误。科学家们发现, 出现这些错误的根源在于这5个程序产生的数串其实并不随机, 它们实际上隐藏了一些相互关系和样式, 这一点只是在这种微小的非随机性歪曲了晶体模型的已知特 性时才表露出来。贝尔实验室的里德博士告诫人们记住伟大的诺伊曼的忠告:“任何人如果相信计算机能够产生出真正的随机的数序组都是疯子。” 蒙特卡罗方法（MC） 蒙特卡罗（Monte Carlo）方法： 蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是我们采用该方法的原因。 蒙特卡罗方法的基本原理及思想如下： 当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。 蒙特卡罗解题三个主要步骤： 构造或描述概率过程： 对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 实现从已知概率分布抽样： 构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样
蒙特卡洛习题 1.利用蒙特卡洛计算数值积分 () ()() 1280ln 1tan x x x xe dx +++? clear all ;clc;close all ; n=1000; count=0; x=0:0.01:1; y=log((1+x).^2+(tan(x).^8)+x.*exp(x)); plot(x,y,'linewidth',2) hold on for i=1:n x1=rand; y1=rand*y(end); plot(x1,y1,'g*') pause(0.01) if y1 2.分别用理论计算和计算机模拟计算，求连续掷两颗骰子，点数之和大于6且第一次掷出的点数大于第二次掷出点数的概率。 clear all;clc;close all; count=0; n=100000; for i=1:n x=floor(rand*6+1); y=ceil(rand*6); if x+y>6&&x>y count=count+1; end end P=count/n 3. clear all;clc;close all; count=0; n=2000; ezplot('x^2/9+y^2/36=1'); hold on ezplot('x^2/36+y^2=1'); hold on ezplot('(x-2)^2+(y+1)^2=9') for i=1:n x=rand*12-6; y=rand*12-6; plot(x,y,'gh','linewidth',2) pause(0.01) if x^2/9+y^2/36<1&&x^2/36+y^2<1&&(x-2)^2+(y+1)^2<9 2.3.2 蒙特卡洛法的基本原理 蒙特卡洛模型的基本原理是模拟单个光子的传输过程，本质上是一系列随机作用和随机过程的计算机模拟，如光子吸收、散射、传输路径、步长等。光子从发射到进入组织再到从组织中逸出要历经许多过程，以单个光子为例，首先是光子发射，即单个光子垂直入射到组织表面，光子质量W 被初始化为1，当组织与周围介质折射率不同时，在入射界面处要考虑镜面反射（界面不光滑时考虑漫折射），其反射比设为RSP ，因此进入介质的能量为1-RSP ，这部分能量就是接下来要进行蒙特卡洛模拟的部分。进入组织后光子继续运动，首先要确定其运动步长s ，根据光子的运动步长和运动方向，可以得到光子与组织发生相互作用的坐标位置，并以此坐标为起点开始下一运动步长的模拟。光子在与组织发生相互作用时有（μ a/μt）W 的能量被吸收，剩余部分能量的光子被散射，并继续重复上述过程，直到光子运动到边界处，此时，它有可能被返回到组织内部或者透过组织进入到周围介质。如果光子被反射，那么它将继续传播，即重复上述运动；如果光子穿透组织，根据其穿透的是前表面还是后表面，则相应被记入透射量和反射量。 由于蒙特卡洛模型的精确性是建立在大量模拟的基础上，因此这一方法耗时长，这与光谱技术的实时特性相矛盾。“查表法”的提出为这一问题提供了一种很好的解决途径，查表法的基本思想在于事先将一系列组织光学特性所对应的模拟结果存储到一个表格中，这样在对每一个光子进行模拟时，能够从这一表格中直接提取最终的模拟结果，从而节省了大量的模拟时间。 对于组织光子传输蒙特卡洛模型的研究已经开展了很多年，目前学术界广为接受和采用的是美国圣路易斯华盛顿大学华人教授Lihong Wang所提出的模型[1]，此模型是前向模型，即在已知组织吸收和散射特性的前提下对光子在组织中的传输分布进行模拟；美国杜克大学助理教授Gregory Palmer等在前向模型的基础上开发出了所谓的后向模型[2]，这一模型是在已知光谱反射特性的基础上，通过多次随机假定光学特性并调用前向模型进行光谱拟合，从而筛选出与实际测量结果最为匹配的一组假定数据作为组织的光学特性参数。后向模型的提出使得蒙特卡洛模型能够从真正意义上对组织的光学参数进行检测，并定量得出组织的各组分参数。目前蒙特卡洛模型已被广泛用于多种肿瘤的离体及临床在体研究，并取得了令人满意的结果，最终应用于临床检测的相关仪器也已得到开发，并预计将在未来的十几年甚至是十年之内推向临床应用。 当然目前关于这一模型仍有一定的发展提升空间，难点主要集中于如何进一步提高其精确性，这主要体现在两个方面：（1）如何进一步优化模型来提高精确性，目前这一模型对于仿体吸收散射特性的提取检测已经能够达到10%以内的误差精度，但最近的研究发现，将这一模型应用于仿体荧光检测时，其精确性仍有较大提升空间[3]。仿体荧光检测主要是为了研究模型提取固有荧光的能力，由于吸收和散射的存在，我们所检测的荧光并不是荧光物质本身的固有荧光，其光谱形状和强度均受到一定程度的改变，模型通过反射信号首先提取仿体的吸收和散射特性，进而用于对荧光信号进行矫正从而得到固有荧光光谱。研究发现，蒙特卡洛模型能够对荧光光谱形状进行良好恢复，但对于荧光光强的恢复其精确度仍有待提高。（2）如何提高用于人体组织检测的精确性，人体组织的情况往往是极为复杂的，这就需要开发精确的光子蒙特卡洛多层介质传输模型。目前关于这方面的研究已经取得一定的成果[1]，但仍需要开展更多的工作。 参考文献： [1] Wang L，Jacques SL，Zheng L. MCMLMonte Carlo Modeling of Light Transport in Multi-layered Tissues[J]. Comput Methods Programs Biomed，1995，47（2）：131-146. [2] Palmer GM，Ramanujam N. Monte Carlobased Inverse Model for Calculating Tissue Optical 浅析蒙特卡洛方法原理及应用 于希明 （英才学院1236103班测控技术与仪器专业6120110304） 摘要：本文概述了蒙特卡洛方法产生的历史及基本原理，介绍了蒙特卡洛方法的最初应用——蒲丰投针问题求圆周率，并介绍了蒙特卡洛方法在数学及生活中的一些简单应用，最后总结了蒙特卡洛方法的特点。 关键词：蒙特卡洛方法蒲丰投针生活应用 蒙特卡洛方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。它是以概率统计理论为基础, 依据大数定律( 样本均值代替总体均值) , 利用电子计算机数字模拟技术, 解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。蒙特卡洛方法在金融工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。 一、蒙特卡洛方法的产生及原理 蒙特卡洛方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡洛方法就已经存在。1777年，法国数学家蒲丰（Georges Louis Leclere de Buffon，1707—1788）提出用投针实验的方法求圆周率π。这被认为是蒙特卡洛方法的起源。 其基本原理如下：由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡洛法正是基于此思路进行分析的。 设有统计独立的随机变量Xi(i=1，2，3，…，k)，其对应的概率密度函数分别为fx1，fx2，…，fxk，功能函数式为Z=g(x1，x2，…，xk)。首先根据各随机变量的相应分布，产生N组随机数x1，x2，…，xk值，计算功能函数值Zi=g(x1，x2，…，xk)(i=1，2，…，N)，若其中有L组随机数对应的功能函数值Zi≤0，则当N→∞时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。 二、蒲丰投针问题 作为蒙特卡洛方法的最初应用, 是解决蒲丰投针问题。1777 年, 法国数学家蒲丰提出利用投针实验求解圆周率的问题。设平面上等距离( 如为2a) 画有一些平行线, 将一根长度为2l( l< a) 的针任意投掷到平面上, 针与任一平行线相交的频率为p 。针的位置可以用针的中心坐标x 和针与平行线的夹角θ来决定。任意方向投针, 便意味着x与θ可以任意取一值, 只是0≤x ≤a, 0≤θ≤π。那么, 投针与任意平行线相交的条件为x ≤ l sinθ。相交频率p 便可用下式求 二、蒙特卡洛模拟原理及步骤 （一）蒙特卡洛模拟原理：经济生活中存在大量的不确定与风险问题，很多确定性问题实际上是不确定与风险型问题的特例与简化，财务管理、管理会计中同样也存在大量的不确定与风险型问题，由于该问题比较复杂，一般教材对此问题涉及较少，但利用蒙特卡洛模拟可以揭示不确定与风险型问题的统计规律，还原一个真实的经济与管理客观面貌。 与常用确定性的数值计算方法不同，蒙特卡洛模拟是用来解决工程和经济中的非确定性问题，通过成千上万次的模拟，涵盖相应的可能概率分布空间，从而获得一定概率下的不同数据和频度分布，通过对大量样本值的统计分析，得到满足一定精度的结果，因此蒙特卡洛模拟是进行不确定与风险型问题的有力武器。 1、由于蒙特卡洛模拟是以实验为基础的，因此可以成为财务人员进行风险分析的“实验库”，获得大量有关财务风险等方面的信息，弥补确定型分析手段的不足，避免对不确定与风险决策问题的误导； 2、财务管理、管理会计中存在大量的不确定与风险型问题，目前大多数教材很少涉及这类问题，通过蒙特卡洛模拟，可以对其进行有效分析，解决常用决策方法所无法解决的难题，更加全面深入地分析不确定与风险型问题。 （二）蒙特卡洛模拟步骤以概率型量本利分析为例，蒙特卡洛模拟的分析步骤如下： 1、分析评价参数的特征，如企业经营中的销售数量、销售价格、产品生产的变动成本以及固定成本等，并根据历史资料或专家意见，确定随机变量的某些统计参数； 2、按照一定的参数分布规律，在计算机上产生随机数，如利用EXCEL提供的RAND函数，模拟量本利分析的概率分布，并利用VLOOKUP寻找对应概率分布下的销售数量、销售价格、产品生产的变动成本以及固定成本等参数； 3、建立管理会计的数学模型，对于概率型量本利分析有如下关系式，产品利润=产品销售数量×（产品单位销售价格-单位变动成本）-固定成本，这里需要说明的是以上分析参数不是确定型的，是依据某些概率分布存在的； 4、通过足够数量的计算机仿真，如文章利用RAND、VLOOKUP等函数进行30000次的模拟，得到30000组不同概率分布的各参数的排列与组合，由于模拟的数量比较大，所取得的实验数据具有一定的规律性； 5、根据计算机仿真的参数样本值，利用函数MAX、MIN、A VERAGE等，求出概率型量本利分析评价需要的指标值，通过对大量的评价指标值的样本分析，得到量本利分析中的利润点可能的概率分布，从而掌握企业经营与财务中的风险，为财务决策提供重要的参考。三、概率型量本利分析与比较 （一）期望值分析方法假设某企业为生产与销售单一产品的企业，经过全面分析与研究，预计未来年度的单位销售价格、销售数量、单位变动成本和固定成本的估计值及相应的概率如表1，其中销售数量单位为件，其余反映价值的指标单位为元，试计算该企业的生产利润。表1概率型量本利分析参数 项目概率数值 单位销售价格0.3 40 0.4 43 0.3 45 单位变动成本0.4 16 0.2 18 0.4 20 固定成本0.6 28000 0.4 30000 计算材料学之蒙特卡洛方法 一、计算材料学要紧内容 计算材料学涉及材料的各个方面，如不同层次的结构、各种性能等等，因此，有专门多相应的计算方法。在进行材料计算时，首先要依照所要计算的对象、条件、要求等因素选择适当的方法。要想做好选择，必须了解材料计算方法的分类。目前，要紧有两种分类方法：一是按理论模型和方法分类，二是按材料计算的特征空间尺寸(Characteristic space scale)分类。材料的性能在专门大程度上取决于材料的微结构，材料的用途不同，决定其性能的微结构尺度会有专门大的差不。例如，对结构材料来讲，阻碍其力学性能的结构尺度在微米以上，而关于电、光、磁等功能材料来讲可能要小到纳米，甚至是电子结构。因此，计算材料学的研究对象的特征空间尺度从埃到米。时刻是计算材料学的另一个重要的参量。关于不同的研究对象或计算方法，材料计算的时刻尺度可从10-15秒（如分子动力学方法等）到年（如关 下面要紧介绍蒙特卡罗方法： 蒙特卡罗方法： 一、方法的简介 蒙特·卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的进展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类特不重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决专门多计算问题的方法。与它对应的是确定性算法这种方法作为一种独立的方法被提出来，并首先在核武器的试验与研制中得到了应用。蒙特卡罗方法是一种计算方法，但与一般数值计算方法有专门大区不。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。蒙特·卡罗方法在金融工程学，宏观经济学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。 蒙特卡洛方法 1、蒙特卡洛方法的由来 蒙特卡罗分析法（Monte Carlo method），又称为统计模拟法，是一种采用随机抽样（Random Sampling）统计来估算结果的计算方法。由于计算结果的精确度很大程度上取决于抽取样本的数量，一般需要大量的样本数据，因此在没有计算机的时代并没有受到重视。 第二次世界大战时期，美国曼哈顿原子弹计划的主要科学家之一，匈牙利美藉数学家约翰·冯·诺伊曼（现代电子计算机创始人之一）在研究物质裂变时中子扩散的实验中采用了随机抽样统计的手法，因为当时随机数的想法来自掷色子及轮盘等赌博用具，因此他采用摩洛哥著名赌城蒙特卡罗来命名这种计算方法，为这种算法增加了一层神秘色彩。 蒙特卡罗方法提出的初衷是用于物理数值模拟问题, 后来随着计算机的快速发展, 这一方法很快在函数值极小化、计算几何、组合计数等方面得到应用, 于是它作为一种独立的方法被提出来, 并发展成为一门新兴的计算科学, 属于计算数学的一个分支。如今MC方法已是求解科学、工程和科学技术领域大量应用问题的常用数值方法。 2、蒙特卡洛方法的核心—随机数 蒙特卡洛方法的基本理论就是通过对大量的随机数样本进行统计分析，从而得到我们所需要的变量。因此蒙特卡洛方法的核心就是随机数，只有样本中的随机数具有随机性，所得到的变量值才具有可信性和科学性。 在连续型随机变量的分布中, 最基本的分布是[0, 1]区间上的均匀分布, 也称单位均匀分布。由该分布抽取的简单子样ξ1，ξ2ξ3……称为随机数序列, 其中每一个体称为随机数, 有时称为标准随机数或真随机数, 独立性和均匀性是其必备的两个特点。真随机数是数学上的抽象, 真随机数序列是不可预计的, 因而也不可能重复产生两个相同的真随机数序列。真随机数只能用某些随机物理过程来产生, 如放射性衰变、电子设备的热噪音、宇宙射线的触发时间等。 实际使用的随机数通常都是采用某些数学公式产生的，称为伪随机数。真随机数只是一种数学的理想化概念，实际中我们所接触到的和使用的都是伪随机数。要把伪随机数当成真随机数来使用, 必须要通过随机数的一系列的统计检验。 无论伪随机数用什么方法产生,它的局限性都在于这些随机数总是一个有限长的循环集合, 而且序列偏差的上确界达到最大值。所以若能产生低偏差的确定性序列是很有用的,产生的序列应该具有这样的性质, 即任意长的子序列都能均匀地填充函数空间。 人们已经产生了若干种满足这个要求的序列,如Halton序列、Faure序列、Sobol序列和Niederreiter序列等。称这些序列为拟随机数序列。伪随机序列是为了模拟随机性, 而拟随机序列更致力于均匀性。 3、蒙特卡洛方法的原理 当问题可以抽象为某个确定的数学问题时，应当首先建立一个恰当的概率模型，即确定某个随机事件A或随机变量X，使得待求的解等 统计计算课程设计 题目基于蒙特卡洛方法求数值积分 中文摘要 蒙特卡洛方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。 1 利用随机投点法，平均值法，重要性采样法，分层抽样法，控制变量法，对偶变量法，运 用R软件求 1 d x e x θ- =?，42d x e x θ- =?和12 d 1 x e x x θ - = + ?数值积分。计算以上各种估计的方差， 给出精度与样本量的关系，比较各种方法的效率， 关键字蒙特卡洛随机投点法平均值法 R软件 2 1 绪论 蒙特卡洛的基本思想是，当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。 蒙特卡洛方法解题过程的三个主要步骤： （1）构造或描述概率过程 对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 （2）实现从已知概率分布抽样 构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡洛方法模拟实验的基本手段，这也是蒙特卡洛方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡洛模拟的基本工具。 （3）建立各种估计量 一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。 1 期权定价中的蒙特卡洛模拟方法 期权作为最基础的金融衍生产品之一，为其定价一直是金融工程的重要研究领域，主要使用的定价方法有偏微分方程法、鞅方法和数值方法。而数值方法又包括了二叉树方法、有限差分法和蒙特卡洛模拟方法。 蒙特卡洛方法的理论基础是概率论与数理统计，其实质是通过模拟标的资产价格路径预测期权的平均回报并得到期权价格估计值。蒙特卡洛方法的最大优势是误差收敛率不依赖于问题的维数，从而非常适宜为高维期权定价。 §1. 预备知识 ◆两个重要的定理：柯尔莫哥洛夫(Kolmogorov)强大数定律和莱维一林德贝格(Levy-Lindeberg)中心极限定理。 大数定律是概率论中用以说明大量随机现象平均结果稳定性的一系列极限定律。在蒙特卡洛方法中用到的是随机变量序列同分布的Kolmogorov 强大数定律： 设12,,ξξL 为独立同分布的随机变量序列，若 [],1,2,k E k ξμ=<∞=L 则有1 1(lim )1n k n k p n ξμ→∞===∑ 显然，若12,,,n ξξξL 是由同一总体中得到的抽样，那么由 此大数定律可知样本均值1 1n k k n ξ=∑当 n 很大时以概率1收敛于 总体均值μ。 中心极限定理是研究随机变量之和的极限分布在何种情形下是正态的，并由此应用正态分布的良好性质解决实际问题。 设12,,ξξL 为独立同分布的随机变量序列，若 2 [],[],1,2,k k E D k ξμξσ=<∞=<∞=L (0,1)n k d n N ξ μ -??→∑ 其等价形式为2 1 1lim ()exp(),2n x k k n t n P x dt x ξμσ =→∞ -∞ -≤= --∞<<∞∑?。 ◆Black-Scholes 期权定价模型 模型的假设条件： 1、标的证券的价格遵循几何布朗运动 dS dt dW S μσ=+ 其中，标的资产的价格S 是时间t 的函数，μ为标的资产 的瞬时期望收益率，σ为标的资产的波动率，dW 是维纳过程。 2、证券允许卖空、证券交易连续和证券高度可分。 3、不考虑交易费用或税收等交易成本。 4、在衍生证券的存续期内不支付红利。 5、市场上不存在无风险的套利机会。 6、无风险利率r 为一个固定的常数。 下面，通过构造标的资产与期权的资产组合并根据无套利定价原理建立期权定价模型。首先，为了得到期权的微分形式，先介绍随机微积分中的最重要的伊藤公式。 图1-1 蒙特卡罗方法学习总结 核工程与核技术2014级3班张振华20144530317 一、蒙特卡罗方法概述 1.1蒙特卡罗方法的基本思想 1.1.1基本思想 蒙特卡罗方的基本思想就是，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。 1.1.2计算机模拟打靶游戏 为了能更为深刻地理解蒙特卡罗方法的基本思想，我们学习了蒲丰氏问题和打靶游戏两大经典例子。下面主要对打靶游戏进行剖析、计算机模拟（MATLAB 程序）。 设某射击运动员的弹着点分布如表1-1 所示， 首先用一维数轴刻画出已知该运动员的弹 着点的分布如图1-1所示。研究打靶游戏，我 们不用考察子弹的运动轨迹，只需研究每次“扣动扳机”后的子弹弹着点。每一环数对应唯一确定的概率，且注意到概率分布函数有单调不减和归一化的性质。首先我们产生一个在(0,1)上均匀分布的随机数（模拟扣动扳机），然后将该随机数代表的点投到P 轴上（模拟子弹射向靶上的一个确定点），得到对应的环数（即子弹的弹着点），模拟打靶完成。反复进行N 次试验，统计出试验结果的样本均值。样本均值应当等于数学期望值，但允许存在一定的偏差，即理论计算值应该约等于模拟试验结果。 clear all;clc; N=100000;s=0; for n=1:N %step 4.重复N 次打靶游戏试验 x=rand(); %step 1.产生在(0,1)上均匀分布的随机数if(x<=0.1) %step 2.若随机数落在(0.0,0.1)上，则代表弹着点在7环g=7; s=s+g; %step 3.统计总环数elseif(x<=0.2) %step 2.若随机数落在(0.1,0.2)上，则代表弹着点在8环g=8;s=s+g; elseif(x<=0.5) %step 2.若随机数落在(0.2,0.5)上，则代表弹着点在9环g=9;s=s+g; else %step 2.若随机数落在(0.5,1.0)上，则代表弹着点在10环 g=10;s=s+g; end end gn_th=7*0.1+8*0.1+9*0.3+10*0.5; %step 5.计算、输出理论值fprintf('理论值：%f\n',gn_th); gn=s/N; %step 6.计算、输出试验结果 fprintf('试验结果：%f\n',gn);1.2蒙特卡罗方法的收敛性与误差 1.2.1收敛性 由大数定律可知，应用蒙特卡罗方法求近似解，当随机变量Z 的简单子样数N 趋向于无穷大（N 充分大）时，其均值依概率收敛于它的数学期望。 1.2.2误差 由中心极限定理可知，近似值与真值的误差为N Z E Z N αλ<-)(?。式中的αλ的值可以根据给出的置信水平，查阅标准正态分布表来确定。 1.2.3收敛性与误差的关系 在一般情况下，求具有有限r 阶原点矩()∞ Monte Carlo 法 §8.1 概述 Monte Carlo 法不同于前面几章所介绍的确定性数值方法，它是用来解决数学和物理问题的非确定性的（概率统计的或随机的）数值方法。Monte Carlo 方法（MCM ），也称为统计试验方法，是理论物理学两大主要学科的合并：即随机过程的概率统计理论（用于处理布朗运动或随机游动实验）和位势理论，主要是研究均匀介质的稳定状态[1]。它是用一系列随机数来近似解决问题的一种方法，是通过寻找一个概率统计的相似体并用实验取样过程来获得该相似体的近似解的处理数学问题的一种手段。运用该近似方法所获得的问题的解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 计算方法需要有可得的、服从特定概率分布的、随机选取的数值序列。 §8.2 随机数和随机变量的产生 [5]－[10]全面的论述了产生随机数的各类方法。其中较为普遍应用的产生随机数的方法是选取一个函数)(x g ，使其将整数变换为随机数。以某种方法选取 0x ，并按照)(1k k x g x =+产生下一个随机数。最一般的方程)(x g 具有如下形式： m c ax x g mod )()(+= （8.1） 其中 =0x 初始值或种子（00>x ） =a 乘法器（0≥a ） =c 增值（0≥c ） =m 模数 蒙特卡洛模拟法 蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程，反复生成时间序列，计算参数估计量和统计量，进而研究其分布特征的方法。具体的，当系统中各个单元的可靠性特征量已知，但系统的可靠性过于复杂，难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时，可用随机模拟法近似计算出系统可靠性的预计值；随着模拟次数的增多，其预计精度也逐渐增高。由于涉及到时间序列的反复生成，蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的，因此只是在近些年才得到广泛推广。 这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。 蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时，可以用计算机模拟的方法产生抽样结果，根据抽样计算统计量或者参数的值；随着模拟次数的增多，可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。 蒙特卡洛模拟法的应用领域 蒙特卡洛模拟法的应用领域主要有： 1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统，得到某些参数或重要指标。 2.蒙特卡洛积分:利用随机数列计算积分，维数越高，积分效率越高。 3.MCMC:这是直接应用蒙特卡洛模拟方法的推广，该方法中随机数的产生是采用的马尔科夫链形式。 （也叫随机模拟法）当系统中各个单元的可靠性特征量已知，但系统的可靠性过于复杂，难以建立可靠性预计的精确数学模型或模型太复杂而不便应用则可用随机模拟法近似计算出系统可靠性的预计值。随着模拟次数的增多，其预计精度也逐渐增高。由于需要大量反复的计算，一般均用计算机来完成。 应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。解题步骤如下： 1.根据提出的问题构造一个简单、适用的概率模型或随机模型，使问题的解对应于该模型中随机变量的某些特征（如概率、均值和方差等），所构造的模型在主要特征参量方面要与实际问题或系统相一致 2 .根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，方可进行随机模拟试验。 引言 最近在和同学讨论研究Six Sigma（六西格玛）软件开发方法及CMMI相关问题时，遇到了需要使用Monte-Carlo算法模拟分布未知的多元一次概率密度分布问题。于是花了几天时间，通过查询相关文献资料，深入研究了一下Monte-Carl o算法，并以实际应用为背景进行了一些实验。 在研究和实验过程中，发现Monte-Carlo算法是一个非常有用的算法，在许多实际问题中，都有用武之地。目前，这个算法已经在金融学、经济学、工程学、物理学、计算科学及计算机科学等多个领域广泛应用。而且这个算法本身并不复杂，只要掌握概率论及数理统计的基本知识，就可以学会并加以应用。由于这种算法与传统的确定性算法在解决问题的思路方面截然不同，作为计算机科学与技术相关人员以及程序员，掌握此算法，可以开阔思维，为解决问题增加一条新的思路。 基于以上原因，我有了写这篇文章的打算，一是回顾总结这几天的研究和实验，加深印象，二是和朋友们分享此算法以及我的一些经验。 这篇文章将首先从直观的角度，介绍Monte-Carlo算法，然后介绍算法基本原理及数理基础，最后将会和大家分享几个基于Monte-Carlo方法的有意思的实验。所以程序将使用C#实现。 阅读本文需要有一些概率论、数理统计、微积分和计算复杂性的基本知识，不过不用太担心，我将尽量避免过多的数学描述，并在适当的地方对于用到的数学知识进行简要的说明。 Monte-Carlo算法引导 首先，我们来看一个有意思的问题：在一个1平方米的正方形木板上，随意画一个圈，求这个圈的面积。 我们知道，如果圆圈是标准的，我们可以通过测量半径r，然后用S = pi * r^2 来求出面积。可是，我们画的圈一般是不标准的，有时还特别不规则，如下图是我画的巨难看的圆圈。 图1、不规则圆圈 蒙特卡罗方法 为了验证蒙特卡罗方法，我们考虑一个简化的模型：通过一个冷涡轮叶片的热传递。下图是叶片的横截面： 内部冷却通道沿着虚线一截，可以给出： 金属 注意：MM T 金属叶片热的一边的温度 MC T 金属叶片凉的一边的温度 对这个问题，一维热传导模型可以写为： ()()gas gas TBC TBC MH TBC TBC q h T T k q T T L =?=? ()()M MH MC M MC cool cool k q T T L q h T T =?=? 在一个确定的问题里，有四个未知量：，TBC T MM T ，MC T 和q ，我们可以利用阻力求解。同样我们也可以写出下列的线性方程组： 001010440010 01gas gas gas TBC TBC TBC TBC TBC MH M M MC M M cool cool cool h h T T k k L L T k k T L L q h T h ?????????????????????????=?×??????????????????????????? 的线性方程组 输入量是：gas h ，，TBC k M k ， cool h gas T ，，TBC L M L ， cool T 对于一个确定的模拟，我们通常使用标称设计的参数初始值，假设是如下数据： 23000gas W h m = 21000cool W h m = 1300gas T =℃ 200cool T =℃ 1TBC k W m =K 21.5M k W mK = 0.0005TBC L m = 0.003M L m = 模拟的结果如下： 835MH T =℃ 1114TBC T =℃ 758MC T =℃ 525.5810W q m =× 蒙特卡洛模型方法 蒙特卡罗方法（Monte Carlo method） 蒙特卡罗方法概述 蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。 蒙特卡罗方法的提出 蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S．M．乌拉姆和J．冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已经存在。1777年，法国Buffon提出用投针实验的方 样调查来确定可能的优胜者。其基本思想是一样的。 科技计算中的问题比这要复杂得多。比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”（Curse of Dimensionality），传统的数值方法难以对付（即使使用速度最快的计算机）。Monte Carlo 方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。 另一类形式与Monte Carlo方法相似，但理论基础不同的方法—“拟蒙特卡罗方法”（Quasi －Monte Carlo方法）—近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华—王”方法即是其中的一例。这种方法的基本思想是“用确 蒙特卡洛方法及其应用 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的概率分布和数字特征。 蒙特卡洛法的特点是预测结果给出了预测值的最大值,最小值和最可能值，给出了预测 值的区间范围及分布规律。 １.２风险评估概述。 风险表现为损损益的不确定性，说明风险产生的结果可能带来损失、获利或是无损失也无获利，属于广义风险。正是因为未来的不确定性使得每一个项目都存在风险。对于一个公司而言，各种投资项目通常会具有不同程度的风险，这些风险对于一个公司的影响不可小视，小到一个项目投资资本的按时回收，大到公司的总风险、公司正常运营。因此，对于风险的 测量以及控制是非常重要的一个环节。 风险评估就是量化测评某一事件或事物带来的影响的可能程度。根据“经济人”假设，收益最大化是投资者的主要追求目标，面对不可避免的风险时，降低风险，防止或减少损失， 以实现预期最佳是投资的目标。 当评价风险大小时，常有两种评价方式：定性分析与定量分析法。定性分析一般是根据风险度或风险大小等指标对风险因素进行优先级排序，为进一步分析或处理风险提供参考。这种方法适用于对比不同项目的风险程度，但这种方法最大的缺陷是在于，在多个项目中风险最小者也有可能亏损。而定量分析法则是将一些风险指标量化得到一系列的量化指标。通过这些简单易懂的指标，才能使公司的经营者、投资者对于项目分风险有正确的评估与判断， 运用蒙特卡罗模拟进行风险分析 蒙特卡罗模拟由著名的摩纳哥赌城而得名，他是一种非常强有力的方法学。对专业人员来说，这种模拟为方便的解决困难而复杂的实际问题开启了一扇大门。估计蒙特卡罗模拟最著名的早期使用是诺贝尔奖物理学家Enrico Fermi（有时也说是原子弹之父）在1930年的应用，那时他用一种随机方法来计算刚发现的中子的性质。蒙特卡罗模拟是曼哈顿计划所用到的模拟的核心部分，在20世纪50年代蒙特卡罗模拟就用在Los Alamos国家实验室发展氢弹的早期工作中，并流行于物理学和运筹学研究领域。兰德公司和美国空军是这个时期主要的两个负责资助和传播蒙特卡罗方法的组织，今天蒙特卡罗模拟也被广泛应用于不同的领域，包括工程，物理学，研发，商业和金融。 简而言之，蒙特卡罗模拟创造了一种假设的未来，它是通过产生数以千计甚至成千上万的样本结果并分析他们的共性实现的。在实践中，蒙特卡罗模拟法用于风险分析，风险鉴定，敏感度分析和预测。模拟的一个替代方法是极其复杂的随机闭合数学模型。对一个公司的分析，使用研究生层次的高等数学和统计学显然不合逻辑和实际。一个出色的分析家会使用所有他或她可得的工具以最简单和最实际的方式去得到相同的结果。任何情况下，建模正确时，蒙特卡罗模拟可以提供与更完美的数学方法相似的答案。此外，有许多实际生活应用中不存在闭合模型并且唯一的途径就是应用模拟法。那么，到底什么是蒙特卡罗模拟以及它是怎么工作的？ 什么是蒙特卡罗模拟？ 今天，高速计算机使许多过去看来棘手的复杂计算成为可能。对科学家，工程师，统计学家，管理者，商业分析家和其他人来说，计算机使创建一个模拟现实的模型成为可能，这有助于做出预测，其中一种方法应用于模拟真实系统，它通过调查数以百计甚至数以千计的可能情况来解释随机性和未来不确定性。结果通过编译后用于决策。这就是蒙特卡罗模拟的全部内容。 形式最简单的蒙特卡罗模拟是一个随机数字生成器，它对预测，估计和风险分析都很有用。一个模拟计算模型的许多情况，这通过反复地从预先定义的特定变量概率分布中采集数据并将之应用于模型来实现。因为所有的情况都产生相应的结果，每种情况都可以蕴含一种预测。预测的是你定义为重要模型结果的事项（通常含有公式或函数）。 将蒙特卡罗模拟法想象为从一个大篮子里可放回的反复拿出高尔夫球。拦在的大小和形蒙特卡洛法的基本原理
浅析蒙特卡洛方法原理及应用
蒙特卡洛模拟原理及步骤
计算材料学之蒙特卡洛方法论述
蒙特卡洛方法
基于蒙特卡洛方法求数值积分与R
(定价策略)期权定价中的蒙特卡洛模拟方法
蒙特卡罗方法学习总结
(完整版)蒙特卡洛算法详讲
蒙特卡洛模拟法
蒙特卡洛算法
蒙特卡洛模拟——【数学建模 蒙特卡罗算法】
蒙特卡洛模型方法
蒙特卡洛方法及其在风险评估中的应用(1)
运用蒙特卡罗模拟进行风险分析
相关主题
文本预览