蒙特卡罗方法介绍及其建模应用---视频课程课件-2016-03-31

Monte-Carlo方法介绍及其建模应用朱连华南京信息工程大学2020年4月22日2020/4/2216:27南京信息工程大学主要内容蒙特卡洛方法应用实例22009-B眼科病床安排应用3蒙特卡洛方法介绍12020/4/2216:27南京信息工程大学蒙特卡洛方法介绍蒙特卡洛的优缺点及其适用范围2随机数的生成3蒙特卡洛方法概述1MonteCarlo的起源设总计投了M个点，落入阴影部分N个，则不规则图形的面积为MNS蒙特卡罗方法数学建模的解有两类•精确解•近似解自然现象有两类•确定性现象•不确定性现象•随机现象•模糊现象12020/4/2216:27南京信息工程大学随机投点试验求近似解引例---计算面积在一定条件下必然发生现象2020/4/2216:27南京信息工程大学MonteCarlo的起源MonteCarlo方法：–又称随机模拟方法，对研究的系统进行随机观察抽样,通过对样本值的统计分析，求得所研究系统的某些参数–不同于一般数值计算方法，它是以概率统计理论为基础的一种方法，由于能够比较逼真地描述事物的特点及物理实验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛–它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来，“曼哈顿计划”主持人之一、数学家：冯·诺伊曼用驰名世界的赌城—摩纳哥最大的城市MonteCarlo—来命名这种方法Monte-CarloMonacoJohnVonNeumann(1903-1957)2020/4/2216:27南京信息工程大学MonteCarlo方法的应用物理：核物理，热力学与统计物理，粒子输运问题等数学：多重积分、解微分方程、非线性方程组求解等工程领域：真空技术，水力学，激光技术等经济学领域：期权定价、项目管理、投资风险决策等其他领域：化学、医学，生物，生产管理、系统科学、公用事业等方面，随着科学技术的发展，其应用范围将更加广泛。2020/4/2216:27南京信息工程大学MonteCarlo方法的基本思想蒲丰（Buffon）投针实验:为了求得圆周率𝜋值，将长为2𝑙的一根针任意投到地面上，用针与一组相间距离为2𝑎（𝑙𝑎）的平行线相交的频率𝑓替代概率𝑝，再利用准确的关系式：𝑝=2𝑙𝜋𝑎求出𝜋值:𝜋=2𝑙𝑎𝑝≈2𝑙𝑎(𝑁𝑛)其中𝑁为投针次数，𝑛为针与平行线相交次数，这就是古典概率论中著名的蒲丰氏问题！法国数学家ComtedeBuffon例1.（蒲丰投针问题）为了求得圆周率π值，在十九世纪后期，有很多人作了这样的试验：将长为2l的一根针任意投到地面上，用针与一组相间距离为2a（l＜a）的平行线相交的频率代替概率P，再利用准确的关系式：alP2求出π值)(22nNalaPl其中Ｎ为投计次数，n为针与平行线相交次数。这就是古典概率论中著名的蒲丰氏问题。一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929解：设针投到地面上的位置可以用一组参数（x,θ）来描述，x为针中心的坐标，θ为针与平行线的夹角，如图所示。任意投针，就是意味着x与θ都是任意取的，但x的范围限于［0，a］，夹角θ的范围限于［0，π］。在此情况下，针与平行线相交的数学条件是sinlx如何产生任意的（x,θ）？x在［0，a］上任意取值，表示x在［0，a］上是均匀分布的，其分布密度函数为：其他,00,/1)(1axaxf类似地，θ的分布密度函数为：2020/4/2216:27南京信息工程大学MonteCarlo方法的基本思想解：针投到地面上的位置用一组参数𝑥,𝜃来刻画，𝑥为针中心坐标，𝜃为针与平行线的夹角（如图），由于为任意投针，𝑥,𝜃均匀落在0,𝑎,0,𝜋,针与平行线相交的条件为：𝑥≤𝑙𝑠𝑖𝑛𝜃sin00sin12lpPxldxdala每次投针试验实际上变成从两个均匀分布的随机变量中抽样，从而针线相交的概率为:𝜋=2𝑙𝑎𝑝≈2𝑙𝑎(𝑁𝑛)2020/4/2216:27南京信息工程大学计算机模拟程序实现：functionpiguji=buffon(llength1,llength2,mm)%llength1,llength2分别表示1/2线的宽度和针的长度%mm是随机实验次数frq=0;xrandnum=unifrnd(0,llength1,1,mm);theta=unifrnd(0,pi,1,mm);forii=1:mmif(xrandnum(1,ii)=(llength2*sin(theta(1,ii))))frq=frq+1;endendpiguji=(2*llength2/llength1)/(frq/mm)endbuffon(1,.6,1000)piguji=3.14622020/4/2216:27南京信息工程大学①建立统计模型，主要特征参量方面要与实际问题或系统相一致，问题的解对应于模型中随机变量的概率分布或其某些数字特征②根据模型中各个随机变量的分布，在计算机上产生随机数，实现一次模拟过程所需的足够数量的随机数，进而进行随机模拟实验③根据概率模型的特点和随机变量的分布特性，设计和选取合适的抽样方法，并对每个随机变量进行抽样（包括直接抽样、分层抽样、相关抽样、重要抽样等）④按照所建立模型进行仿真试验、计算，求出问题的随机解⑤统计分析模拟试验结果，给出问题的估计以及其精度估计。必要时，还应改进模型以降低估计方差和减少试验费用，提高模拟计算的效率。蒙特卡洛方法计算机模拟基本步骤2020/4/2216:27南京信息工程大学蒙特卡洛模拟的理论基础大数定律---贝努里（Bernoulli）大数定律中心极限定理1limpnnPAn()PAnpnn1~(0,1)()nkkXnNnn~(0,1)/nXNn2020/4/2216:27南京信息工程大学蒙特卡洛模拟的误差分析由中心极限定理可知：这表明，不等式近似地以概率1成立。上式也表明，收敛到的阶为O(n-1/2)。通常，蒙特卡罗方法的误差ε定义为：nuXnnX1)(unuXPNnu2020/4/2216:27南京信息工程大学蒙特卡洛方法优缺点及其适用范围MonteCarlo方法及其程序结构简单–产生随机数，通过大量简单重复抽样和简单计算计算相应的值收敛速度与问题维数无关–MonteCarlo方法的收敛速度为O(n-1/2)，与一般数值方法相比很慢。因此，用MonteCarlo方法不能解决精确度要求很高的问题–MonteCarlo方法误差只与标准差和样本容量n有关，而与样本所在空间无关，即MonteCarlo方法的收敛速度与问题维数无关，而其他数值方法则不然。MonteCarlo方法的适用性强–MonteCarlo方法对多维问题的适用性–在解题时受问题条件限制的影响较小例如：要计算s维空间中的任一区域Ds上的积分2020/4/2216:27南京信息工程大学蒙特卡洛方法优缺点及其适用范围优点1.能够比较逼真地描述具有随机性质的事物的特点及物理实验过程2.受几何条件限制小3.收敛速度与问题的维数无关4.误差容易确定5.程序结构简单，易于实现缺点1.收敛速度慢。2.误差具有概率性。3.进行模拟的前提是各输入变量是相互独立的蒙特卡罗方法适应于随机性问题，对于非随机性问题，必须将它转化为随机问题2020/4/2216:27南京信息工程大学Matlab中的随机数生成函数rand(m,n):(0,1)均匀分布unifrnd(a,b,m,n):(a,b)均匀分布randn(m,n):标准正态normrnd(mu,sigma,m,n):正态分布exprnd(theta,m,n):指数分布poissrnd(lamda,m,n):泊松分布binornd(n,p,m,n):二项分布hist(poissrnd(3,1,1000),50)2020/4/2216:27南京信息工程大学Matlab中的随机数生成函数random('name',A1,A2,A3,M,N)name的取值可以是'norm'or'Normal''unif'or'Uniform''poiss'or'Poisson''beta'or'Beta''exp'or'Exponential''gam'or'Gamma''geo'or'Geometric''unid'or'DiscreteUniform'......2020/4/2216:27南京信息工程大学一般分布随机数产生方法基本方法有如下三种：–逆变换法–复合抽样方法–筛选法定理设随机变量U服从(0,1)上的均匀分布，则X=F-1(U)的分布函数为F(x)因此，要产生来自F(x)的随机数，只要先产生来自U(0,1)的随机数，然后计算F-1(u)即可。其步骤为uFxuU11,0计算，抽取由2020/4/2216:27南京信息工程大学连续分布对于连续型分布，如果分布函数F(x)的反函数F－1(x)能够解析表示，则直接抽样方法是:1(2)()FxFu(1)0,1Uu由抽取2020/4/2216:27南京信息工程大学指数分布为连续型分布，其一般形式如下：其分布函数为：则(1)由U(0,1)抽取u因为1-u也是(0,1)上均匀随机数，可将上式简化为0,1)(xexFx0,)(xexfx)1ln(1uxF所以/lnuxF(2)(),1ln(1)FxFFFxueuxu令得，即例.指数分布2020/4/2216:27南京信息工程大学Randnum=（-2）*log(unifrnd(0,1,1,1000));cdfplot(Randnum)0,21)(21xexfx例.产生指数分布的随机数)ln(2uxF例：病人视网膜疾病术后观察时间一般为[5,15]天，其人数分布如直方图所示，求满足该分布的术后观察天数序列。离散分布天数概率分布函数值区间53/1013/101162/1015/101277/10112/1013816/10128/1014911/10139/10151018/10157/10161115/10172/10171211/10183/1018139/10192/1019145/10197/10110154/1011112020/4/2216:27南京信息工程大学离散分布x=[5:15];X=[];y=[3,2,7,16,11,18,15,11,9,5,4];yy=y(1);fori=2:length(y)yy=[yy,y(i)+yy(i-1)];endyy=yy/yy(end);fori=1:10000r=rand(1);d=size(find(r=yy),2)+1;X(i)=x(d);Endhist(X,x)2020/4/2216:27南京信息工程大学2020/4/2216:27南京信息工程大学蒙特卡洛方法应用实例2020/4/2216:27南京信息工程大学定积分的MC计算事实上，不少的统计问题最后都归结为定积分的近似计算问题!相对于其它方法，用MC方法比一般的数值方法有优点，主要体现在它的误差与维数m无关！下面考虑一个简单的定积分为了说明问题，这里介绍两种求的简单的MC方法。dxxfba2020/4/2216:27南京信息工程大学方法简述：设a，b有限，0f(x)M，={(x,y):axb，0yM},并设(X,Y)是在上均匀分布的二维随机变量，其联合密度函数为MybxaIabM0,1则易见是中y=f(x)曲线下方的面积假设我们向中进