蒙特卡洛法-算法

程序设计方法之蒙特卡罗(MonteCarlo)法MonteCarlo法蒙特卡罗方法（MonteCarlomethod），也称统计模拟方法一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法基本思想：当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解MonteCarlo法简介MonteCarlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定圆周率π。本世纪40年代计算机的出现、特别是近年来高速计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。使用蒙特卡罗方法估算π值.放置30000个随机点后,π的估算值与真实值相差0.07%.蒙特卡洛方法应用用蒙特卡洛方法模拟某一过程时，需要产生各种概率分布的随机变量。用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。应用举例考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？MonteCarlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷M个点，其中有N个点落于“图形”内，则该“图形”的面积近似为：N/M。科技计算中的问题比这要复杂得多。比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”，传统的数值方法难以对付（即使使用速度最快的计算机）。MonteCarlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。C语言中伪随机数生成函数rand()产生随机数的函数，它可生成0~RAND_MAX内的整数RAND_MAX在cstdlib(stdlib.h)中声明，32,767srand(unsigndeint)设置产生随机数rand函数产生随机数的种子在调用rand函数之前调用srand设置种子，如果不调用rand使用一个默认值作为种子举例：产生10个随机整数#includeiostream//预编译命令#includecstdlib//预编译命令usingnamespacestd;intmain(){intk=0;//定义整型变量kfor(k=0;k10;k++)//循环输出随机数coutrand()““;coutendl;cout“rand能产生的最大随机数为RAND_MAXendl;//输出最大随机数return0;}//主函数结束使用srand/rand产生随机数过程1、给srand()提供一个“种子”，它是一个unsignedint类型，其取值范围是从0到65,535;2、调用rand()，它会根据提供给srand()的“种子”值返回一个随机数(在0到32，767之间)；3、根据需要多次调用rand()，从而不断地得到新的随机数；4、无论什么时候，你都可以给srand()提供一个新的“种子”，从而进一步“随机化rand()的输出结果。问题：如果你每次调用srand()时都提供相同的“种子”值，那么你将会得到相同的“随机”数序列。一种简单而有效的办法来产生一个相当随机的“种子”值——当前的时间值。举例：产生10个随机整数#includeiostream//预编译命令#includecstdlib//预编译命令#includectime//预编译命令usingnamespacestd;intmain(){intk=0;//定义整型变量ksrand((unsignedint)time(NULL));//设置种子for(k=0;k10;k++)//循环输出随机数coutrand()““;coutendl;cout“rand能产生的最大随机数为RAND_MAXendl;//输出最大随机数return0;}//主函数结束C++11中随机数函数库random最新的2011C++标准中有随机数库函数，对产生随机数有更好的支持注意：旧版本的编译器不支持includerandom举例：std::default_random_enginestd::normal_distribution可以产生正态分布的随机小数，可指定均值和方差课后练习产生n个随机小数（题库作业）蒙特卡洛法举例：求π的近似值如右图，正方形的面积A=1；1/4圆的面积B=π/4。我们想象有一个容器在正方形中夹有一个极薄的圆弧隔板。下小雨时搬至屋外，经一定时间后，称1/4圆的容器内的水重C，与作为一个整体的正方形中的水重D。C与D之比应该等于B与A之比，即可得𝐶𝐷=𝐵𝐴𝜋=4×𝐶𝐷求π的近似值(2)让计算机来模拟雨点落下：产生伪随机数x和y，让x的值的范围在0～1之间；让y的值的范围也在0～1之间，模拟雨点落在正方形中，当然会有的雨点落在1/4圆中。数以百万计雨点可以累计得到C和D，从而上述公式算出π的近似值。关键点：落入扇形区的判据𝑥2+𝑦2≤1代码#includeiostream//预编译命令#includecstdlib//预编译命令#includectime//预编译命令#includecmath//预编译命令usingnamespacestd;intmain()//主函数{longk=0,c=0,d=0;//定义长整型变量floatpai=0.0,x=0.0,y=0.0;//定义浮点类型变量srand((unsignedint)time(NULL));//设置种子for(k=1;k=10000000;k++){//循环体开始d=d+1;//累加正方形中落入的一个雨点x=(float)rand()/32767;//雨点在x方向的位置y=(float)rand()/32767;//雨点在y方向的位置if(sqrt(x*x+y*y)=1)c=c+1;//累加扇形中落入的一个雨点}pai=4.0f*(float)c/d;//计算pai的值coutpai=paiendl;//输出pai的值return0;}//主函数结束代码（续）举例：计算阴影部分面积乍看似乎有些难，但有了前面的基础就不难了。思路：考虑落在阴影这块面积上的雨点数g与落在正方形整体上的雨点数d的比就是阴影部分面积的近似值。𝑥2+𝑦2=1(𝑥−1)2+𝑦2=1代码#includeiostream//预编译命令#includecstdlib//预编译命令#includectime//预编译命令#includecmath//预编译命令usingnamespacestd;intmain()//主函数{//主函数开始longk=0,d=0,g=0;//定义长整型变量floats=0.0,x=0.0,y=0.0;//定义浮点类型变量srand((unsignedint)time(NULL));//设置种子for(k=1;k=10000000;k++){//循环体开始d=d+1;//累加正方形中落入的一个雨点x=(float)rand()/32767;//雨点在x方向的位置y=(float)rand()/32767;//雨点在y方向的位置if((sqrt(x*x+y*y)1)&&(sqrt((x-1)*(x-1)+y*y)1))g=g+1;//累加阴影面积中落入的一个雨点}代码（续）s=(float)g/d;//计算s的近似值couts的近似值为sendl;//输出s的近似值//输出s的精确值，比较用couts的精确值为1.0-(3.14159/6.0+1.73205/4.0)endl;return0;}为了比较所得结果是否有效，程序中给出了s的精确值，这个值是用几何的解法计算出的，见下图无阴影线部分可分解为ABA三块AAB30°30°