蒙特卡罗方法.docx

蒙特卡罗方法（Monte Carlo method） 蒙特卡罗方法概述蒙特卡罗方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。 蒙特卡罗方法的提出蒙特卡罗方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯诺伊曼首先提出。数学家冯诺伊曼用驰名世界的赌城摩纳哥的Monte Carlo来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡罗方法就已经存在。1777年，法国Buffon提出用投针实验的方法求圆周率。这被认为是蒙特卡罗方法的起源。 蒙特卡罗方法的基本思想 Monte Carlo方法的基本思想很早以前就被人们所发现和利用。早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针试验的方法来决定圆周率。本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。 考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？Monte Carlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点，有M个点落于“图形”内，则该“图形”的面积近似为M/N。可用民意测验来作一个不严格的比喻。民意测验的人不是征询每一个登记选民的意见，而是通过对选民进行小规模的抽样调查来确定可能的优胜者。其基本思想是一样的。 科技计算中的问题比这要复杂得多。比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”(Curse of Dimensionality)，传统的数值方法难以对付（即使使用速度最快的计算机）。Monte Carlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。以前那些本来是无法计算的问题现在也能够计算量。为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。 另一类形式与Monte Carlo方法相似，但理论基础不同的方法“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)近年来也获得迅速发展。我国数学家华罗庚、王元提出的“华王”方法即是其中的一例。这种方法的基本思想是“用确定性的超均匀分布序列(数学上称为Low Discrepancy Sequences)代替Monte Carlo方法中的随机数序列。对某些问题该方法的实际速度一般可比Monte Carlo方法提出高数百倍，并可计算精确度。 蒙特卡罗方法的基本原理由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡罗法正是基于此思路进行分析的。 设有统计独立的随机变量Xi(i=1，2，3，k)，其对应的概率密度函数分别为fx1，fx2，fxk，功能函数式为Z=g(x1，x2，xk)。 首先根据各随机变量的相应分布，产生N组随机数x1，x2，xk值，计算功能函数值 Zi=g(x1，x2，xk)(i=1，2，N)，若其中有L组随机数对应的功能函数值Zi0，则当N时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。 从蒙特卡罗方法的思路可看出，该方法回避了结构可靠度分析中的数学困难，不管状态函数是否非线性、随机变量是否非正态，只要模拟的次数足够多，就可得到一个比较精确的失效概率和可靠度指标。特别在岩土体分析中，变异系数往往较大，与JC法计算的可靠指标相比，结果更为精确，并且由于思路简单易于编制程序。 蒙特卡罗方法在数学中的应用通常蒙特卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特卡罗方法是一种有效的求出数值解的方法。一般蒙特卡罗方法在数学中最常见的应用就是蒙特卡罗积分。 蒙特卡罗方法的应用领域蒙特卡罗方法在金融工程学，宏观经济学，生物医学，计算物理学(如粒子输运计算、量子热力学计算、空气动力学计算)等领域应用广泛。 蒙特卡罗方法的工作过程在解决实际问题的时候应用蒙特卡罗方法主要有两部分工作： 1. 用蒙特卡罗方法模拟某一过程时，需要产生各种概率分布的随机变量。 2. 用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。 蒙特卡罗方法分子模拟计算的步骤使用蒙特卡罗方法进行分子模拟计算是按照以下步骤进行的： 1. 使用随机数发生器产生一个随机的分子构型。 2. 对此分子构型的其中粒子坐标做无规则的改变，产生一个新的分子构型。 3. 计算新的分子构型的能量。 4. 比较新的分子构型于改变前的分子构型的能量变化，判断是否接受该构型。 l 若新的分子构型能量低于原分子构型的能量，则接受新的构型，使用这个构型重复再做下一次迭代。 l 若新的分子构型能量高于原分子构型的能量，则計算玻尔兹曼因子，并产生一个随机数。 l 若这个随机数大于所计算出的玻尔兹曼因子，则放弃这个构型，重新计算。 l 若这个随机数小于所计算出的玻尔兹曼因子，则接受这个构型，使用这个构型重复再做下一次迭代。 5. 如此进行迭代计算，直至最后搜索出低于所给能量条件的分子构型结束。 蒙特卡罗模型的发展运用 从理论上来说，蒙特卡罗方法需要大量的实验。实验次数越多，所得到的结果才越精确。以上Buffon的投针实验为例、历史上的记录如下表1。 从表中数据可以看到，一直到公元20世纪初期，尽管实验次数数以千计，利用蒙特卡罗方法所得到的圆周率值，还是达不到公元5世纪祖冲之的推算精度。这可能是传统蒙特卡罗方法长期得不到推广的主要原因。 计算机技术的发展，使得蒙特卡罗方法在最近10年得到快速的普及。现代的蒙特卡罗方法，已经不必亲自动手做实验，而是借助计算机的高速运转能力，使得原本费时费力的实验过程，变成了快速和轻而易举的事情。它不但用于解决许多复杂的科学方面的问题，也被项目管理人员经常使用。 借助计算机技术，蒙特卡罗方法实现了两大优点： 一是简单，省却了繁复的数学报导和演算过程，使得一般人也能够理解和掌握； 二是快速。简单和快速，是蒙特卡罗方法在现代项目管理中获得应用的技术基础。 蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省，这些是用该方法处理大型复杂问题时的优势。因此，随着电子计算机的发展和科学技术问题的日趋复杂，蒙特卡罗方法的应用也越来越广泛。它不仅较好地解决了多重积分计算、微分方程求解、积分方程求解、特征值计算和非线性方程组求解等高难度和复杂的数学计算问题，而且在统计物理、核物理、真空技术、系统科学 、信息科学 、公用事业、地质、医学，可靠性及计算机科学等广泛的领域都得到成功的应用。 项目管理中蒙特卡罗模拟方法的一般步骤 项目管理中蒙特卡罗模拟方法的一般步骤是： 1、对每一项活动，输入最小、最大和最可能估计数据，并为其选择一种合适的先验分布模型； 2、计算机根据上述输入，利用给定的某种规则，快速实施充分大量的随机抽样； 3、对随机抽样的数据进行必要的数学计算，求出结果； 4、对求出的结果进行统计学处理，求出最小值、最大值以及数学期望值和单位标准偏差； 5、根据求出的统计学处理数据，让计算机自动生成概率分布曲线和累积概率曲线(通常是基于正态分布的概率累积S曲线)； 6、依据累积概率曲线进行项目风险分析。 非权重蒙特卡罗积分非权重蒙特卡罗积分，也称确定性抽样，是对被积函数变量区间进行随机均匀抽样，然后对被抽样点的函数值求平均，从而可以得到函数积分的近似值。此种方法的正确性是基于概率论的中心极限定理。当抽样点数为m时，使用此种方法所得近似解的统计误差恒为 1除于根号M,不随积分维数的改变而改变。因此当积分维度较高时，蒙特卡罗方法相对于其他数值解法更优。 运用SAS进行MonteCarlo蒙特卡罗模拟简要介绍运用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 each of which has a calculated probability of being the solution (Merriam-Webster, Inc., 1994, pp. 754-755)蒙特卡罗模拟的基本思想当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。（资料来源：/wiki/%E8%92%99%E7%89%B9%C2%B7%E5%8D%A1%E7%BD%97%E6%96%B9%E6%B3%95）蒙特卡罗模拟的试验过程：计算机模拟试验过程，就是将试验过程转化为数学问题，在计算机上实现。在解决实际问题的时候应用蒙特卡罗模拟主要有两部分工作：用蒙特卡罗模拟某一过程时，需要产生各种概率分布的随机变量，然后用统计方法把模型的数字特征估计出来，从而得到实际问题的数值解。（资料来源：/wiki/%E8%92%99%E7%89%B9%C2%B7%E5%8D%A1%E7%BD%97%E6%96%B9%E6%B3%95）2 举例：用蒙特卡罗模拟掷骰子一个均匀正方形的六面体，它在每一面刻有1到6的六个数码，如果要想知道掷骰子时某一面向上机会的多少，或者掷两次，求两次的和的值出现的概率，既可以用数学理论来证明，也可以用多次投掷方法来验证。只要投掷的数目足够多，就一定能证明每一面朝上的机会是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。如果我们不用概率理论，我们也可以用经验的方式得到结果，即用蒙特卡罗模拟每次掷骰子的结果。用蒙特卡罗模拟掷骰子的程序如下：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; * 保存两次掷骰子的结果的和; OUTCOME=OUTCOME1; OUTPUT OUTCOMES; * 保存第一次掷骰子的结果; OUTCOME=OUTCOME2; OUTPUT OUTCOMES; * 保存第二次掷骰子的结果; END;RUN;这里我们模拟10000次掷骰子的结果，得到这些结果的分布：PROC FREQ DATA=DICE; * 得到两次掷骰子的结果的和的分布; TABLE SUM;RUN;结果如下： Cumulative Cumulative SUM Frequency Percent Frequency Percent - 2 299 2.99 299 2.99 3 534 5.34 833 8.33 4 811 8.11 1644 16.44 5 1177 11.77 2821 28.21 6 1374 13.74 4195 41.95 7 1685 16.85 5880 58.50 8 1361 13.61 7241 72.41 9 1083 10.83 8324 83.24 10 852 8.52 9176 91.76 11 540 5.40 9716 97.10 12 284 2.84 10000 100.00这里我们可以看到，其和为7的概率为16.85，与概率理论计算的结果1/6=16.67相差无几。PROC FREQ DATA=OUTCOMES; *得到每一次掷骰子的结果的分布; TABLE OUTCOME;RUN;结果如下： Cumulative CumulativeOUTCOME Frequency Percent Frequency Percent- 1 3298 16.49 3298 16.49 2 3367 16.84 6665 33.33 3 3362 16.81 10027 50.14 4 3372 16.86 13399 67.00 5 3341 16.71 16740 83.70 6 3260 16.30 20000 100.00这里我们可以看到每个outcome的Percent均在1/6=16.67附近，也即验证了骰子每一面朝上的机会是1/6。3 蒙特卡罗模拟主要适用于哪些情况刚才的例子中，我们可能会有疑问：用概率理论已经能很有效地得到结果了，为什么还要需要蒙特卡罗模拟呢。原因是刚才的例子只是一个例子，蒙特卡罗模拟可以做很多事情。比如一般的概率理论都是有一些假设前提的，特别是对样本数据分布的假设，另一方面，如果我们手上的数据并不满足概率理论的假设条件，那么我们将得到有偏的结论。这时，我们就可以利用蒙特卡罗模拟，它是基于对样本分布特征的经验预测而非理论预测的，因此，当我们手上的数据并不能满足统计理论的假设时，蒙特卡罗模拟就是理论预测的一个很好的替代。那么蒙特卡罗模拟到底可以应用于哪些地方呢，一个主要应用是评估估计结果是否违反最初的假设，另一个应用是当一个样本没有理论分布时，我们可以用蒙特卡罗模拟来决定样本分布。蒙特卡罗方法的主要应用范围：蒙特卡罗方法所特有的优点，使得它的应用范围越来越广。它的主要应用范围包括：粒子输运问题，统计物理，典型数学问题，真空技术，激光技术以及医学，生物，探矿等方面。4 蒙特卡罗模拟的优缺点优点：能够比较逼真地描述具有随机性质的事物的特点及物理实验过程：从这个意义上讲，蒙特卡罗方法可以部分代替物理实验，甚至可以得到物理实验难以得到的结果。用蒙特卡罗方法解决实际问题，可以直接从实际问题本身出发，而不从方程或数学表达式出发。它有直观、形象的特点。受几何条件限制小：在具有随机性质的问题中，如考虑的系统形状很复杂，难以用一般数值方法求解，而使用蒙特卡罗方法，不会有原则上的困难。收敛速度与问题的维数无关：使用蒙特卡罗方法时，抽取的子样总数N与维数s无关。维数的增加，除了增加相应的计算量外，不影响问题的误差。这一特点，决定了蒙特卡罗方法对多维问题的适应性。而一般数值方法，比如计算定积分时，计算时间随维数的幂次方而增加，而且，由于分点数与维数的幂次方成正比，需占用相当数量的计算机内存，这些都是一般数值方法计算高维积分时难以克服的问题。具有同时计算多个方案与多个未知量的能力：对于那些需要计算多个方案的问题，使用蒙特卡罗方法有时不需要像常规方法那样逐个计算，而可以同时计算所有的方案，其全部计算量几乎与计算一个方案的计算量相当。例如，对于屏蔽层为均匀介质的平板几何，要计算若干种厚度的穿透概率时，只需计算最厚的一种情况，其他厚度的穿透概率在计算最厚一种情况时稍加处理便可同时得到。另外，使用蒙特卡罗方法还可以同时得到若干个所求量。例如，在模拟粒子过程中，可以同时得到不同区域的通量、能谱、角分布等，而不像常规方法那样，需要逐一计算所求量。误差容易确定：对于一般计算方法，要给出计算结果与真值的误差并不是一件容易的事情，而蒙特卡罗方法则不然。根据蒙特卡罗方法的误差公式，可以在计算所求量的同时计算出误差。对干很复杂的蒙特卡罗方法计算问题，也是容易确定的。一般计算方法常存在着有效位数损失问题，而要解决这一问题有时相当困难，蒙特卡罗方法则不存在这一问题。程序结构简单，易于实现：在计算机上进行蒙特卡罗方法计算时，程序结构简单，分块性强，易于实现。 缺点：收敛速度慢：蒙特卡罗模拟一般不容易得到精确度较高的近似结果。对于维数少（三维以下）的问题，不如其他方法好。误差具有概率性：由于蒙特卡罗方法的误差是在一定置信水平下估计的，所以它的误差具有概率性，而不是一般意义下的误差。对于蒙特卡罗模拟更多的内容请查阅相关文献。参考资料Xitao Fan, etc.Monte Carlo Studies: A Guide for Quantitative Researchers. SAS Institute Inc.,2002蒙特卡罗方法概述：SAS进行MonteCarlo蒙特卡罗模拟的基本步骤运用SAS进行Monte Carlo蒙特卡罗模拟（第二弹）：SAS进行Monte Carlo蒙特卡罗模拟的基本步骤本文未经作者同意严禁转载1 将要解决的事项转换成蒙特卡罗模拟可以解决的问题蒙特卡罗模拟适用于解决哪些问题？一般来说，与抽样分布相关的统计都适合蒙特卡罗模拟，特别是当这些问题没有可靠的理论解，例如理论假设有偏差，或者与统计量相关的理论不强，或者理论根本不存在时。2设计蒙特卡罗模拟我们通过一个实例来说明如何应用蒙特卡罗模拟。在本文中，我们分析一下影响样本中变量间相关系数变化的原因，在这里，我们为了简化过程，仅考虑样本的大小这一个因素，而不考虑其它的因素，例如数据是否正态分布等等。在决定只考虑样本大小对相关系数的影响后，我们可以模拟一些样本大小分别为n=10, 20, 40, 100的数据集，这些数据集是从一个相关性为零的总体里分别随机抽取10, 20, 40, 100个组成的，然后再计算这些样本数据集中的变量的相关性。这里，我们还有一个重要的问题：就是对于每一类指定大小的样本，我们需要计算多少个相关系数值来验证我们的假设（即总体中两个变量的相关系数为零）。如果相关系数的数量太少，则可能产生有偏的结论，这里我们取的数量为2000，即：每个指定大小的样本随机产生2000次，然后计算每一次产生的样本的变量的相关系数。因此，蒙特卡罗模拟中的相关系数的样本分布主要有以下设计：四类大小的样本：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 &A=1 %THEN %DO; %LET SMPLSIZE=10; %END; %IF &A=2 %THEN %DO; %LET SMPLSIZE=20; %END; %IF &A=3 %THEN %DO; %LET SMPLSIZE=50; %END; %IF &A=4 %THEN %DO; %LET SMPLSIZE=100; %END;%DO B=1 %TO &NO_SMPL; *每类样本大小情况下，生成&NO_SMPL个随机样本;DATA DAT; *生成两个随机变量,其相关性为零; DO I=1 TO &SMPLSIZE; 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=&SMPLSIZE; IF _NAME_=X; CORR=Y; KEEP 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里已有分布的数据，例如标准正态分布（RANNOR）。II 将第I步生成的数据按需要的形状进行转换，例如将标准正态分布转换成其它正态分布III 已知变量相关性的多变量分布数据3.1生成一个SAS里已有分布的数据这里，我们需要生成两个相互独立的变量的数据集，因此，我们可以生成两个独立的正态分布随机变量，通过SAS的RANNOR函数可以实现随机正态变量的生成。在上一节的程序中，我们生成了两个变量，分别是X和Y，两者的均值为零，方差为1，即都符合标准正态分布N(0,1)。SAS数据生成我们将在以后的文章中更详细地介绍。3.2 将第I步生成的数据按需要的形状进行转换在很多情况下，我们需要将生成的符合某些标准分布的样本数据转换成其它分布的数据，此类转换主要有两类：一类是将样本数据转换成特定的均值和方差的数据，另一类是将数据转换成指定的形状（例如特定的峰度和偏度等）。对于第一类转换，我们可以直接通过公式进行线性转换即可，因为此类转换并不改变样本分布的正态性，公式为：X(new)=X*SD(new)+Mean(new)，这里X(new)即为转换后的变量，SD(new)为转换后变量的新的标准差，Mean(new)为转换后变量的新的均值。转换后的变量X(new)和原变量X有着同样的峰度和偏度。第二类转换则相对来说比较复杂，它需要一些统计模拟来生成特定峰度和偏度的非正态分布数据。此类转换将在以后的文章中更详细地介绍。3.3 已知变量相关性的多变量分布数据前面两种转换只针对单个变量，而实现情况中，我们需要生成有着相关性的多变量的数据此类变量的生成我们将在后面的文章中介绍。4 运用统计手段生成我们需要的统计量在第2节中，我们通过PROC CORR过程步计算每个样本的相关性，并将其保存在PEARSON数据集中，并最终汇总到数据集CORR.COR_RDM中。这样，我们就得到了4类样本大小各2000个相关系数的数据集，共4*2000=8000条数据。5 分析得到的统计量我们通过PROC MEANS来分析上一步得到的相关系数数据集，计算出每一类样本大小的相关系数的均值、标准差、最大值、最小值等统计量。程序如下：DATA A; SET CORR.COR_RDM; *生成一个临时数据集来存放CORR.COR_RDM，这样做的好处是保护好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 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.0 FONT=TRIPLEX) MINOR=NONE;PATTERN COLOR=BLACK VALUE=X2;PROC GCHART DATA=A; BY SMPLSIZE; *更多PROC GCHART的内容请参考前面的文章; VBAR CORR/ TYPE=PERCENT MIDPOINTS= -.9 TO .9 BY .05 RAXIS=AXIS1 MAXIS=AXIS2 WIDTH=1 SPACE=1;RUN; QUIT;结果如下：样本大小为10的相关系数结果：样本大小为20的相关系数结果：样本大小为50的相关系数结果：样本大小为100的相关系数结果：有关PROC GCHART的内容请参考前面的文章。6 基于蒙特卡罗模拟的结果得出结论我们从均值的结果可以看到，四类样本大小的相关系数均值均在0附近，说明它们的总体均值为零，但是四类样本大小的相关系数的方差却有很大的差别，样本越大，方差越小，反之则方差越大。从图形结果我们能更明显地看出这种趋势。样本越大，相关系数越集中在0附近，样本越小，相关系数分布就越分散。因此，我们在计算样本相关性的时候，其样本的大小一定要达到一定的数量，其结果才可靠，否则，如果样本大小太小，很容易产生一些有偏的结论。参考资料：Xitao Fan, etc.Monte Carlo Studies: 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 distribution：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, size, prob),quantile function 分位数函数。分位数：若概率0pZa)=的实数。如t分布的分位数表，自由度f=20和=0.05时的分位数为1.7247。 这个定义指的是上侧分位数。分位数：实数满足0 1 时，分位数是使PX x=F(x)=的数x；双侧分位数是使PX2=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掷骰子，掷到一即视为成功。则每次掷骰的成功率是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是希腊字母的读音)3.几何分布Geometric 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时，这是一个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 ()，泊松分布的参数是单位时间(或单位面积)内随机事件的平均发生率.泊松分布适合于描述单位时间内随机事件发生的次数。如某一服务设施在一定时间内到达的人数，电话交换机接到呼叫的次数，汽车站台的候客人数，机器出现的故障数，自然灾害发生的次数等等。rpois(n, lambda)，dpois(x,lambda)连续型6.均匀分布 Uniform Distribution，uniff(x) = 1/(max-min) for min = x = max。runif(n,min,max)。生成16位数的随机数：as.character(runif(1,1000000000000000,9999999999999999)。dunif(x,min,max)=1,恒定等于1/(max-min).对于连续变量，dfunction的值是x去特定值代入概率密度函数得到的函数值。7.正态分布Normal Distribution，normf(x) = 1/(sqrt(2 pi) sigma) e-(x mu)2/(2 sigma2)。其中mu是均值，sigma是standard deviation标准差。理论上可以证明如果把许多小作用加起来看做一个变量,那么这个变量服从正态分布rnorm(n,mean=0,sd=1)后两个参数如果不填则默认为0,1。dnorm(x,mean,sd),sd是标准差。画出正态分布概率密度函数的大致图形：x-seq(-3,3,0.1)，plot(x,dnorm(x) plot中的x,y要有相关关系才会形成函数图。qnorm(p,mean,sd),这个还是上侧分位数，如qnorm(0.05)=-1.644854,即x= 0, a 0 and s 0。Gamma分布中的参数，称为形状参数（shape parameter），即上式中的s，称为尺度参数（scale parameter）上式中的a。E(x)=s*a, Var(x)=s*a2. 当shape=1/2,scale=2时，这样的gamma分布是自由度为1的开方分布/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)pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = 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是分布的一个参数，常被称为率参数（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/(rexp10,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分钟以内来。按照以上分析一个小时出现的公交车次数应该不符合指数分布