统计模拟介绍..

专业必修课教学目的：掌握统计模拟在实际中的应用36理论课+18实验课(PBL)考试形式：闭卷考试(60%)+小论文(40%)1.统计模拟，Ross著，王兆军，陈广雷，邹长亮译2007年7月由人民邮电出版社出版2.统计建模与R软件，薛毅,陈广萍编著，2007年4月清华大学出版社MonteCarlo方法：亦称统计模拟方法，statisticalsimulationmethod利用随机数进行数值模拟的方法MonteCarlo名字的由来：•是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；•MonteCarlo是摩纳哥（monaco)的首都，该城以赌博闻名NicholasMetropolis(1915-1999)Monte-Carlo,MonacoMonteCarlo方法简史简单地介绍一下MonteCarlo方法的发展历史1、Buffon投针实验：1768年，法国数学家ComtedeBuffon利用投针实验估计的值dLp2dLProblemofBuffon’sneedle:Ifaneedleoflengthlisdroppedatrandomonthemiddleofahorizontalsurfaceruledwithparallellinesadistancedlapart,whatistheprobabilitythattheneedlewillcrossoneofthelines?Solution:Thepositioningoftheneedlerelativetonearbylinescanbedescribedwitharandomvectorwhichhascomponents:),0[),0[dATherandomvectorisuniformlydistributedontheregion[0,d)×[0,).Accordingly,ithasprobabilitydensityfunction1/d.TheprobabilitythattheneedlewillcrossoneofthelinesisgivenbytheintegraldldAdpld20sin011777年，古稀之年的蒲丰在家中请来好些客人玩投针游戏（针长是线距之半），他事先没有给客人讲与π有关的事。客人们虽然不知道主人的用意，但是都参加了游戏。他们共投针2212次，其中704次相交。蒲丰说，2212/704=3.142，这就是π值。这着实让人们惊喜不已。2、1930年，EnricoFermi利用MonteCarlo方法研究中子的扩散，并设计了一个MonteCarlo机械装置，Fermiac,用于计算核反应堆的临界状态3、VonNeumann是MonteCarlo方法的正式奠基者,他与StanislawUlam合作建立了概率密度函数、反累积分布函数的数学基础，以及伪随机数产生器。在这些工作中，StanislawUlam意识到了数字计算机的重要性合作起源于Manhattan工程：利用ENIAC(ElectronicNumericalIntegratorandComputer)计算产额一些人进行了实验，其结果列于下表：实验者年份投计次数π的实验值沃尔弗(Wolf)185050003.1596斯密思(Smith)185532043.1553福克斯(Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929蒙特卡罗方法又称计算机随机模拟方法。它是以概率统计理论为基础的一种方法。由蒲丰试验可以看出，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。这就是蒙特卡罗方法的基本思想。l=1;d=2;m=0;n=10000fork=1:n;x=unifrnd(0,d);y=unifrnd(0,pi);ifx1*sin(y)m=m+1elseendendp=m/npi_m=1/p关系式成立产生随机数验证模型成立次数k=k+1否是计算估计结果k/n成立次数不变试验次数是否达到n次是否编写R程序MonteCarlo模拟的应用：自然现象的模拟：宇宙射线在地球大气中的传输过程；高能物理实验中的核相互作用过程；实验探测器的模拟数值分析：利用MonteCarlo方法求积分金融工程：股票期权的模拟定价物理化学生物环境工程医学金融交通教育心里卫生数学语言军事历史经济天文。。。MonteCarlo模拟在物理研究中的作用MonteCarlo模拟的步骤：1.根据欲研究的物理系统的性质，建立能够描述该系统特性的理论模型，导出该模型的某些特征量的概率密度函数；2.从概率密度函数出发进行随机抽样，得到特征量的一些模拟结果；3.对模拟结果进行分析总结，预言物理系统的某些特性。注意以下两点：•MonteCarlo方法与数值解法的不同:MonteCarlo方法利用随机抽样的方法来求解物理问题;数值解法:从一个物理系统的数学模型出发,通过求解一系列的微分方程来的导出系统的未知状态;•MonteCarlo方法并非只能用来解决包含随机的过程的问题:许多利用MonteCarlo方法进行求解的问题中并不包含随机过程例如:用MonteCarlo方法计算定积分.对这样的问题可将其转换成相关的随机过程,然后用MonteCarlo方法进行求解MonteCarlo算法的主要组成部分概率密度函数(pdf)—必须给出描述一个物理系统的一组概率密度函数;随机数产生器—能够产生在区间[0,1]上均匀分布的随机数抽样规则—如何从在区间[0,1]上均匀分布的随机数出发,随机抽取服从给定的pdf的随机变量;模拟结果记录—记录一些感兴趣的量的模拟结果误差估计—必须确定统计误差（或方差）随模拟次数以及其它一些量的变化；减少方差的技术—利用该技术可减少模拟过程中计算的次数；并行和矢量化—可以在先进的并行计算机上运行的有效算法1.具有统计功能的软件Excel,Matab,C,Fortran2.专业的统计软件SPSS,SAS,S-Plus,R,Gauss,minitab优点1)能够比较逼真地描述具有随机性质的事物的特点及物理实验过程。2)受几何条件限制小。3)收敛速度与问题的维数无关。4)具有同时计算多个方案与多个未知量的能力。5)误差容易确定。6)程序结构简单，易于实现。缺点1)收敛速度慢。2)误差具有概率性。3)在粒子输运问题中，计算结果与系统大小有关。火车离站时刻13:0013:0513:10概率0.70.20.1一列列车从A站开往B站，某人每天赶往B站上车。他已经了解到火车从A站到B站的运行时间是服从均值为30min，标准差为2min的正态随机变量。火车大约下午13：00离开A站，此人大约13:30到达B站。火车离开A站的时刻及概率如表1所示，此人到达B站的时刻及概率如表2所示。问此人能赶上火车的概率有多大？表1：火车离开A站的时刻及概率表2：某人到达B站的时刻及概率人到站时刻13:2813:3013:3213:34概率0.30.40.20.1——问题的分析——这个问题用概率论的方法求解十分困难，它涉及此人到达时刻、火车离开站的时刻、火车运行时间几个随机变量，而且火车运行时间是服从正态分布的随机变量，没有有效的解析方法来进行概率计算。在这种情况下可以用计算机模拟的方法来解决。：火车从A站出发的时刻；：火车从A站到B站的运行时间；：某人到达B站的时刻；：随机变量服从正态分布的均值；：随机变量服从正态分布的标准差；进行计算机统计模拟的基础是抽象现实系统的数学模型为了便于建模，对模型中使用的变量作出如下假定：1T2T3T2T2T此人能及时赶上火车的充分必要条件为：，所以此人能赶上火车的概率模型为：。123TTT123{}pTTT为了分析简化，假定13时为时刻t=0，则变量、的分布律为：1T3T05100.70.20.1283032340.30.40.20.11/minT()Pt3/minT()PtR软件求解的总算法：关系式成立产生随机数验证模型成立次数k=k+1否是计算估计结果k/n成立次数不变试验次数是否达到n次是否编写R程序①借助区间(0,1)分布产生的随机数，对变量、概率分布进行统计模拟；1T3T②根据变量、、概率分布及模拟程序、命令产生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(28,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=0while(i=r){i=i+1segments(x,1,x+i,1)if(x+i=y)points(y,1,pch=19)Sys.sleep(0.1)}points(y,1,pch=19)title(ifelse(x+r=y,poor...missedthetrain!,Bingo!catchedthetrain!))Sys.sleep(0.5)x+ry})mean(prb)1.1矩母函数和生成函数定义1.1设随机变量X的密度函数为p(x),称为X的矩母函数，记为性质：(1)(2)若X和Y相互独立，则Xgt0nnEXgXYXYgtgtgtexpEtX定义1.2若Ｘ为离散型随机变量，称为其概率生成函数，记为性质(1)(2)若X和Y独立，XEss2'1,''1'1EXEXXYXYsss1.条件分布其中为X的边际分布例子：令X和Y的联合密度函数为试求,|XpxypyxpxXpx,01,01,0xyxypxyelse11|43PXY2.条件数学期望命题：,||,yfxydyEYXxypyxdyfxydy|EEYXEY3.条件方差条件方差公式222|||||VarYXEYEYXXEYXEYX||VarYEVarYXVarEYX例1.3从某大学任意挑选一个学院，然后从此学院中任意挑选n个学生，令X表示这些学生中来自武汉市的人数，令Q代表该学院来自武汉市的人数所占的比例，因为学院之间的比例不相同，因此Q也是一个随机变量。若Q~U(0,1),X|Q=q~B(n,q)，求X的方差。随机过程是一个随机变量集合，状态空间S是随机变量取值集合，集合T为指标集。指标集可以是离散的也可以是连续的。例：天气变化情况是一个随机过程晴，晴，多云，雨，雨，…,tXtTtXPoisson计数过程定义2-4：计数过程｛N(t)，t≥0｝被称为强度参数为λ的Poisson过程，若：①独立增量，平稳；②对任意h＞0，有P｛N(h)=1｝=λh+ο(h)；③对任意h＞0，有P｛N(h)≥2｝=ο(h)。命题：设为Poisson过程中事件发生的间隔时间序列，则为独立同分布的随机变量序列，且共同分布为具有参数的指数分布1234,,,...XXXX1234,,,...XXXX例1.5顾客依