数据分析团队策划书

此处看到的幻灯片标题2011.12.28蒙特卡罗方法在可靠性分析中的应用（一）何劼2014.3内容提要蒙特卡罗方法概述随机数的生成随机变量的抽样MC在PSA不确定性分析中的应用MC在HRA不确定性分析中的应用小结22121NNSS1S2S蒙特卡罗方法概述如何测量下面不规则图形的面积？3蒙特卡罗方法概述如何近似计算圆周率π？42121NN414/SS21NN4当N2=30000时，π的估计值与真实值只相差0.07%蒙特卡罗方法概述(续)蒙特卡罗(MonteCarlo,MC)方法又名随机模拟法或统计试验法。它是在第二次世界大战期间兴起和发展起来的。它的奠基人是冯.诺依曼(VonNeumann)。其主要思想是在计算机上模拟实际概率过程，然后加以统计处理。该方法与传统数学方法相比，具有直观性强、简便易行的优点，它能处理一些其它方法所不能的复杂问题，并且容易在计算机上实现。可以采用蒙特卡罗方法解决的问题分为两类，一类是问题本身具有概率或随机事件背景，可以直接用蒙特卡罗随机试验模拟；另一类是问题本身是确定性问题，需要通过构造随机事件来逼近真实情景。有记载的最早的蒙卡方法应用实例----蒲丰(Buffon)投针试验（1777）5优点1)能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。2)受几何条件限制小。3)收敛速度与问题的维数无关。4)具有同时计算多个方案与多个未知量的能力。5)误差容易确定。6)程序结构简单，易于实现。缺点1)收敛速度慢。2)误差具有概率性。蒙特卡罗方法概述(续)6蒙特卡罗方法的优点和缺点蒙特卡罗方法概述(续)PSA技术要素蒙特卡罗方法应用范围始发事件分析始发事件频率估计及不确定性分析，特殊始发事件建模及系统失效仿真系统分析非能动系统可靠性分析，系统失效概率计算仪控系统可靠性建模，Markov方法中方程组求解数据分析设备可靠性数据的参数估计，似然方程求解人员可靠性分析人员操作失误概率的不确定性传播定量化CDF/LRF的参数不确定性分析，PSA计算软件开发外部事件地震脆弱性分析火灾起火频率、水淹始发事件频率不确定性分析外部灾害发生频率估算（外部水淹、强风）PSA应用在役检查管道劣化机理分析，结构可靠性分析(SRA)，管道失效频率定量化蒙特卡罗方法在PSA技术中的应用7蒙特卡罗方法概述随机数的生成随机变量的抽样MC在PSA不确定性分析中的应用MC在HRA不确定性分析中的应用小结8随机数的生成在用蒙特卡罗方法进行系统可靠性分析和不确定性分析时，采用随机数代替设备可靠性参数作为PSA模型的输入，随机数质量的好坏直接影响定量化结果。通常所说的随机数是指［0，1］上均匀分布的连续型随机变量的子样。然而由于受到递推（迭代）公式和计算机字长的限制，在计算机上不可能生成“真正”的随机数，计算机生成的随机数只是近似［0，1］上均匀分布的随机数，称为（伪）随机数。常用的生成伪随机数的方法是线性同余法(LCG)，例如有乘同余法、乘加同余法、加同余法等。)(1nnT随机数与伪随机数(pseudorandomnumber)9随机数的生成(续)RiskSpectrum中不确定性分析设置界面10乘加同余法随机数的生成(续)•M通常取计算机中可以表示的最大的数，例如对于32位的计算机，M=232=4,294,967,296•a和c为（1，M-1）区间内的正整数，均为常数。•当i=0时，初始值x0称为种子，一般可任取一个正整数。•例如，当M=16，a=5，c=13，种子x0=12345时，计算过程如下：表示axi+c除以M之后的余数caxxii1，Mxii11)(modM,2,1,0i11随机数的生成(续)ixiξi012345—1(5*12345+13)/16=3858余1010/16=0.6252(5*10+13)/16=3余1515/16=0.93753(5*15+13)/16=5余88/16=0.54(5*8+13)/16=3余55/16=0.31255(5*5+13)/16=2余66/16=0.375M=16，a=5，c=13，种子x0=123451351iixx，1611iix)16(mod,2,1,0i得到随机数序列为：0.625，0.9375，0.5，0.3125，0.375，…12随机数的生成(续)M=16，a=5，c=13，种子x0=54321ixiξi054321—1(5*54321+13)/16=16976余22/16=0.1252(5*2+13)/16=1余77/16=0.43753(5*7+13)/16=3余00/16=04(5*0+13)/16=0余1313/16=0.81255(5*13+13)/16=4余1414/16=0.875得到随机数序列为：0.125，0.4375，0，0.8125，0.875，…种子唯一确定了随机数序列，当种子的取值给定之后整个随机数序列就不再“随机”了；在进行蒙特卡罗模拟时，不应当刻意选定某个数值作为种子，否则就违背了“随机试验”的理念。任何蒙特卡罗试验的结论都不应当对种子的取值具有依赖性。13随机数的生成(续)生成伪随机数时需要考虑的问题：•伪随机数的质量（独立性、均匀性）；•伪随机数序列的周期和容量；•方法的经济性。针对上面的问题，生成伪随机数的策略是：•选取合适的递推公式T，设置合适的参数值。例如对于乘加同余法，M的所有质数因子P，应当满足a≡1(modP)。当M=2s(s为正整数)时，取a≡1(mod4)，c≡1(mod2)时乘加同余法生成的伪随机数质量最好；•用蒙特卡罗方法进行有限次试验时，尽量避免试验次数不要超过伪随机数序列的周期；•如果调用系统自带的随机数生成函数，如果不知道其后台算法，可以先随机生成种子再根据随机的种子生成随机数，避免在同一个随机数序列中多次取随机数。14随机数的生成(续)ixiξi+1012345—1(7*12345+13)/16=5401余1212/16=0.752(7*12+13)/16=6余11/16=0.06253(7*1+13)/16=1余44/16=0.254(7*4+13)/16=2余99/16=0.56255(7*9+13)/16=4余1210/16=0.75M=16，a=7，c=13，种子x0=12345循环13x7xi1i，1611iix)16(mod,2,1,0i15蒙特卡罗方法概述随机数的生成随机变量的抽样MC在PSA不确定性分析中的应用MC在HRA不确定性分析中的应用小结16随机变量的抽样由已知分布的随机抽样指的是由己知分布的随机变量总体中抽取简单子样。随机数序列是由单位均匀分布的总体中抽取的简单子样，属于一种特殊的由已知分布的随机抽样问题。为方便起见，用ξn表示[0,1]区间上均匀分布的随机数，用XF表示由己知分布F(x)中产生的简单子样的个体。对于连续型分布，常用分布密度函数f(x)表示总体的己知分布，用Xf表示由己知分布密度函数f(x)产生的简单子样的个体。常用的抽样方法有直接抽样法、舍选法等。目前在PSA中常用的Gamma、Beta和Lognormal分布都已有较成熟的抽样方法。17随机变量的抽样(续)直接抽样法对于任意给定的分布函数F(x)，直接抽样方法如下：其中，ξ1，ξ2，…，ξn为随机数序列。为方便起见，将上式简化为：直接抽样法的优点是抽样效率高、费用低，缺点是需要F(x)和F-1(x)都能写出解析的表达式。因此直接抽样法主要适用于离散型随机变量抽样。N,,2,1n,tinfXn)t(FntinfX)t(FF使得F(t)≥ξ的所有t值当中最小的18随机变量的抽样(续)直接抽样法（图示）190.000.100.200.300.400.500.600.700.800.901.001.00E+001.00E+011.00E+021.00E+031.00E+041.00E+051.00E+06ξ1ξ2ξnX1X2X3F(t)随机变量的抽样(续)舍选法（Acceptance-RejectionMethod,AR）舍选法是VonNeumann提出的，是应用最广泛也是研究成果最多的抽样方法。设密度函数f(x)只在有限区间上为正，且最大值有限，即当a≤x≤b时独立地生成两个均匀随机数ξ，η~U(0,1)，令Y=a+(b-a)ξ如果Mη≤f(ξ)则取XF=ξ，否则重新生成ξ和η。M)x(f020随机变量的抽样(续)舍选法的几何意义f(x)21ab随机变量的抽样(续)PSA常用分布的随机抽样分布类抽样方案正态分布U，V为[0,1]区间上的随机数对数正态分布将正态分布随机数取指数Gamma分布根据a1和0a1的情况分别采用舍选法，可得到Г(a,1)，根据尺度特性可得到Г(a,b)Beta分布由Gamma分布得到，若X1~Г(a,1)，X2~Г(b,1)，且二者独立，则X1/(X1+X2)~Beta(a,b))U2cos(Vln2X)U2sin(Vln2Y22随机变量的抽样(续)蒙特卡罗方法的收敛性•蒙特卡罗方法是由随机变量X的简单子样X1，X2，…，XN的算术平均值：作为所求解的近似值。由大数定律可知，如X1，X2，…，XN独立同分布，且具有有限期望值（E(X)∞），则即随机变量X的简单子样的算术平均值，当子样数Ｎ充分大时，以概率1收敛于它的期望值E(X)。NiiNXNX111)(limXEXPNN23随机变量的抽样(续)蒙特卡罗方法的误差•根据中心极限定理，有如下近似式：•蒙特卡罗方法的误差ε定义为上式中λα表示置信度为α的标准正态分布的分位点，σ表示随机变量X的标准差，N为样本容量。误差收敛速度的阶为。显然，当给定置信度α后，误差ε由σ和N决定。要想减小ε，要么增大N，要么减小σ。在σ固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。所以，在做PSA不确定性分析时，如果将模拟次数由5000增加一倍到10000，误差并不会随之减小一倍。N)(2/1NO122)(02/2dteNXEXPtN24蒙特卡罗方法概述随机数的生成随机变量的抽样MC在PSA不确定性分析中的应用MC在HRA不确定性分析中的应用小结25MC在PSA不确定性分析中的应用不确定性分析算法流程：26MC在PSA不确定性分析中的应用(续)失效概率分布类型EFA3E-4LogN5B1E-3LogN10例1.基本事件A、B的可靠性参数如下表，如何求顶事件不可用度的不确定性区间和概率密度函数（PDF）曲线？27MC在PSA不确定性分析中的应用(续)根据对数正态分布的性质，将基本事件A和B的失效概率转换为分别服从LogN(-8.59，0.98)和LogN(-7.89，1.40)分布。生成5000个分别服从LogN(-8.59，0.98)和LogN(-7.89，1.40)的随机数。由于顶门是与门，因此将A和B两个基本事件的5000个抽样数据分别相乘，得到顶事件不可用度的5000个样点值。28MC在PSA不确定性分析中的应用(续)将顶门的5000个样点值从小到大排序，分别取出第50、100、150、…、5000个点，共100个点值，对应于1%、2%、3%、…、100%分位点。以这100个分位点值作为横坐标，以1%、2%、3%、…、100%作为纵坐标，做出的曲线就是顶门的不可用度的累积密度函数（cdf，分布函数）图像。29MC在PSA不确定性分析中的应用(续)0.000.200.400.600.801.001.201.00E-101.00E-091.00E-081.00E-071.00E-061.00E-051.00E-04顶事件不可用度的累计密度函数图像30MC在PSA不确定性分析中的应用(续)由于概率