统计模拟介绍

1、专业必修课教学目的：掌握统计模拟在实际中的应用36理论课+18实验课(PBL)考试形式：闭卷考试(60%)+小论文(40%)1.统计模拟，Ross著，王兆军，陈广雷，邹长亮译 2007年7月由人民邮电出版社出版2.统计建模与R软件，薛毅,陈广萍编著，2007年4月清华大学出版社Monte CarloMonte Carlo方法：方法：亦称统计模拟方法，亦称统计模拟方法，statistical simulation methodstatistical simulation method 利用随机数进行数值模拟的方法利用随机数进行数值模拟的方法Monte CarloMonte Carlo名字的由来：

2、名字的由来： 是由是由MetropolisMetropolis在二次世界大战期间提出的：在二次世界大战期间提出的：ManhattanManhattan计划，研究与原子计划，研究与原子弹有关的中子输运过程；弹有关的中子输运过程； Monte CarloMonte Carlo是摩纳哥（是摩纳哥（monaco)monaco)的首都，该城以赌博闻名的首都，该城以赌博闻名Nicholas Metropolis (1915-1999)Monte-Carlo, MonacoMonte CarloMonte Carlo方法简史方法简史简单地介绍一下简单地介绍一下Monte CarloMonte Carlo方法

3、的发展历史方法的发展历史1 1、BuffonBuffon投针实验：投针实验：17681768年，法国数学家年，法国数学家Comte de Buffon利用投针实验估计利用投针实验估计 的的值值dLp2dLProblem of Buffons needle:Problem of Buffons needle:If a needle of lengthIf a needle of length l l is dropped at random on the is dropped at random on the middle of a horizontal surface ruled with p

4、arallel middle of a horizontal surface ruled with parallel lines a distance lines a distance d d l l apart, what is the probability apart, what is the probability that the needle will cross one of the lines?that the needle will cross one of the lines?SolutionSolution:The positioning of the needle Th

5、e positioning of the needle relative to nearby lines can be relative to nearby lines can be described with a random vector described with a random vector which has components:which has components:),0),0dAThe random vector is uniformly distributed on the region The random vector is uniformly distribu

6、ted on the region 0,0,d d) )0,0, ). Accordingly, it has probability density function 1/d). Accordingly, it has probability density function 1/d . .The probability that the needle will cross one of the lines is given The probability that the needle will cross one of the lines is given by the integral

7、by the integraldldAdpld20sin01 1777年，古稀之年的蒲丰在家中请来好些客人玩投针游戏（针长是线距之半），他事先没有给客人讲与有关的事。客人们虽然不知道主人的用意，但是都参加了游戏。他们共投针2212次，其中704次相交。蒲丰说，2212/704=3.142，这就是值。这着实让人们惊喜不已。2 2、19301930年，年，Enrico FermiEnrico Fermi利用利用Monte CarloMonte Carlo方法研究中子方法研究中子的扩散，并设计了一个的扩散，并设计了一个Monte CarloMonte Carlo机械装置，机械装置，Fermiac,F

8、ermiac,用于计算核反应堆的临界状态用于计算核反应堆的临界状态3 3、Von NeumannVon Neumann是是Monte CarloMonte Carlo方法的正式奠基者方法的正式奠基者, ,他与他与Stanislaw UlamStanislaw Ulam合作建立了概率密度函数、反累积分布函数合作建立了概率密度函数、反累积分布函数的数学基础，以及伪随机数产生器。在这些工作中，的数学基础，以及伪随机数产生器。在这些工作中， Stanislaw UlamStanislaw Ulam意识到了数字计算机的重要性意识到了数字计算机的重要性合作起源于合作起源于ManhattanManhatta

9、n工程：利用工程：利用ENIAC(Electronic ENIAC(Electronic Numerical Integrator and Computer)Numerical Integrator and Computer)计算产额计算产额一些人进行了实验，其结果列于下表 ：实验者年份投计次数的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929蒙特卡罗方法又称计算机随机模拟方法。它是以概率统计理论为基础的一种方法。 由蒲丰试验可以看出，

10、当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。 l=1;d=2;m=0;n=10000for k=1:n;x=unifrnd(0,d);y=unifrnd(0,pi);if x1*sin(y)m=m+1elseendendp=m/npi_m=1/p关系式成立产生随机数验证模型成立次数k=k+1否是计算估计结果k/n成立次数不变试验次数是否达到n次是否编写R程序Monte CarloMonte Carlo模拟的应

11、用：模拟的应用：自然现象的模拟：自然现象的模拟：宇宙射线在地球大气中的传输过程；宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；高能物理实验中的核相互作用过程；实验探测器的模拟实验探测器的模拟数值分析：数值分析：利用利用Monte CarloMonte Carlo方法求积分方法求积分金融工程：金融工程：股票期权的模拟定价股票期权的模拟定价物理 化学 生物 环境工程 医学 金融 交通教育 心里 卫生 数学语言 军事 历史 经济天文。Monte CarloMonte Carlo模拟在物理研究中的作用模拟在物理研究中的作用Monte CarloMonte Carlo模拟的步骤：模拟的

12、步骤：根据欲研究的物理系统的性质，建立能够描述该系统特性根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；的理论模型，导出该模型的某些特征量的概率密度函数；从概率密度函数出发进行随机抽样，得到特征量的一些模从概率密度函数出发进行随机抽样，得到特征量的一些模拟结果；拟结果；对模拟结果进行分析总结，预言物理系统的某些特性。对模拟结果进行分析总结，预言物理系统的某些特性。注意以下两点：注意以下两点：Monte CarloMonte Carlo方法与数值解法的不同方法与数值解法的不同: :Monte CarloMonte Carlo方法利用随机抽样的方

13、法来求解物理问方法利用随机抽样的方法来求解物理问题题; ;数值解法数值解法: :从一个物理系统的数学模型出发从一个物理系统的数学模型出发, ,通过求解一通过求解一系列的微分方程来的导出系统的未知状态系列的微分方程来的导出系统的未知状态; ;Monte CarloMonte Carlo方法并非只能用来解决包含随机的过程的问题方法并非只能用来解决包含随机的过程的问题: :许多利用许多利用Monte CarloMonte Carlo方法进行求解的问题中并不包含方法进行求解的问题中并不包含随机过程随机过程 例如例如: :用用Monte CarloMonte Carlo方法计算定积分方法计算定积分. .

14、 对这样的问题可将其转换成相关的随机过程对这样的问题可将其转换成相关的随机过程, , 然后用然后用Monte CarloMonte Carlo方法进行求解方法进行求解Monte CarloMonte Carlo算法的主要组成部分算法的主要组成部分概率密度函数概率密度函数(pdf)(pdf) 必须给出描述一个物理系统的一组概必须给出描述一个物理系统的一组概率密度函数率密度函数; ;随机数产生器随机数产生器能够产生在区间能够产生在区间0,10,1上均匀分布的随机数上均匀分布的随机数抽样规则抽样规则如何从在区间如何从在区间0,10,1上均匀分布的随机数出发上均匀分布的随机数出发, ,随随机抽取服从给

15、定的机抽取服从给定的pdfpdf的随机变量的随机变量; ;模拟结果记录模拟结果记录记录一些感兴趣的量的模拟结果记录一些感兴趣的量的模拟结果误差估计误差估计必须确定统计误差（或方差）随模拟次数以及必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；其它一些量的变化；减少方差的技术减少方差的技术利用该技术可减少模拟过程中计算的次数；利用该技术可减少模拟过程中计算的次数；并行和矢量化并行和矢量化可以在先进的并行计算机上运行的有效算法可以在先进的并行计算机上运行的有效算法1.具有统计功能的软件 Excel, Matab, C, Fortran 2. 专业的统计软件 SPSS, SAS, S-Pl

16、us, R, Gauss, minitab优点 能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。 受几何条件限制小。 收敛速度与问题的维数无关。 具有同时计算多个方案与多个未知量的能力。 误差容易确定。 程序结构简单，易于实现。 缺点 收敛速度慢。 误差具有概率性。 在粒子输运问题中，计算结果与系统大小有关。火车离站时刻火车离站时刻13:0013:0513:10概率概率0.70.20.1 一列列车从一列列车从A A站开往站开往B B站，某人每天赶往站，某人每天赶往B B站上车。他站上车。他已经了解到火车从已经了解到火车从A A站到站到B B站的运行时间是服从均值为站的运行时间是服从均

17、值为30min30min，标准差为，标准差为2min2min的正态随机变量。火车大约下午的正态随机变量。火车大约下午1313：0000离开离开A A站，此人大约站，此人大约13:3013:30到达到达B B站。火车离开站。火车离开A A站的时刻站的时刻及概率如表及概率如表1 1所示，此人到达所示，此人到达B B站的时刻及概率如表站的时刻及概率如表2 2所示。所示。问此人能赶上火车的概率有多大？问此人能赶上火车的概率有多大？表1：火车离开A站的时刻及概率 表2：某人到达B站的时刻及概率 人到站时刻人到站时刻13:2813:3013:3213:34概率概率0.30.40.20.1这个问题用概率论的

18、方法求解十分困难，它涉及此人到达时刻、火车离开站的时刻、火车运行时间几个随机变量，而且火车运行时间是服从正态分布的随机变量，没有有效的解析方法来进行概率计算。在这种情况下可以用计算机模拟的方法来解决。：火车从：火车从A A站出发的时刻；站出发的时刻；：火车从：火车从A A站到站到B B站的运行时间；站的运行时间；：某人到达：某人到达B B站的时刻；站的时刻；：随机变量：随机变量 服从正态分布的均值；服从正态分布的均值；：随机变量：随机变量 服从正态分布的标准差服从正态分布的标准差；1T2T3T2T2T此人能及时赶上火车的充分必要条件为：此人能及时赶上火车的充分必要条件为： ，所以此人能赶上火车

19、的概率模型为：所以此人能赶上火车的概率模型为： 。123TTT123p TTT为了分析简化，假定为了分析简化，假定1313时为时刻时为时刻t=0，则变量，则变量 、 的分布律为：的分布律为：1T3T05100.70.20.1283032340.30.40.20.11/minT( )P t3/minT( )P t关系式成立产生随机数验证模型成立次数k=k+1否是计算估计结果k/n成立次数不变试验次数是否达到n次是否编写R程序借助区间借助区间(0,1)(0,1)分布产生的随机数，对分布产生的随机数，对变量变量 、 概率分布进行统计模拟；概率分布进行统计模拟；1T3T根据变量根据变量 、 、 概率分

20、布及模拟概率分布及模拟程序、命令产生程序、命令产生n 个随机分布数；个随机分布数；1T2T3T使用随机产生的使用随机产生的n 组随机数验证模型中组随机数验证模型中的关系表达式是否成立；的关系表达式是否成立；计算计算n 次模拟实验中，使得关系表达次模拟实验中，使得关系表达式成立的次数式成立的次数k ；当当 时，以时，以 作为此人能赶作为此人能赶上火车的概率上火车的概率p 的近似估计；的近似估计；nknwindows(7, 3)prb = replicate(100, x = sample(c(0, 5, 10), 1, prob = c(0.7, 0.2, 0.1) y = sample(c(2

21、8, 30, 32, 34), 1, prob = c(0.3, 0.4, 0.2, 0.1) plot(0:40, rep(1, 41), type = n, xlab = time, ylab = ,axes = FALSE) axis(1, 0:40) r = rnorm(1, 30, 2) points(x, 1, pch = 15) i = 0 while (i = y) points(y, 1, pch = 19) Sys.sleep(0.1) points(y, 1, pch = 19) title(ifelse(x + r y )mean(prb)1.1 矩母函数和生成函数 定

22、义1.1 设随机变量X的密度函数为p(x),称 为X的矩母函数，记为性质：(1) (2) 若X和Y相互独立，则 Xgt 0nnE Xg X YXYgtgt gtexpEtX定义1.2若为离散型随机变量，称为其概率生成函数，记为性质(1)(2)若X和Y独立， XE s s 2 1 , 1 1E XE X X YXYsss 1.条件分布其中 为X的边际分布例子：令X和Y的联合密度函数为试求 ,|Xp x yp y xpx Xpx,01, 01,0 xyxyp x yelse11|43P XY2.条件数学期望命题： ,|,yfx y dyE YXxyp y x dyfx y dy |E E YXE

23、Y 3.条件方差条件方差公式222|Var YXEYE YXXE YXE YX |Var YE Var YXVar E YX例1.3 从某大学任意挑选一个学院，然后从此学院中任意挑选n个学生，令X表示这些学生中来自武汉市的人数，令Q 代表该学院来自武汉市的人数所占的比例，因为学院之间的比例不相同，因此Q也是一个随机变量。若QU(0,1), X|Q=qB(n,q)，求X的方差。随机过程 是一个随机变量集合，状态空间S 是随机变量 取值集合，集合T 为指标集。指标集可以是离散的也可以是连续的。例：天气变化情况是一个随机过程 晴，晴，多云，雨，雨， ,tX tTtX命题：设 为Poisson过程中事

24、件发生的间隔时间序列，则 为独立同分布的随机变量序列，且共同分布为具有参数 的指数分布1234,.XXXX1234,.XXXX 例1.5 顾客依Poisson过程到达商店，速率为 人/小时。已知商店上午9:00开门。试求：到9:30时仅到一位顾客，而到11:30时总计已到5位顾客的概率。.4 2.非齐次Posson过程一、Markov链的定义1Markov链,) 1(,)(|) 1(1ninXinXjnXP) 0 (,) 1 (,01iXiX)(|) 1(inXjnXP有限马氏链状态空间是有限集I=0,1,2,，k2一步转移概率马氏链在时刻n处于状态 i 的条件下，到时刻n+1转移到状态 j

25、的条件概率，即|1iXjXPnn称为在时刻n的一步转移概率，记 作)(npij注:马氏链由 和条件概率 决定00P Xi11|nnnnP Xi Xi注:由于概率是非负的，且过程从一状态出发，经过一步转移后，必到达状态空间中的某个状态一步转移概率满足3一步转移矩阵称为在时刻n的一步转移矩阵 随机 矩阵即有有限马氏链状态空间I=0，1，2，k)()()()()()(10111001001npnpnpnpnpnpPnn)()()()()()()()()(1011110001001npnpnpnpnpnpnpnpnpPkkkkkk4齐次马氏链即则称此马氏链为齐次马氏链（即关于时间为齐次）如果马氏链的一步转移概率)(npij与 n 无关，ijnnpiXjXP|15初始分布设)(00iXPip，Ii，如果对一切Ii都有0)(0ip1)(0ipIi称)(0ip为马氏链的初始分布注马氏链在初始时刻有可能处于I中任意状态，初始分布就是马氏链在初始时刻的概率分布。6绝对分布概率分布)(iXPipnn，Ii，0n称为马氏链的绝对分布或称绝对概率13452. 随例1 （机游动）12 3 4 55),Iii 设一醉汉在， ， ， ，作随机游动：如果现在位于点（1则下一时刻各以1/3概率向左或向右移动一