MonteCarlo(蒙特卡洛算法)算法

MonteCarloSimulation简介概述蒙特卡罗(MonteCarlo)方法，或称计算机随机模拟方法或随机抽样方法或统计试验方法，属于计算数学的一个分支。是一种基于“随机数”的计算方法。起源MonteCarlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定圆周率π。成型这一方法成型于美国在第一次世界大战进研制原子弹的“曼哈顿计划”。该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的MonteCarlo—来命名这种方法，为它蒙上了一层神秘色彩。发展本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。实质MonteCarlo方法也称为统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。与它对应的是确定性算法。把一些复杂的东西用大量的模拟实验来做，最后得到一些结论。基本思想和原理基本思想：当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。原理：抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。。步骤可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量构造或描述概率过程对于本身就具有随机性质的问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。实现从已知概率分布抽样构造了概率模型以后，按照这个概率分布抽取随机变量（或随机向量），这一般可以直接由软件包调用，或抽取均匀分布的随机数构造。这样，就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。建立各种估计量一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。例子考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？MonteCarlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点落于“图形”内，则该“图形”的面积近似为M/N。比喻可用民意测验来作一个不严格的比喻。民意测验的人不是征询每一个登记选民的意见，而是通过对选民进行小规模的抽样调查来确定可能的民意。其基本思想是一样的。应用科技计算中的问题比这要复杂得多。但MonteCarlo方法广泛地应用于许多应用领域，如计算物理学、粒子输运计算、量子热力学计算、量子化学、分子动力学与。特别在金融计算中，各方法有不可取代的优势。金融中的应用金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”(CourseDimensionality)，传统的数值方法难以对付（即使使用速度最快的计算机）。MonteCarlo方法的优势MonteCarlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。MonteCarlo模拟适用于研究复杂体系。研究具有多得数不清的结构、状态的体系，对此我们可以采用蒙特卡洛模拟，以统计的方法寻找出现几率最高的结构、状态，或相应的有关数据。MonteCarlo方法处理的问题MonteCarlo方法处理的问题可以分两类确定性的数学问题多重积分、求逆矩阵、解线性代数方程组、解积分方程、解某些偏微分方程边值问题和计算代数方程组、计算微分算子的特征值等等随机性问题方法在解决实际问题的时候应用MonteCarlo方法主要有两部分工作：1、用此方法模拟某一过程时，需要产生各种概率分布的随机变量。2、用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。用MonteCarlo计算定积分考虑积分假定随机变量具有密度函数则10,0.xIxedx(),xXfxe1().IEX用MonteCarlo计算定积分-抽取密度为e^{-x}的随机数X_1,…X_n构造统计数则111ˆ.niiIXn111111111ˆ()1.nniiiiniEIEXEXnnnEXEXInn用MonteCarlo计算定积分--且即1121111211ˆ()1.nniiiiVarIVarXVarXnnnVarXVarXnn11.IXn用MonteCarlo计算定积分---例如α=1.9取Γ(1.9)=0.96176–模拟结果不好！如果要达到0.001的精确度，要4X530^2=1123600计算！0.91.90.xIxedx1234ln,(0,1):0.0587,0.0961,0.9019,0.3095,ˆ1.497iiiXRRURRRRI用MonteCarlo计算定积分----例子说明分析和设计是重要的。重写积分取两个随机数1.90.101.()xyYIxedxfyyex密度函数为12121.90.11,(0,1),ln,RRUYRRIEY令算用MonteCarlo计算定积分-----取8个随机数大大改善了结果！123456781.90.0078,0.93250.10800.00630.5490,0.85560.97710.2783ˆ0.9187RRRRRRRRI，，，，随机数的产生随机数是我们实现蒙特卡罗模拟的基本工具。随机数的产生就是抽样问题。可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。随机数的取得如果你对随机数有更高的要求，需要自己编辑“随机数生成器”最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）例如在Matlab中，命令“rand()”将产生一个(0,1)中均匀分布的随机数你可以根据需要给随机数一个“种子”，以求不同的数Matlab的随机数函数均匀分布R=unidrnd(N)，-产生1到N间的均匀分布随机数R=unidrnd(N,n,m),产生1到N间的均匀分布随机数矩阵连续均匀分布R=unifrnd(A,B)-产生(A,B)间的均匀分布随机数R=unifrnd(A,B,m,n)产生(A,B)间的均匀分布随机数矩阵Matlab的随机数函数-正态分布随机数R=normrnd(mu,sigma)R=normrnd(mu,sigma,m)R=normrnd(mu,sigma,m,n)特定分布随机数发生器R=random(‘name’,A1,A2,A3,m,n)例a=random(‘Normal’,0,1,3,2)a=.-0.43260.2877-1.6656-1.14650.12531.1909精确性由于MonteCarlo方法的随机性，精确性建立在大量的重复模拟上，最后去平均值。对确定值的计算，要估计出样本的个数与精确度之间的关系。对随机过程的模拟，有置信区域的估算等方差削减技术对偶变量技术(适用正态分布函数）取一组随机数Z_i，可得模拟值C_i,i=1,2,..n估计值为期平均C^再取Z_i的对偶Z’_i=-Z_i,再生成估计值C’^然后去新的平均值C*=（C^+C’^)/2则varC*=1/2varC^+1/2cov(C^,C’^)1/2varC^+该技术使计算更稳定一个例子问题下图是一个中子穿过用于中子屏蔽的铅墙示意图。铅墙的高度远大于左右厚度。设中子是垂直由左端进入铅墙，在铅墙中运行一个单位距离然后与一个铅原子碰撞。碰撞后，任意改变方向，并继续运行一个单位后与另一个铅原子碰撞。这样下去，如果中子在铅墙里消耗掉所有的能量或者从左端逸出就被视为中子被铅墙挡住，如果中子穿过铅墙由右端逸出就视为中子逸出。如果铅墙厚度为5个单位，中子运行7个单位后能量耗尽，求中子逸出的几率。x这个问题并不复杂，但不容易找到一个解析表达式。而用模拟的方法求解却可以有满意的结果。一个例子-这个问题并不复杂，但不容易找到一个解析表达式。而用模拟的方法求解却可以有满意的结果。一个例子--建模下面我们给出这个问题的模拟程序。我们关心的是一次碰撞后，中子在x轴方向行进了多少，所以行进方向是正负θ的结果是一样的，我们就只考虑θ是正的情形。由于中子运行的方向θ是随机的，我们用计算机抽取在0到π间均衡分布的随机数，模拟1000000个中子在铅墙里行进的情形，看看这些中子与铅原子碰撞7次后，有多少超过了铅墙的右端。一个例子---n=1000000;m=0;t=1;fori=1:nx=1;fork=1:7ang=pi*rand;x=x+cos(ang);ifx0l=0;t=0;endendifx5&t==1l=1;elsel=0;endm=m+l;t=1;endm/n一个例子----解模我们运行程序得出逸出铅墙的中子的可能性约为1.5%。应用有了这个数字，我们可以报告安全部门，如果数字不能达到安全要求，我们则要加厚铅墙。MonteCarlo模拟二叉树期权定价二叉树定价模型是从构造好的二叉树中随机选择一条路径样本，从二叉树的末端开始倒推计算出衍生证券的价格，但是采用了MonteCarlo后，是顺着二叉树往后计算的。基本方法：在第一个节点（根节点），随机产生一个0到1间的随机数，如果这个数小于p，就选择当前的上升分支，反之选择下降分支。这样就产生了一个新节点，继续上面的过程，直到二叉树的末端。一条路径产生了，衍生证券的最终价值就可以计算出来了（可以看作是全部可能终值集合中的一个随机样本），这样完成了第一次模拟。更多的样本路径得到更多的样本终值。进行多次模拟，用平均值来估计衍生证券的价格MonteCarlo模拟连续过程的欧式期权定价欧式期权定价的期望公式为如果标的资产服从几何布朗运动则风险中性测度下，标的资产的过程为ˆ()rTTVeEVdSSdtSdW20exp2TSSrTTMonteCarlo模拟连续过程的欧式期权定价-所以CallOption到期日的现金流为抽一个正态分布的随机数ε，算得上式再贴现，就可以得到一个结果重复上面的过程（如计算10000次）将所有的结果取平均2max0,(0)exp2SrTTK附表:Matlab随机数生成器betarnd贝塔分布binornd二项分布chi2rnd卡方分布exprnd指数分布frndf分布gamrnd伽玛分布geornd几何分布hygernd超几何分布lognrnd对数正态分