一、蒙特卡洛随机模拟

系列一 蒙特卡洛随机模拟

实验目的：学会用计算机随机模拟方法来解决随机性问题

蒙特卡洛模拟法简介

蒙特卡洛(Monte Carlo)方法是一种应用随机数来进行计算机摸你的方法。此方法对研究对象进行随机抽样，通过对样本值的观察统计，求得所研究系统的某些参数。作为随机模拟方法，起源可追溯到18世纪下半叶蒲峰实验。 蒙特卡洛模拟法的应用领域

蒙特卡洛模拟法的应用领域主要有：

1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统，得到某些参数或重要指标。

2.蒙特卡洛积分:利用随机数列计算积分，维数越高，积分效率越高。 蒙特卡洛模拟法求解步骤

应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。解题步骤如下：

1.根据提出的问题构造一个简单、适用的概率模型或随机模型，使问题的解对应于该模型中随机变量的某些特征（如概率、均值和方差等），所构造的模型在主要特征参量方面要与实际问题或系统相一致

2 .根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数，然后生成服从某一分布的随机数，方可进行随机模拟试验。

3. 根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）。

4.按照所建立的模型进行仿真试验、计算，求出问题的随机解。

5. 统计分析模拟试验结果，给出问题的概率解以及解的精度估计。

在可靠性分析和设计中，用蒙特卡洛模拟法可以确定复杂随机变量的概率分布和数字特征，可以通过随机模拟估算系统和零件的可靠度，也可以模拟随机过程、寻求系统最优参数等。

一. 预备知识:

随机数的产生

提示：均匀分布(0, 1)U 的随机数可由C 语言或Matlab 自动产生，在此基础上可产生其他分布的随机数. 1．逆变换法：

设随机变量U 服从(0，1)上的均匀分布，则)(1U F X -=的分布函数为)(x F . 步骤：(1) 产生)1,0(U 的随机数U ；(2) 计算)(1

U F X -=， 则X 服从)(x F 分布. 问题：练习用此方法产生常见分布随机数.例如“指数分布，均匀分布),(b a U ”.还有其它哪种常见分布的随机数可用此方法方便产生？

2．产生离散分布随机数

已知离散随机变量X 的概率分布:)2,1(,)( ===K P x X P k k ，

产生随机变量X 的随机数可采用如下算法：

a) 将区间[0.1]依次分为长度为 ,,21p p 的小区间 ,,21I I ；

b) 产生[0，1]均匀分布随机数R ，若k I R ∈则令k x X =，重复(b)，即得离散随机变量X 的随机数序列.

问题：(1) 下表给出了离散分布X 的概率分布表，试产生100个随机数.

X 的概率分布表

(2) 用此方法给出100个二项分布(20, 0.1)B 的随机数及10个泊松分布P(1)的随机数. 3. 正态分布的抽样

提示：设21,U U 是独立同分布的)1,0(U 变量，令

)

2sin()

ln 2()2cos()ln 2(22

/11222/111U U X U U X ππ-=-=

则1X 与2X 独立 ，均服从标准正态分布. 步骤：(1) 由)1,0(U 独立抽取1122,U u U u ==

(2) 用(*)式计算21,x x .

用此方法可同时产生两个标准正态分布的随机数. 问题： 有关随机数产生方法很多，查阅相关材料进行系统总结. 也可在matlab 中自行产生。

习题：每组做四个题目，其中第四题任选一个。 1. 蒙特卡罗计算π值

思路：图1表示一个内接于正方形的圆的半径R.圆的面积是2

R π，正方形的面积是2(2)R .

圆和正方形的面积的比值就是

2

2

=

(2)4

R R ππ

，将上述比值乘以4，就能获得π值.

图1

2. 蒙特卡罗方法求一元函数的根

根据数学分析相关知识，求一元函数()0f x =的根有很多方法，如一般迭代法，牛顿切线法等等，这些方法讨论根的收敛性，与初始迭代值0x 密切相关.通常要给出一个合适的初始迭代值0x 是比较困难的，利用蒙特卡罗方法可以摆脱根的收敛性对初始值0x 的依赖性.具体方法如下：

要解方程()0f x =，[,]x a b ∈，其中函数()f x 为连续函数，ξ为指定精度. 令down X a =，up X b =，1k =，min 0(1)F F =，步骤如下：

（1）令1k k =+，在[,]down up X X 内产生n 个随机数(1,2,...)i x i n =，计算并比较出这n 个随机数的函数绝对值的最小值1()min{()}nk i i n

f x f x ≤≤=，nk 为i 的某个取值，令

min min ()min{(1),()}x x nk F k F k f x =-.

（2）若min ()x F k ξ<成立，则终止计算，令xroot nk f x =，根就是xroot f ；

若min ()x F k ξ>且min min ()(1)x x F k F k =-，则令1k k =-，转至（1）； 若min ()x F k ξ>且min min ()(1)x x F k F k <-，则令xroot nk f x =，转至（3）. （3）令0/d d k =，down xroot X f d =-，up xroot X f d =+，转至（1）. 说明：a) 0F 是人为给定的一个很大的正数，0d b a <<-且00d >.

b) 1k k =+表示重新赋值给k ，使k 的值增加1，对1k k =-同理. C ）区间[,]down up X X 一定在定义域[,]a b 之内. 此迭代步骤能使函数值序列m i n m i n

m i n

(1)

(2)

...()...

x x x

F F F k >>>>，最终使

min ()x F k ξ<成立，得出函数()0f x =达到精度要求的根xroot f .

例 求3

()tan 8000x f x e x -=-+=在(0,

)2

π

上的实根， (假定-5|()|<10f x ，则认为x 为根)

分析：对上面的方程，如取初值0.5，0.7，0.9，1.1，1.3，1.5，用牛顿迭代公式在区间内找不到根，若用蒙特卡罗方法，则不需要给定初值. 3. 蒙特卡罗方法解线性规划

用蒙特卡罗方法可解约束规划问题：

min (), ()0,1,2,...,.. {,1,2,...,n

X E

i j j j f X g X i m s t

a x

b j n

∈≥=≤≤=

基本思想：在估计的区域{(1x ，2x ，…，n x |j x ∈[,]j j a b ，1,2,...j n =)}内随机取若干个试验点， 然后从试验点中找出可行点，再从可行点中选择出使函数值最小的点.

符号假设和说明：设试验点的第j 个分量j x 服从[,]j j a b 内的均匀分布；P ：试验点总数；

{}MAX P ：最大试验点总数；K ：可行点总数；{}MAX K ：最大可行点数；*X ：迭代产

生的最优点；R ：在[0， 1 ]上的均匀随机数；Q ：迭代产生的最小值*()f X ，其初始值为计算机所能表示的最大数.

求解过程：先产生一个随机数作为初始试验点，以后将上一个试验点的第j 个分量随机产生，其它分量不变而产生一新的试验点.这样，每产生一个新试验点只需一个新的随机数分量.当

{}K MAX K >或{}P MAX P >时停止迭代.

例 求解 22

121212min 283z x x x x x x =--+++，

..s t 12310x x +=，120,0x x ≥≥

蒙特卡罗解线性规划不受问题维数的约束，具有较强的应用性，对于较大维数的问题，只需增加迭代次数就可以获得最优解. 4（1）一重定积分的蒙特卡罗算法

问题描述：假设函数()f x 在[,]a b 内有界连续，且()0f x ≥，求解定积分()b

a

I f x dx =

?

.

为计算出其值，可构造概率模型如下：取一个边长分别为b a -和c 的矩形D ，使曲边梯形在矩形域之内，如图2，并在矩形内随机投点，假设随机点均匀地落在整个矩形之内，则落在图中灰色区域内的随机点数k 与投点总数N 之比k/N 就近似地等于曲线下方面积（即阴影面积）与矩形面积之比，从而得出近似积分()k

I b a c N

≈

-.

图2

例 求2

1

1

x e

--?

由于2

x e -是非初等函数，我们很难求出其原函数，所以用牛顿-莱布尼茨公式无法求解，但可以运用蒙特卡罗方法求出其近似值.

将上述方法推广到一般情况：假设函数()f x 在[a ，b]内有界连续，对于定积分

()b

a

I f x dx =?，为计算出其值，可构造如下概率模型：

取一个边长分别为b a -和c d -的矩形D ，使曲线[,]a b 段的值在矩形域之内，如图3，并在矩形内随机投点，假设随机点均匀地落在整个矩形之内，则落在图中x 轴上下灰色区域内的随机点数m 与n 的差与投点总数p 之比(m-n)/P 就近似地等于曲线上下方面积之差（即阴影面积之差）与矩形面积之比，从而得出近似积分()()m n

I b a c d P

-≈

--.

图3

（2）. 二重积分的蒙特卡罗算法 问题描述：实际计算中常常要遇到如

(,)D

f x y dxdy ??的二重积分，发现被积函数的原函数往

往很难求出，或者原函数根本就不是初等函数，对于这样的重积分，蒙特卡罗方法也有成熟

的计算方法. 方法1： 步骤：

1，取一个包含D 的矩形区域Ω：,a x b c y d ≤≤≤≤，面积()()A b a d c =--；

2，(,), 1,2,,i i x y i n = ，为Ω上的均匀分布随机数列，不妨设(,),1,2,i i x y i n = （）为落在D 中的n 个随机数，则n 充分大时，有1

(,)(,)k

i i i D

A f x y dxdy f x y n =≈∑??

.

方法2：

对二重积分(,)A

I f x y dxdy =

??

，假设(,)f x y 为区域A 上的有界函数，且(,)0f x y ≥，

几何意义对应的是以(,)f x y 为曲面顶， A 为底的曲顶柱体C 的体积.因此，用均匀随机数计算二重积分的蒙特卡罗方法基本思路为：假设曲顶柱体C 包含在己知体积为D

V

的几何体

D 的内部，在D 内产生N 个均匀随机点，统计出在C 内部的随机点数目C N ，则D

C V I N N

=.

例：计算

(1A

dxdy +

??，其中22{(,)|1}A x y x y =+≤.

分析：该二重积分可以看作以1+曲顶柱体在一个边长为2的立方体内，用数学分析方法可计算出其精确值为π. （3）用蒙特卡洛方法计算体积

,

冰激凌锥含于体积=8的六面体