北京物流信息联盟

解析Monte-Carlo算法(基本原理,理论基础,应用实践)

2021-06-16 12:49:39

限时干货下载:

回复【69】免费获取【精选R语言大礼包

来源:http://www.cnblogs.com/leoo2sk/archive/2009/05/29/1491526.html

作者:张洋


1

引言


最近在和同学讨论研究Six Sigma(六西格玛)软件开发方法及CMMI相关问题时,遇到了需要使用Monte-Carlo算法模拟分布未知的多元一次概率密度分布问题。于是花了几天时间,通过查询相关文献资料,深入研究了一下Monte-Carlo算法,并以实际应用为背景进行了一些实验。
      

在研究和实验过程中,发现Monte-Carlo算法是一个非常有用的算法,在许多实际问题中,都有用武之地。目前,这个算法已经在金融学、经济学、工程学、物理学、计算科学及计算机科学等多个领域广泛应用。而且这个算法本身并不复杂,只要掌握概率论及数理统计的基本知识,就可以学会并加以应用。由于这种算法与传统的确定性算法在解决问题的思路方面截然不同,作为计算机科学与技术相关人员以及程序员,掌握此算法,可以开阔思维,为解决问题增加一条新的思路。
      

基于以上原因,我有了写这篇文章的打算,一是回顾总结这几天的研究和实验,加深印象,二是和朋友们分享此算法以及我的一些经验。
      

这篇文章将首先从直观的角度,介绍Monte-Carlo算法,然后介绍算法基本原理及数理基础,最后将会和大家分享几个基于Monte-Carlo方法的有意思的实验。所以程序将使用C#实现。
      

阅读本文需要有一些概率论、数理统计、微积分和计算复杂性的基本知识,不过不用太担心,我将尽量避免过多的数学描述,并在适当的地方对于用到的数学知识进行简要的说明。



2

Monte-Carlo算法引导


      首先,我们来看一个有意思的问题:在一个1平方米的正方形木板上,随意画一个圈,求这个圈的面积。
      我们知道,如果圆圈是标准的,我们可以通过测量半径r,然后用 S = pi * r^2 来求出面积。可是,我们画的圈一般是不标准的,有时还特别不规则,如下图是我画的巨难看的圆圈。

图1、不规则圆圈

      

显然,这个图形不太可能有面积公式可以套用,也不太可能用解析的方法给出准确解。不过,我们可以用如下方法求这个图形的面积:
      

假设我手里有一支飞镖,我将飞镖掷向木板。并且,我们假定每一次都能掷在木板上,不会偏出木板,但每一次掷在木板的什么地方,是完全随机的。即,每一次掷飞镖,飞镖扎进木板的任何一点的概率的相等的。这样,我们投掷多次,例如100次,然后我们统计这100次中,扎入不规则图形内部的次数,假设为k,那么,我们就可以用 k/100 * 1 近似估计不规则图形的面积,例如100次有32次掷入图形内,我们就可以估计图形的面积为0.32平方米。

      

以上这个过程,就是Monte-Carlo算法直观应用例子。
      

非形式化地说,Monte-Carlo算法泛指一类算法。在这些算法中,要求解的问题是某随机事件的概率或某随机变量的期望。这时,通过“实验”方法,用频率代替概率或得到随机变量的某些数字特征,以此作为问题的解。
      

上述问题中,如果将“投掷一次飞镖并掷入不规则图形内部”作为事件,那么图形的面积在数学上等价于这个事件发生的概率(稍后证明),为了估计这个概率,我们用多次重复实验的方法,得到事件发生的频率 k/100 ,以此频率估计概率,从而得到问题的解。

      

从上述可以看出,Monte-Carlo算法区别于确定性算法,它的解不一定是准确或正确的,其准确或正确性依赖于概率和统计,但在某些问题上,当重复实验次数足够大时,可以从很大概率上(这个概率是可以在数学上证明的,但依赖于具体问题)确保解的准确或正确性,所以,我们可以根据具体的概率分析,设定实验的次数,从而将误差或错误率降到一个可容忍的程度。
      

上述问题中,设总面积为S,不规则图形面积为s,共投掷n次,其中掷在不规则图形内部的次数为k。根据伯努利大数定理,当试验次数增多时,k/n依概率收敛于事件的概率s/S。下面给出严格证明:


      

上述证明从数学上说明用频率估计不规则图形面积的合理性,进一步可以给出误差分析,从而选择合适的实验次数n,以将误差控制在可以容忍的范围内,此处从略。

      

从上面的分析可以看出,Monte-Carlo算法虽然不能保证解一定是准确和正确,但并不是“撞大运”,其正确性和准确性依赖概率论,有严格的数学基础,并且通过数学分析手段对实验加以控制,可以将误差和错误率降至可容忍范围。



3

Monte-Carlo算法的数理基础


这一节讨论Monte-Carlo算法的数理基础。
      

首先给出三个定义:优势,一致,偏真。这三个定义在后面会经常用到。

      1) 设p为一个实数,且0.5<p<1。如果一个Monte-Carlo方法对问题任一实例的得到正确解的概率不小于p,则该算法是p正确的,且p-0.5叫做此算法的优势。

      2) 如果对于同一实例,某Monte-Carlo算法不会给出不同的解,则认为该算法时一致的。

      3) 如果某个解判定问题的Monte-Carlo算法,当返回true时是一定正确的。则这个算法时偏真的。注意,这里没有定义“偏假”,因为“偏假”和偏真是等价的。因为只要互换算法返回的true和false,“偏假”就变成偏真了。

      

下面,我们讨论Monte-Carlo算法的可靠性和误差分析。
      

总体来说,适用于Monte-Carlo算法的问题,比较常见的有两类。一类是问题的解等价于某事件概率,如上述求不规则图形面积的问题;另一类是判定问题,即判定某个命题是否为真,如主元素存在性判定和素数测试问题。

      

先来分析第一类。对于这类问题,通常的方法是通过大量重复性实验,用事件发生的频率估计概率。之所以能这样做的数学基础,是伯努利大数法则:事件发生的频率依概率收敛于事件的概率p。这个法则从数学生严格描述了频率的稳定性,直观意义就是当实验次数很大时,频率与概率偏差很大的概率非常小。此类问题的误差分析比较繁杂,此处从略。有兴趣的朋友可以参考相关资料。

     

接着,我们分析第二类问题。这里,我们只关心一致且偏真的判定问题。下面给出这类问题的正确率分析:

      

由以上分析可以看到,对于一致偏真的Monte-Carlo算法,即使调用一次得到正确解的概率非常小,通过多次调用,其正确率会迅速提高,得到的结果非常可靠。例如,对一个q为0.5的问题,假设p仅为0.01,通过调用1000次,其正确率约为0.9999784,几乎可以认为是绝对准确的。重要的是,使用Monte-Carlo算法解判定问题,其正确率不随问题规模而改变,这就使得仅需要损失微乎其微的正确性,就可以将算法复杂度降低一个数量级,在后面中可以看到具体的例子。



4

应用实例一

   使用Monte-Carlo算法计算定积分


计算定积分是金融、经济、工程等领域实践中经常遇到的问题。通常,计算定积分的经典方法是使用Newton-Leibniz公式:


这个公式虽然能方便计算出定积分的精确值,但是有一个局限就是要首先通过不定积分得到被积函数的原函数。有的时候,求原函数是非常困难的,而有的函数,如f(x) = (sinx)/x,已经被证明不存在初等原函数,这样,就无法用Newton-Leibniz公式,只能另想办法。
      

下面就以f(x) = (sinx)/x为例介绍使用Monte-Carlo算法计算定积分的方法。首先需要声明,f(x) = (sinx)/x在整个实数域是可积的,但不连续,在x = 0这一点没有定义。但是,当x趋近于0其左右极限都是1。为了严格起见,我们补充定义当x = 0时f(x) = 1。另外为了需要,这里不加证明地给出f(x)的一些性质:补充x = 0定义后,f(x)在负无穷到正无穷上连续、可积,并且有界,其界为1,即|f(x)| <= 1,当且仅当x = 0时f(x) = 1。

      

下面开始介绍Monte-Carlo积分法。为了便于比较,在本节我们除了介绍使用Monte-Carlo方法计算定积分外,同时也探讨和实现数值计算中常用的插值积分法,并通过实验结果数据对两者的效率和精确性进行比较。


1.插值积分法


我们知道,对于连续可积函数,定积分的直观意义就是函数曲线与x轴围成的图形中,y>0的面积减掉y<0的面积。那么一种直观的数值积分方法是通过插值方法,其中最简单的是梯形法则:用以f(a)和f(b)为底,x轴和f(a)、f(b)连线为腰组成的梯形面积来近似估计积分。如下图所示。

图2、梯形插值

      

如图2所示,蓝色部分是x1到x2积分的精确面积,而在梯形插值中,用橙色框所示的梯形面积近似估计积分值。
      

显然,梯形法则的效果一般,而且某些情况下偏差很大,于是,有人提出了一种改进的方法:首先将积分区间分段,然后对每段计算梯形插值再加起来,这样精度就大大提高了。并且分段越多,精度越高。这就是复化梯形法则。

      

除了梯形插值外,还有许多插值积分法,比较常见的有Sinpson法则,当然对应的也有复化Sinpson法则。下面给出四种插值积分的公式:



下面是四种插值积分法的程序代码,用C#编写。



2.Monte-Carlo积分法


我们知道,求定积分的直观意义就是求面积,所以,用Monte-Carlo求积分的原理就是通过模拟统计方法求解面积。即通过向特定区域随机产生大量点,然后统计点落在函数区域内的频率,以此频率估计面积,从而得到积分值。下面给出Monte-Carlo求取积分的算法程序。


3.积分法的测试与比较


 下面对各种积分方法进行测试,对sinx/x在[1,2]区间上进行定积分。其中,我们分别对复化梯形和复化Sinpson法则做分段为10,10000,和10000000的积分测试。另外,对Monte-Carlo法也投点数也分为10,10000,和10000000。测试结果如下:

图3、积分法测试结果

      

为了分析偏差,我们必须给出一个精确值。但是现在我手头没有这个积分的精确值,不过1000万次的梯形法则和Sinpson法则已经精确度很高了,所以这里就以0.65932985作为基本,进行误差分析。下面给出分析结果:

表1、积分方法实验结果

      

首先看时间效率。当频度较低时,各种方法没有太多差别,但在1000万级别上复化梯形和复化Sinpson相差不大,而Monte-Carlo算法的效率快一倍。
      

而从准确率分析,当频度较低时,几种方法的误差都很大,而随着频度提高,插值法要远远优于Monte-Carlo算法,特别在1000万级别时,Monte-Carlo法的相对误差是插值法的的近万倍。总体来说,在数值积分方面,Monte-Carlo方法效率高,但准确率不如插值法。



5

应用实例二

  在O(n)复杂度内判定主元素

这次,我们看一个判定问题。问题是这样的:在一个长度为n的数组中,如果有超过[n/2]的元素具有相同的值,那么具有这个值的元素叫做数组的主元素。现在要求给出一种算法,在O(n)时间内判定给定数组是否存在主元素。

      

如果采用确定性算法,由于最坏情况下要搜索n/2次,而每次要比较的次数为O(n)量级,这样,算法的复杂度就是O(n^2),不可能在O(n)时间内完成。所以我们只好换一种思路:不是要一个一定正确的结果,而只需要结果在很大概率上正确就行。我们可以这样做:

图4、Monte-Carlo法判定主元素

      

上述算法,就是用Monte-Carlo思想求解主元素判定问题的过程。由于阈值N是一个给定的常数,不随规模变化而变化,所以这个算法的时间复杂度为O(n),符合题设要求。但这个算法给出的解并不是100%正确的,正确率和N有关。N设得过大,影响效率,N太小,正确率太低,那么到底N设多大合适呢。这就要对算法进行概率分析。

      

首先,这个算法是一致且偏真的,证明很简单,这里从略。所以,如果数组中不存在主元素,则结果一定正确,而如果存在,调用一次得到正确结果的概率不低于1/2。由于偏真,在N次调用中只要返回一次True,就可以认为得到正确结果,所以,调用N此得到正确结果的概率不低于1 – (1/2)^N,可以看到,随着N的增大,这个概率增加很快,只要调用10次,正确率就可以达到99.9%,重要的是,这个正确率和规模无关,即使数组的元素有1千万亿,只需调用10次,正确率依然是99.9%,这就体现出在数组很大时,Monte-Carlo方法的优势。

      

下面是使用Monte-Carlo算法进行主元素测试的C#程序示例。


程序很简单,不做赘述。下面测试这个算法。我们分别将阈值设为1、3、10,并且在每个阈值下测试100次,看看这个算法的准确率如何。测试数组是[ 4, 5, 8, 1, 8, 4, 9, 2, 2, 2, 2, 2, 5, 7, 8, 2, 2, 2, 2, 2, 1, 0, 9, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 4, 7, 8, 2, 2, 2, 2, 2, 0, 1, 2, 2, 2, 2, 2 ],其中存在主元素2。下面是测试结果:

图5、Monte-Carlo算法判定主元素实验结果

 

测试数组有49个元素,主元素2有29个,比率为59%。从测试结果可以看出,即使阈值为1,正确率也高达84%,而仅仅为3的阈值就使正确率升高到98%,阈值为10时,100次测试全部正确。虽然理论上来说,阈值为10时有0.41^10=0.013%的概率给出错误判断,但是笔者多次试验,还没有在阈值为10时得到错误结果。所以,Monte-Carlo方法求解判定问题,不论从理论上还是实践中,都是不错的方法。

      

另外一个与判定主元素类似的应用是素数判定问题,我们知道,对于寻找上百位的大素数,完全测试在时间效率上时不允许的。于是,结合费马小定理使用Monte-Carlo法进行素数判定,是广泛使用的方法。具体这里不再详述,感兴趣的朋友可以参考相关资料。



6

应用实例三

   分布未知的概率密度函数模拟


现在我们来看看Monte-Carlo算法的第三种应用:模拟。在这种应用中,不再是用Monte-Carlo算法求解问题,而是用来模拟难以解析描述的东西。问题是这样的:

      

这个问题是实验室一个师兄在开发Six Sigma软件开发过程管理工具时遇到的一个实际需求,最终Y的概率密度函数将被用来计算分位点,从而进行过程控制。其中X可能是正态分布(高斯分布)、泊松分布、均匀分布或指数分布等。将多个不同分布的概率密度函数相加,得到的Y的分布式很难解析表示出来的,但如果是为了计算分位点,我们可以采取这样一个策略:对于每一个X,产生若干符合其分布的点,带入公式就得到若干符合Y分布的点,然后分段计算频率,从而模拟出Y的分布,这些模拟点也可以用于分位点计算。这就是Monte-Carlo模拟的思想。

      

下面我们实现这个算法,这里的X我们仅给出最常用的正态分布,如果要实现其他分布,只要编写相应的随机点发生器就可以了。由于C#中只能产生符合均匀分布的随机数,所以我们需要一种算法,将均匀分布的随机数转为正态分布随机数。这种算法很多,Marc Brysbaert在1991年发表的《Algorithms for randomness in the behavioral sciences: A tutorial》一文中,共总结了5种将均匀分布随机数转为正态分布的随机数的算法,这里笔者用到的是Knuth在1981年提出的一种算法。这个算法是将符合u(0,1)均匀分布的随机点转换为符合N(0,1)标准正态分布的随机点p,由概率知识可知,要转为符合N(e,v)的一般正态分布,只需进行p*v+e即可。下面是这个算法:

      

下面是根据这个算法,使用C#编写的正态分布随机点发生器:

   

接着是利用这个正态分布发生器获得X的随机值,并计算出Y的随机值的代码。也就是Y的随机点发生器:


这样,我们就可以产生任意多个符合Y分布的随机点,从而借此模拟Y的概率密度分布。

      

接着,我们测试一下这个模拟程序的效果,首先我们将初始值设为仅有一个符合标准正态分布的X,这样Y=X,我们看看直接模拟一个标准正态分布的效果。这里,我们产生100万个随机点。

图6、使用Monte-Carlo算法模拟标准正态分布

      

可以看到,模拟效果基本令人满意。接下来,我们实际应用这个程序模拟一个分布未知的Y,其中Y = 15 + 2*N(2,8) + 5*N(-10,9) + 7*N(0,0.5)。模拟结果如下:

图7、使用Monte-Carlo算法模拟未知分布

      

有了符合Y分布的大量随机点以及频率统计,就可以随心所欲绘出分布模拟图,并进行分位点计算。这样就用Monte-Carlo算法解决了本节开头提到的问题。



7

总结


本文首先通过一个不规则图形面积计算的例子直观介绍了Monte-Carlo算法,然后给出了Monte-Carlo算法在应用过程中需要了解的数理基础。然后大篇幅介绍了三个应用:计算、判定和模拟。
      

总体来说,当需要求解的问题依赖概率时,Monte-Carlo方法是一个不错的选择。但这个算法毕竟不是确定性算法,在应用过程中需要冒一定“风险”,这就要求不能滥用这个算法,在应用过程中,需要对其准确率或正确率进行数理分析,合理设计实验,从而得到良好的结果,并将风险控制在可容忍的范围内。
      

其实,不确定性算法不只Monte-Carlo一种,Sherwood算法、Las Vegas算法和遗传算法等也是经典的不确定算法。在很多问题上,不确定性算法具有很好大的应用价值。有兴趣的朋友可以参考相关资料。


参考文献

      [1] 孙海燕,周梦等 著,应用数理统计。北京航空航天大学出版社,2008.8
      [2] 盛骤,谢式千,潘承毅 著,概率论与数理统计。高等教育出版社,2006.12
      [3] David Kincaid,WardCheney 著,王国荣等 译,数值分析(原书第三版)。机械工业出版社,2005.9
      [4] Thomas H. Cormen等 著,算法导论(第二版,影印版)。高等教育出版社,2002.5
      [5] 王晓东 著,计算机算法设计与分析。电子工业出版社,2001.1
      [6] Marc Brysbaert,Algorithms for randomness in the behavioral sciences: A tutorial。Behavior Research Methods, Instruments, & Computers 1991, 23 (1) 45-60
      [7] Patrick Smacchia 著,施凡等 译,C#和.NET2.0 平台、语言与框架。2008.1
      [8] Google。www.google.com
      [9] Wikipedia。www.wikipedia.org


Copyright © 2023 All Rights Reserved 版权所有 北京物流信息联盟