实验十五：MATLAB的蒙特卡洛仿真

实验十五： MATLAB 的蒙特卡洛仿真

一、实验目的

1. 了解蒙特卡洛仿真的基本概念。

2. 了解蒙特卡洛仿真的某些应用

二．实验内容与步骤

1． 蒙特卡洛（Monte Carlo ）仿真的简介

随机模拟方法，也称为Monte Carlo 方法，是一种基于“随机数”的计算方法。这一方法源于美国在第一次世界大战进行的研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo 来命名这种方法，为它蒙上了一层神秘色彩。冯·诺伊曼是公理化方法和计算机体系的领袖人物，Monte Carlo 方法也是他的功劳。

事实上，Monte Carlo 方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。18世纪下半叶的法国学者Buffon 提出用投点试验的方法来确定圆周率π的值。这个著名的Buffon 试验是Monte Carlo 方法的最早的尝试！

历史上曾有几位学者相继做过这样的试验。不过他们的试验是费时费力的，同时精度不够高，实施起来也很困难。然而，随着计算机技术的飞速发展，人们不需要具体实施这些试验，而只要在计算机上进行大量的、快速的模拟试验就可以了。Monte Carlo 方法是现代计算技术的最为杰出的成果之一，它在工程领域的作用是不可比拟的。

蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程，反复生成时间序列，计算参数估计量和统计量，进而研究其分布特征的方法。具体的，当系统中各个单元的可靠性特征量已知，但系统的可靠性过于复杂，难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时，可用随机模拟法近似计算出系统可靠性的预计值；随着模拟次数的增多，其预计精度也逐渐增高。由于涉及到时间序列的反复生成，蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的，因此只是在近些年才得到广泛推广。

2． MC 的原理

针对实际问题建立一个简单且便于实现的概率统计模型，使问题的解对应于该模型中随机变量的概率分布或其某些数字特征，比如，均值和方差等。所构造的模型在主要特征参量方面要与实际问题或系统相一致的。

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

收敛性: 由大数定律, Monte-Carlo 模拟的收敛是以概率而言的．

误差: 用频率估计概率时误差的估计，可由中心极限定理，给定置信水平 的条件下，有： 模拟次数:由误差公式得 3． 定积分的MC 计算原理

事实上，不少的统计问题，如计算概率、各阶距等，最后都归结为定积分的近似计算问题。设 a,b ，有限， , (){}M y b x a y x ≤≤≤≤=Ω0,:,并设()Y X ,是在Ω

N U σεα2

/1||-≤))((X g Var =σ()M x f ≤≤0

上均匀分布的二维随机变量，其联合密度函数为

()()M y b x a I a b M ≤≤≤≤-0,1 。则易见 ()dx x f b a

⎰=θ 是 Ω 中 ()x f y = 曲线下方的面积。假设我们向Ω 中进行随机投点，若点落在 ()x f y =下方，(即 ()x f y <称为中的，否则称为不中，则点中的概率为 ()a b M p -=θ

。若我们进行了 n 次投点，其中0n 次中的，则用频率来估计概率p 。即n

n p 0= 。 4． 蒙特卡洛仿真应用举例

例1 计算定积分

事实上，其精确解为

用随机投点法求解如下：

sjtdf(0,4,4,1000000) result = 7.2336

function result=sjtdf(a,b,m,mm)

%a 是积分的下限

%b 是积分的上限

%m 是函数的上界

%mm 是随机实验次数

frq=0;

xrandnum = unifrnd(a,b,1,mm);

yrandnum = unifrnd(0,m,1,mm);

for ii=1:mm

if (cos(xrandnum(1,ii))+2>=yrandnum(1,ii))

frq=frq+1;

end

end

result=frq*m*(b-a)/mm

例2 π的计算

（单位圆的面积等于 π）

>> sjtdf_pi1(1000) piguji = 3.0520

>> sjtdf_pi1(10000) piguji = 3.1204

>> sjtdf_pi1(100000) piguji = 3.1296

function piguji=sjtdf_pi1(mm)

%mm 是随机实验次数

frq=0;xrandnum = unifrnd(0,1,1,mm);

yrandnum = unifrnd(0,1,1,mm);

for ii=1:mm

if xrandnum(1,ii)^2+yrandnum(1,ii)^2<=1

frq=frq+1;

end

end

piguji=4*frq/mm

思考题

运用定积分的MC计算原理，用随机投点法计算二重定积分