系列一 蒙特卡洛随机模拟实验目的：学会用计算机随机模拟方法来解决随机性问题 蒙特卡洛模拟法简介蒙特卡洛(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) 计算)(1U 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.圆的面积是2R π，正方形的面积是2(2)R .圆和正方形的面积的比值就是22=(2)4R R ππ，将上述比值乘以4，就能获得π值.图12. 蒙特卡罗方法求一元函数的根根据数学分析相关知识，求一元函数()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 nf 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 之内.此迭代步骤能使函数值序列min min min (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 Ei 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 >时停止迭代.例 求解 22121212min 283z x x x x x x =--+++，..s t 12310x x +=，120,0x x ≥≥蒙特卡罗解线性规划不受问题维数的约束，具有较强的应用性，对于较大维数的问题，只需增加迭代次数就可以获得最优解.4（1）一重定积分的蒙特卡罗算法问题描述：假设函数()f x 在[,]a b 内有界连续，且()0f x ≥，求解定积分()ba I f x dx =⎰.为计算出其值，可构造概率模型如下：取一个边长分别为b a -和c 的矩形D ，使曲边梯形在矩形域之内，如图2，并在矩形内随机投点，假设随机点均匀地落在整个矩形之内，则落在图中灰色区域内的随机点数k 与投点总数N 之比k/N 就近似地等于曲线下方面积（即阴影面积）与矩形面积之比，从而得出近似积分()k I b a c N≈-.图2例 求211x e --⎰由于2x e -是非初等函数，我们很难求出其原函数，所以用牛顿-莱布尼茨公式无法求解，但可以运用蒙特卡罗方法求出其近似值.将上述方法推广到一般情况：假设函数()f x 在[a ，b]内有界连续，对于定积分()b aI 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）. 二重积分的蒙特卡罗算法问题描述：实际计算中常常要遇到如(,)Df 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=. 例：计算2222(11)Ax y x y dxdy +---+⎰⎰，其中22{(,)|1}A x y x y =+≤. 分析：该二重积分可以看作以222211x y x y +---+为顶的曲顶主体的体积，此曲顶柱体在一个边长为2的立方体内，用数学分析方法可计算出其精确值为π.（3）用蒙特卡洛方法计算体积,冰激凌锥含于体积=8的六面体。