1、 1 / 15蒙特卡罗方法（Monte Carlo method） 蒙特卡罗方法概述蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。 蒙特卡罗方法的提出蒙特卡罗方法于 20 世纪 40 年代美国在第二次世界大战中研制原子弹的“曼哈顿计划” 计划的成员 S.M.乌拉姆和 J.冯诺伊曼首先提出。数学家 冯诺伊曼用驰名世界的赌城摩纳哥的
2、 Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已经存在。1777 年，法国 Buffon 提出用投针实验的方法求圆周率 。这被认为是蒙特卡罗方法的起源。 蒙特卡罗方法的基本思想 Monte Carlo 方法的基本思想很早以前就被人们所发现和利用。早在 17 世纪，人们就知道用事件发生的“频率”来决定事件的“ 概率”。19 世纪人们用投针试验的方法来决定圆周率 。本世纪 40 年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。考虑平面上的一个边长为 1 的正方形及其内部的一个形状不规则的“
3、图形”，如何求出这个“ 图形”的面积呢？Monte Carlo 方法是这样一种 “随机化”的方法：向该正方形“随机地”投掷 N 个点，有 M 个点落于“ 图形”内，则该“图形” 的面积近似为 M/N。可用民意测验来作一个不严格的比喻。民意测验的人不是征询每一个登记选民的意见，而是通过对选民进行小规模的抽样调查来确定可能的优胜者。其基本思想是一样的。 科技计算中的问题比这要复杂得多。比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”(Curse of Dimensional
4、ity)，传统的数值方法难以对付（即使使用速度最快的计算机）。Monte Carlo 方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减” 技巧。 另一类形式与 Monte Carlo 方法相似，但理论基础不同的方法 “拟蒙特卡罗方法”(Quasi-Monte Carlo 方法)近年来也获得迅速发展。我国数学家华罗庚、王元提出的“ 华王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为 Low Discrepancy Sequences)代替 Mo
5、nte Carlo 方法中的随机数序列。对某些问题该方法的实际速度一般可比 Monte Carlo 方法提出高数百倍，并可计算精确度。 蒙特卡罗方法的基本原理由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡罗法正是基于此思路进行分析的。 设有统计独立的随机变量 Xi(i=1，2，3，k)，其对应的概率密度函数分别为 fx1，fx2，fxk ，功能函数式为 Z=g(x1，x2
6、，xk)。 首先根据各随机变量的相应分布，产生 N 组随机数 x1，x2，xk 值，计算功能函数值 Zi=g(x1，x2，xk)(i=1，2，N)，若其中有 L 组随机数对应的功能函数值 Zi0，则当 N时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。 2 / 15从蒙特卡罗方法的思路可看出，该方法回避了结构可靠度分析中的数学困难，不管状态函数是否非线性、随机变量是否非正态，只要模拟的次数足够多，就可得到一个比较精确的失效概率和可靠度指标。特别在岩土体分析中，变异系数往往较大，与 JC 法计算的可靠指标相比，结果更为精确，并且由于思路简单易于编制程序。 蒙特卡罗方法在数学
7、中的应用通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一种有效的求出数值解的方法。一般蒙特卡罗方法在数学中最常见的应用就是蒙特卡罗积分。 蒙特卡罗方法的应用领域蒙特卡罗方法在金融工程学，宏观经济学，生物医学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。 蒙特卡罗方法的工作过程在解决实际问题的时候应用蒙特卡罗方法主要有两部分工作： 1. 用蒙特卡罗方法模拟某一过程时，需要产生各种概率分布的 随机变量。 2. 用统计方法把模型的数字特征估计出来，从而得到实际问题的数
8、值解。 蒙特卡罗方法分子模拟计算的步骤使用蒙特卡罗方法进行分子模拟计算是按照以下步骤进行的： 1. 使用随机数发生器产生一个随机的分子构型。 2. 对此分子构型的其中粒子坐标做无规则的改变，产生一个新的分子构型。 3. 计算新的分子构型的能量。 4. 比较新的分子构型于改变前的分子构型的能量变化，判断是否接受该构型。 若新的分子构型能量低于原分子构型的能量，则接受新的构型，使用这个构型重复再做下一次迭代。 若新的分子构型能量高于原分子构型的能量，则計算玻尔兹曼因子，并产生一个随机数。 若这个随机数大于所计算出的玻尔兹曼因子，则放弃这个构型，重新计算。 若这个随机数小于所计算出的玻尔兹曼因子，则
9、接受这个构型，使用这个构型重复再做下一次迭代。 5. 如此进行迭代计算，直至最后搜索出低于所给能量条件的分子构型结束。 蒙特卡罗模型的发展运用 从理论上来说，蒙特卡罗方法需要大量的实验。实验次数越多，所得到的结果才越精确。以上 Buffon 的投针实验为例、历史上的记录如下表 1。 3 / 15从表中数据可以看到，一直到公元 20 世纪初期，尽管实验次数数以千计，利用蒙特卡罗方法所得到的圆周率值，还是达不到公元 5 世纪祖冲之的推算精度。这可能是传统蒙特卡罗方法长期得不到推广的主要原因。计算机技术的发展，使得蒙特卡罗方法在最近 10 年得到快速的普及。现代的蒙特卡罗方法，已经不必亲自动手做实验
10、，而是借助计算机的高速运转能力，使得原本费时费力的实验过程，变成了快速和轻而易举的事情。它不但用于解决许多复杂的科学方面的问题，也被项目管理人员经常使用。 借助计算机技术，蒙特卡罗方法实现了两大优点： 一是简单，省却了繁复的数学报导和演算过程，使得一般人也能够理解和掌握； 二是快速。简单和快速，是蒙特卡罗方法在现代项目管理中获得应用的技术基础。 蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省，这些是用该方法处理大型复杂问题时的优势。因此，随着电子计算机的发展和科学技术问题的日趋复杂
11、，蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学 、信息科学 、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。 项目管理中蒙特卡罗模拟方法的一般步骤 项目管理中蒙特卡罗模拟方法的一般步骤是： 1、对每一项活动，输入最小、最大和最可能估计数据，并为其选择一种合适的先验分布模型； 2、计算机根据上述输入，利用给定的某种规则，快速实施充分大量的随机抽样； 3、对随机抽样的数据进行必要的数学计算，求出结果； 4、对求出的结果进行统计学
12、处理，求出最小值、最大值以及数学期望值和单位标准偏差； 5、根据求出的统计学处理数据，让计算机自动生成概率分布曲线和累积概率曲线(通常是基于正态分布的概率累积 S 曲线)； 6、依据累积概率曲线进行项目风险分析。 非权重蒙特卡罗积分非权重蒙特卡罗积分，也称确定性抽样，是对被积函数变量区间进行随机均匀抽样，然后对被抽样点的函数值求平均，从而可以得到函数积分的近似值。此种方法的正确性是基于概率论的中心极限定理。当抽样点数为 m 时，使用此种方法所得近似解的统计误差恒为 1 除于根号 M,不随积分维数的改变而改变。因此当积分维度较高时，蒙特卡罗方法相对于其他数值解法更优。 4 / 15运用 SAS
13、进行 Monte Carlo 蒙特卡罗模拟简要介绍运用 SAS 进行 Monte Carlo 蒙特卡罗模拟（第一弹）：运用 SAS 进行 Monte Carlo 蒙特卡罗模拟简要介绍1 什么是蒙特卡罗模拟蒙特卡罗模拟的定义：the use of random sampling techniques and often the use of computer simulation to obtain approximate solutions to mathematical or physical problems especially in terms of a range of values
14、 each of which has a calculated probability of being the solution (Merriam-Webster, Inc., 1994, pp. 754-755)蒙特卡罗模拟的基本思想当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验” 的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。（资料来源：http:/zh.wikipedia.org/wiki/%E8%92%99%E7%89%B9%C2%B7%E5%8D%A1%E7%BD%97%E6%96%B
15、9%E6%B3%95）蒙特卡罗模拟的试验过程：计算机模拟试验过程，就是将试验过程转化为数学问题，在计算机上实现。在解决实际问题的时候应用蒙特卡罗模拟主要有两部分工作：用蒙特卡罗模拟某一过程时，需要产生各种概率分布的随机变量，然后用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。（资料来源：http:/zh.wikipedia.org/wiki/%E8%92%99%E7%89%B9%C2%B7%E5%8D%A1%E7%BD%97%E6%96%B9%E6%B3%95）2 举例：用蒙特卡罗模拟掷骰子一个均匀正方形的六面体，它在每一面刻有 1 到 6 的六个数码，如果要想知道掷骰子时某一面
16、向上机会的多少，或者掷两次，求两次的和的值出现的概率，既可以用数学理论来证明，也可以用多次投掷方法来验证。只要投掷的数目足够多，就一定能证明每一面朝上的机会是 1/6。如果根据概率知识，我们知道骰子每一面出现的概率为 1/6，因此，当掷两次时，因为每一次掷都是独立的，根据联合概率我们知道掷两次时的概率为1/6*1/6=1/36，那么当我们要求出掷两次，其和为 7 的概率时，因为 1+6，2+5，3+4 ，4+3，5+2 ，6+1 共六种可能，因为其概率为 6*1/36=1/6。如果我们不用概率理论，我们也可以用经验的方式得到结果，即用蒙特卡罗模拟每次掷骰子的结果。用蒙特卡罗模拟掷骰子的程序如下
17、：DATA DICE(KEEP=SUM) OUTCOMES(KEEP=OUTCOME);DO ROLL=1 TO 10000; * 掷 10000 次骰子;OUTCOME1=1+INT(6*RANUNI(123); * 第一次掷骰子的结果 ;OUTCOME2=1+INT(6*RANUNI(123); *第二次掷骰子的结果 ;SUM=OUTCOME1+OUTCOME2; * 两次掷骰子的结果的和 ;OUTPUT DICE; * 保存两次掷骰子的结果的和;5 / 15OUTCOME=OUTCOME1; OUTPUT OUTCOMES; * 保存第一次掷骰子的结果;OUTCOME=OUTCOME2;
18、 OUTPUT OUTCOMES; * 保存第二次掷骰子的结果;END;RUN;这里我们模拟 10000 次掷骰子的结果，得到这些结果的分布：PROC FREQ DATA=DICE; * 得到两次掷骰子的结果的和的分布;TABLE SUM;RUN;结果如下：Cumulative CumulativeSUM Frequency Percent Frequency Percent-2 299 2.99 299 2.993 534 5.34 833 8.334 811 8.11 1644 16.445 1177 11.77 2821 28.216 1374 13.74 4195 41.957 168
19、5 16.85 5880 58.508 1361 13.61 7241 72.419 1083 10.83 8324 83.2410 852 8.52 9176 91.7611 540 5.40 9716 97.1012 284 2.84 10000 100.00这里我们可以看到，其和为 7 的概率为 16.85，与概率理论计算的结果 1/6=16.67 相差无几。PROC FREQ DATA=OUTCOMES; *得到每一次掷骰子的结果的分布 ;TABLE OUTCOME;RUN;结果如下：Cumulative CumulativeOUTCOME Frequency Percent Freq
20、uency Percent-1 3298 16.49 3298 16.492 3367 16.84 6665 33.333 3362 16.81 10027 50.144 3372 16.86 13399 67.005 3341 16.71 16740 83.706 3260 16.30 20000 100.00这里我们可以看到每个 outcome 的 Percent 均在 1/6=16.67 附近，也即验证了骰子每一面朝上的机会是 1/6。3 蒙特卡罗模拟主要适用于哪些情况6 / 15刚才的例子中，我们可能会有疑问：用概率理论已经能很有效地得到结果了，为什么还要需要蒙特卡罗模拟呢。原因是刚才
21、的例子只是一个例子，蒙特卡罗模拟可以做很多事情。比如一般的概率理论都是有一些假设前提的，特别是对样本数据分布的假设，另一方面，如果我们手上的数据并不满足概率理论的假设条件，那么我们将得到有偏的结论。这时，我们就可以利用蒙特卡罗模拟，它是基于对样本分布特征的经验预测而非理论预测的，因此，当我们手上的数据并不能满足统计理论的假设时，蒙特卡罗模拟就是理论预测的一个很好的替代。那么蒙特卡罗模拟到底可以应用于哪些地方呢，一个主要应用是评估估计结果是否违反最初的假设，另一个应用是当一个样本没有理论分布时，我们可以用蒙特卡罗模拟来决定样本分布。蒙特卡罗方法的主要应用范围：蒙特卡罗方法所特有的优点，使得它的应
22、用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。4 蒙特卡罗模拟的优缺点优点：能够比较逼真地描述具有随机性质的事物的特点及物理实验过程：从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。受几何条件限制小：在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。收敛速度与问题的维数无关：使用蒙特卡罗方法时，抽取的子样总数 N 与
23、维数 s 无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。具有同时计算多个方案与多个未知量的能力：对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率
24、在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。误差容易确定：对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。程序结构简单，易于实现：在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实
25、现。缺点：收敛速度慢：蒙特卡罗模拟一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。误差具有概率性：由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。 对于蒙特卡罗模拟更多的内容请查阅相关文献。 参考资料 Xitao Fan, etcMonte Carlo Studies: A Guide for Quantitative Researchers. SAS Institute Inc.,2002蒙特卡罗方法概述：SAS 进行 Monte Carlo 蒙特卡罗模拟的基本步骤运用 SAS 进行 Monte Carlo 蒙
26、特卡罗模拟（第二弹）：SAS 进行 Monte Carlo 蒙特卡罗模拟的基本步骤本文未经作者同意严禁转载7 / 151 将要解决的事项转换成蒙特卡罗模拟可以解决的问题蒙特卡罗模拟适用于解决哪些问题？一般来说，与抽样分布相关的统计都适合蒙特卡罗模拟，特别是当这些问题没有可靠的理论解，例如理论假设有偏差，或者与统计量相关的理论不强，或者理论根本不存在时。2 设计蒙特卡罗模拟我们通过一个实例来说明如何应用蒙特卡罗模拟。在本文中，我们分析一下影响样本中变量间相关系数变化的原因，在这里，我们为了简化过程，仅考虑样本的大小这一个因素，而不考虑其它的因素，例如数据是否正态分布等等。在决定只考虑样本大小对相
27、关系数的影响后，我们可以模拟一些样本大小分别为 n=10, 20, 40, 100 的数据集，这些数据集是从一个相关性为零的总体里分别随机抽取 10, 20, 40, 100 个组成的，然后再计算这些样本数据集中的变量的相关性。这里，我们还有一个重要的问题：就是对于每一类指定大小的样本，我们需要计算多少个相关系数值来验证我们的假设（即总体中两个变量的相关系数为零） 。如果相关系数的数量太少，则可能产生有偏的结论，这里我们取的数量为 2000，即：每个指定大小的样本随机产生 2000 次，然后计算每一次产生的样本的变量的相关系数。因此，蒙特卡罗模拟中的相关系数的样本分布主要有以下设计：四类大小的
28、样本：10, 20, 40, 100对于每一类大小的样本，都从相关系数为零的总体中随机抽取 2000 次，得到 2000 个样本中的两个变量相关系数。 本节中的 SAS 代码实现如下：LIBNAME CORR C:CORR_EG;%LET NO_SMPL=2000; * 记录每类样本大小情况下，其随机样本个数的宏变量;%MACRO CORR_RDM;%DO A = 1 %TO 4; * 四类样本大小情况: 10, 20, 40, 100;%IF %LET SMPLSIZE=10; %END;%IF %LET SMPLSIZE=20; %END;%IF %LET SMPLSIZE=50; %EN
29、D;%IF %LET SMPLSIZE=100; %END;%DO B=1 %TO *每类样本大小情况下，生成DATA DAT; *生成两个随机变量,其相关性为零;DO I=1 TO X=RANNOR(0);Y=RANNOR(0);OUTPUT;END;PROC CORR DATA=DAT NOPRINT OUTP=PEARSON; *对每个随机样本，计算两个变量的相关系数，并保存在数据集 PEARSON 中;VAR X Y;RUN;DATA PEARSON; SET PEARSON; *将刚才得到的相关系数数据集增加样本大小这个变量;SMPLSIZE=IF _NAME_=X;CORR=Y;8
30、 / 15KEEP CORR SMPLSIZE; *仅保留相关系数和样本大小两个变量;PROC APPEND BASE=CORR.COR_RDM; *将得到的相关系数和样本大小的数据集全部放入CORR.COR_RDM 数据集中;%END;%END;%MEND CORR_RDM;%CORR_RDM;RUN; QUIT;这里，数据集 CORR.COR_RDM 保存了所有的样本大小的 2000 个相关系数。3 生成样本数据这里要强调一下，生成的样本数据是要符合蒙特卡罗模拟需求的，否则得到的结论不足以验证最初的假设。根据蒙特卡罗模拟的复杂程度，生成样本数据主要有三大步：I 生成一个 SAS 里已有分布
31、的数据，例如标准正态分布（RANNOR） 。II 将第 I 步生成的数据按需要的形状进行转换，例如将标准正态分布转换成其它正态分布III 已知变量相关性的多变量分布数据3.1 生成一个 SAS 里已有分布的数据这里，我们需要生成两个相互独立的变量的数据集，因此，我们可以生成两个独立的正态分布随机变量，通过SAS 的 RANNOR 函数可以实现随机正态变量的生成。在上一节的程序中，我们生成了两个变量，分别是 X和 Y，两者的均值为零，方差为 1，即都符合标准正态分布 N(0,1)。SAS 数据生成我们将在以后的文章中更详细地介绍。3.2 将第 I 步生成的数据按需要的形状进行转换在很多情况下，我
32、们需要将生成的符合某些标准分布的样本数据转换成其它分布的数据，此类转换主要有两类：一类是将样本数据转换成特定的均值和方差的数据，另一类是将数据转换成指定的形状（例如特定的峰度和偏度等） 。对于第一类转换，我们可以直接通过公式进行线性转换即可，因为此类转换并不改变样本分布的正态性，公式为：X(new)=X*SD(new)+Mean(new)，这里 X(new)即为转换后的变量，SD(new)为转换后变量的新的标准差，Mean(new) 为转换后变量的新的均值。转换后的变量 X(new)和原变量 X 有着同样的峰度和偏度。第二类转换则相对来说比较复杂，它需要一些统计模拟来生成特定峰度和偏度的非正态
33、分布数据。此类转换将在以后的文章中更详细地介绍。3.3 已知变量相关性的多变量分布数据前面两种转换只针对单个变量，而实现情况中，我们需要生成有着相关性的多变量的数据此类变量的生成我们将在后面的文章中介绍。 4 运用统计手段生成我们需要的统计量在第 2 节中，我们通过 PROC CORR 过程步计算每个样本的相关性，并将其保存在 PEARSON 数据集中，并最终汇总到数据集 CORR.COR_RDM 中。这样，我们就得到了 4 类样本大小各 2000 个相关系数的数据集，共 4*2000=8000 条数据。5 分析得到的统计量我们通过 PROC MEANS 来分析上一步得到的相关系数数据集，计算
34、出每一类样本大小的相关系数的均值、标准差、最大值、最小值等统计量。程序如下：DATA A; SET CORR.COR_RDM; *生成一个临时数据集来存放 CORR.COR_RDM，这样做的好处是保护9 / 15好 CORR.COR_RDM 数据集，使其被误操作等 ;PROC SORT; BY SMPLSIZE;PROC MEANS; BY SMPLSIZE; *更多 PROC MEANS 的相关内容详见前面的文章 ;VAR CORR;TITLE1 DESCRIPTIVE STATS FOR PEARSON RS BETWEEN TWO RANDOM VARIABLES;TITLE2 FOR
35、FOUR DIFFERENT SAMPLE SIZE CONDITIONS;TITLE3 *;RUN; QUIT;结果如下：我们也可以通过画出柱状图等来更直接地查看相关系数的分布情况，程序如下：DATA A; SET CORR.COR_RDM;PROC SORT; BY SMPLSIZE;AXIS1 LABEL=(HEIGHT=1.0 FONT=TRIPLEX) ORDER=(0 TO 20 BY 5)VALUE=(HEIGHT=1.0 FONT=TRIPLEX) MINOR=NONE;AXIS2 LABEL=(HEIGHT=1.0 FONT=TRIPLEX) VALUE=(HEIGHT=1.
36、0 FONT=TRIPLEX)MINOR=NONE;PATTERN COLOR=BLACK VALUE=X2;PROC GCHART DATA=A; BY SMPLSIZE; *更多 PROC GCHART 的内容请参考前面的文章;VBAR CORR/ TYPE=PERCENTMIDPOINTS= -.9 TO .9 BY .05RAXIS=AXIS1 MAXIS=AXIS2WIDTH=110 / 15SPACE=1;RUN; QUIT;结果如下：样本大小为 10 的相关系数结果：样本大小为 20 的相关系数结果：11 / 15样本大小为 50 的相关系数结果：样本大小为 100 的相关系数结
37、果：12 / 15有关 PROC GCHART 的内容请参考前面的文章。6 基于蒙特卡罗模拟的结果得出结论我们从均值的结果可以看到，四类样本大小的相关系数均值均在 0 附近，说明它们的总体均值为零，但是四类样本大小的相关系数的方差却有很大的差别，样本越大，方差越小，反之则方差越大。从图形结果我们能更明显地看出这种趋势。样本越大，相关系数越集中在 0 附近，样本越小，相关系数分布就越分散。因此，我们在计算样本相关性的时候，其样本的大小一定要达到一定的数量，其结果才可靠，否则，如果样本大小太小，很容易产生一些有偏的结论。参考资料：Xitao Fan, etcMonte Carlo Studies:
38、 A Guide for Quantitative Researchers. SAS Institute Inc.,2002常见概率分布及在 R 中的应用常见概率分布及在 R 中的应用R 提供工具来计算累计分布函数 p(cummulative distribution function CDF),概率密度函数 d 和分位数函数q，另外在各种概率分布前加 r 表示产生随机序列（R 这种直接在分布前面加前缀的语法太难读了，pt() 误以为还是一个函数，实际上的含义是 p(t()，为什么不写成这个格式呢？ 不过 t()返回什么好）常见概率分布离散型1.二项分布 Binomial distributi
39、on：binom二项分布指的是 N 重伯努利实验，记为 X b(n,p)，E(x)=np,Var(x)=np(1-p)。pbinom(q,size,prob)， q 是特定取值，比如 pbinom(8,20,0.2)指第 8 次伯努利实验的累计概率。size 指总的实验次数，prob 指每次实验成功发生的概率。dbinom(x,size,prob), x 同上面的 q 同含义。dfunction()对于离散分布来说结果是特定值的概率，对连续变量来说是密度（Density）。rbinom(n, size, prob)，产生 n 个 b(size,prob)的二项分布随机数 qbinom(p, s
40、ize, prob),quantile function 分位数函数。分位数：若概率 0Za)= 的实数。如t 分布的分位数表，自由度 f=20 和 =0.05 时的分位数为 1.7247。 这个定义指的是上侧 分位数。 分位数：实数 满足 0 2=1-F(2)=0.5 的数 2。qbinom 是上侧分位数，如qbinom(0.95,100,0.2)=27,指 27 之后 P(x=27)=0.95。即对于 b(100,0.2)为了达到 0.95 的概率至少需要 27次重复实验。2.负二项分布 negative binomial distribution （帕斯卡分布）nbinom掷骰子，掷到一
41、即视为成功。则每次掷骰的成功率是 1/6。要掷出三次一，所需的掷骰次数属于集合 3, 4, 5, 6, 。掷到三次一的掷骰次数是负二项分布的随机变量。dnbinom(4,3,1/6)=0.0334898，四次连续三次 1 的概率为这个数。概率函数为 f(k;r,p)=choose(k+r-1,r-1)*pr*(1-p)k, 当 r=1 时这个特例分布是几何分布。rnbinom(n,size,prob,mu) 其中 n 是需要产生的随机数个数，size 是概率函数中的 r，即连续成功的次数，prob 是单词成功的概率，mu 未知(mu 是希腊字母 的读音)13 / 153.几何分布 Geomet
42、ric Distribution,geomn 次伯努利试验，前 n-1 次皆失败，第 n 次才成功的机率。dgeom(x,prob),注意这里的 x 取值是 0:n，即dgeom(0,0.2)=0.2,以上的二项分布和负二项分布也是如此。ngeom(n,prob)4.超几何分布 Hypergeometric Distribution，hyper它描述了由有限个(m+n)物件中抽出 k 个物件，成功抽出指定种类的物件的次数（不归还）。概率：p(x) = choose(m, x) choose(n, k-x) / choose(m+n, k) for x = 0, , k。当 n=1 时，这是一个
43、0-1 分布即伯努利分布，当 n 接近无穷大时，超几何分布可视为二项分布。rhyper(nn,m,n,k),nn 是需要产生的随机数个数，m 是白球数（计算目标是取到 x 个白球的概率），n 是黑球数，k 是抽取出的球个数dhyper(x, m, n, k)。5.泊松分布 Poisson Distribution,poisp(x) = lambdax exp(-lambda)/x! for x = 0, 1, 2, . The mean and variance are E(X) = Var(X) = . x ()，泊松分布的参数 是单位时间(或单位面积)内随机事件的平均发生率.泊松分布适合于
44、描述单位时间内随机事件发生的次数。如某一服务设施在一定时间内到达的人数，电话交换机接到呼叫的次数，汽车站台的候客人数，机器出现的故障数，自然灾害发生的次数等等。rpois(n, lambda)，dpois(x,lambda)连续型6.均匀分布 Uniform Distribution，uniff(x) = 1/(max-min) for min = 0, a 0 and s 0。Gamma 分布中的参数 ，称为形状参数（shape parameter），即上式中的s， 称为尺度参数（scale parameter）上式中的 a。E(x)=s*a, Var(x)=s*a2. 当 shape=1/
45、2,scale=2 时，这样的 gamma 分布是自由度为 1 的开方分布http:/zh.wikipedia.org/wiki/File:Gamma_distribution_pdf.pngdgamma(x,shape,rate=1,scale=1/rate), 请注意 R 在这里提供的 rate 是 scale 尺度参数的倒数，如果dgamma(0,1,2)则表示 dgamma(0,shape=1,rate=2),而非 dgamma(0,shape=1,scale=2)14 / 15pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail
46、= TRUE,log.p = FALSE)qgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,log.p = FALSE)rgamma(n, shape, rate = 1, scale = 1/rate)9.指数分布 Exponential Distribution，exp指数分布可以用来表示独立随机事件发生的时间间隔，比如旅客进机场的时间间隔、中文维基百科新条目出现的时间间隔等等。记作 X Exponential（）。f(x) = lambda e(- lambda x) for x = 0.其中 lambda 0 是分
47、布的一个参数，常被称为率参数（rate parameter）. E(x)=1/,Var(x)=1/2dexp(x, rate = 1, log = FALSE)pexp(q, rate = 1, lower.tail = TRUE, log.p = FALSE)qexp(p, rate = 1, lower.tail = TRUE, log.p = FALSE)rexp(n, rate = 1)假设在公交站台等公交车平均 10 分钟有一趟车，那么每小时候有 6 趟车，即每小时出现车的次数 Exponential(1/6)我们可以产生 10 个这些随机数看看 rexp(10,1/6)60/(re
48、xp10,1/6)即为我们在站台等车的随机时间，如下：1 6.443148 24.337131 6.477096 2.824638 15.184945 14.5949037 7.133842 8.222400 42.609784 15.182827可以看见竟然有一个 42.6 分钟的随机数出现，据说这种情况下你可以投诉上海的公交公司。不过 x 符合指数分布，1/x 还符合指数分布吗？pexp(6,1/6)=0.6321206, 也就是说这种情况下只有 37%的可能公交车会 10 分钟以内来。按照以上分析一个小时出现的公交车次数应该不符合指数分布。10.卡方分布(non-central)Chi-Squared Distribution，chisq它广泛的运用于检测数学模型是否适合所得的数据，以及数据间
