基于R软件的统计模拟(精品)

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,基于,R,软件的统计模拟,奚 潭,（南京财经大学统计系,2006,级）,主要内容,1.,统计模拟的基本概念,2.,赶火车问题,3.R,软件的统计模拟功能,4.,应用,R,软件模拟验证大数定律,5.,应用,R,软件模拟验证中心极限定理,一、统计模拟的基本概念,（一）统计模拟的定义,统计模拟即是计算机统计模拟，它实质上是计算机建模，而这里的计算机模型就是计算机方法、统计模型,(,如程序、流程图、算法等,),，它是架于计算机理论和实际问题之间的桥梁。它与统计建模的关系如下图。,实际问题,统计、逻辑 模型,计算机模拟（程序、算法）,统计、计算机解,实际解,一、统计模拟的基本概念,（二）统计模拟方法,一般地，统计模拟分类如下：,若按状态变量的变化性质分为,连续随机模拟,和,离散随机模拟,。,而按变量是否随时间变化又可分为,动态随机模拟,和,静态随机模拟,。,常用的统计模拟方法主要有以下几种：,1.,蒙特卡罗法,2.,系统模拟方法,3.,其它方法：,包括,Bootstrap(,自助法,),、,MCMC,（马氏链蒙特卡罗法）等。,一、统计模拟的基本概念,（三）统计模拟的一般步骤,二、赶火车问题,火车离站时刻,13:00,13:05,13:10,概率,0.7,0.2,0.1,一列列车从,A,站开往,B,站，某人每天赶往,B,站上车。他已经了解到火车从,A,站到,B,站的运行时间是服从均值为,30min,，标准差为,2min,的正态随机变量。火车大约下午,13,：,00,离开,A,站，此人大约,13:30,到达,B,站。火车离开,A,站的时刻及概率如表,1,所示，此人到达,B,站的时刻及概率如表,2,所示。问此人能赶上火车的概率有多大？,表,1,：火车离开,A,站的时刻及概率,表,2,：某人到达,B,站的时刻及概率,人到站时刻,13:28,13:30,13:32,13:34,概率,0.3,0.4,0.2,0.1,二、赶火车问题,问题的分析,这个问题用概率论的方法求解十分困难，它涉及此人到达时刻、火车离开站的时刻、火车运行时间几个随机变量，而且火车运行时间是服从正态分布的随机变量，没有有效的解析方法来进行概率计算。在这种情况下可以用计算机模拟的方法来解决。,：火车从,A,站出发的时刻；,：火车从,A,站到,B,站的运行时间；,：某人到达,B,站的时刻；,：随机变量 服从正态分布的均值；,：随机变量 服从正态分布的标准差；,二、赶火车问题,进行计算机统计模拟的基础是抽象现实系统的数学模型,为了便于建模，对模型中使用的变量作出如下假定：,此人能及时赶上火车的充分必要条件为：，所以此人能赶上火车的概率模型为：。,二、赶火车问题,为了分析简化，假定,13,时为时刻,t=0,，则变量 、的分布律为：,0,5,10,0.7,0.2,0.1,28,30,32,34,0.3,0.4,0.2,0.1,二、赶火车问题,R,软件求解的总算法：,关系式,成立,产生随机数,验证模型,成立次数,k=k+1,否,是,计算估计结果,k/n,成立次数不变,试验次数,是否达到,n,次,是,否,编写,R,程序,借助区间,(0,1),分布产生的随机数，对变量 、概率分布进行统计模拟；,根据变量 、概率分布及模拟程序、命令产生,n,个随机分布数；,使用随机产生的,n,组随机数验证模型中的关系表达式是否成立；,计算,n,次模拟实验中，使得关系表达式成立的次数,k,；,当 时，以 作为此人能赶上火车的概率,p,的近似估计；,进入演示,windows(7,3),prb,=replicate(100,#,括号内程序重复,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=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,),进入模拟,三、,R,软件,的统计模拟功能,1,、,R,软件优秀的随机数模拟功能,分布,产生随机数序列命令,参数设置,binomial,rbinom,(),n,size,prob,chi-squared,rchisq,(),n,df,ncp,exponential,exp(),n,rate,F,F(),n,df1,df2,ncp,normal,norm(),n,mean,sd,Poisson,pois,(),n,lambda,Students t,t(),n,df,ncp,unifom,unif,(),n,min,max,生产某概率分布的随机数是实现统计模拟的前提条件，而使用,R,命令可以生成以下常用分布的随机数：,三、,R,软件,的,统计模拟功能,2,、优良的编程环境和编程语言,R,所拥有的好的兼容性、拓展性和强大的内置函数有利于统计模拟的实现。,3,、高效率的向量运算功能,使用,R,拥有的向量运算功能可以大大减少程序运行的时间，提高程序运行的效率,。,下面以求解,Pi,的程序为例加以说明,未采用,R,向量运算功能的程序为：,mc1-function(n),set.seed(1234579),k-0;,x-,runif(n,);,y-,runif(n,);,for(i in 1:n),if(xi2+yi21),k-k+1;,data.frame(Pi,=4*k/n),引入向量运算功能改进后的程序为：,mc1-function(n),set.seed(1234579),k-0;,x-,runif(n,);,y-,runif(n,);,k-length(xx2+y2,下面用,R,软件,分别执行两个程序，看看有什么差异,程序,1,.,程序,2,三、,R,软件,的统计模拟功能,四、应用,R,软件,模拟验证大数定律,1,、验证的大数定律有：,（,1,）伯努利大数定理,设 是 次独立重复试验中事件 发生的次数。是事件 在每次试验中发生的概率，则对于任意正数,0,有,（,2,）辛钦定理：,设随机变量 相互独立，服从同一分布，且具有数学期望 ，则对于任意正数 ，有,四、,应用,R,软件,模拟验证大数定律,2,、在,R,软件实现的算法思想：,由大数定律可知，当 ，样本的均值趋向与理论分布的期望，因此利用样本容量 逐渐增大这一趋势来模拟 这一趋势，在这种趋势下，样本的均值与理论分布期望的误差 应该呈现出越来越小的趋势，同时，根据上述思想，分别对五种常用分布下的大数定律进行验证。,四、,应用,R,软件,模拟验证大数定律,大数定律模拟算法,设置参数值,产生,m,维序列,绘图,试验次数,是否达到,m,次,是,否,编写,R,程序,选择分布类型,产生随机数,计算样本均值,y,设置循环的跳跃步长 、的第一次抽样的样本容量初始值 和上限值 ；,利用函数 产生由各模拟样本空间大小组成的,m,维序列；,选择随机数 的分布类型，本文中的相关程序仅选择了常用的随机分布：正态分布、指数分布、均匀分布、泊松分布、二项分布、两点分布；,利用,R,软件产生,n,个服从同一分布的随机数 ；,计算 （或 ）的值；,若循环次数,im,是,否,编写,R,程序,产生随机数,计算标准化随机变量,设置参数,j,和,step,中心极限定理模拟算法,选择随机变量 的分布类型，主要分布类型有正态分布、指数分布、均匀分布、泊松分布、二项分布和两点分布；,设置模拟试验总次数,m,及每次模拟试验中随机变量的个数,n,的值,;,利用,R,软件模拟产生,n,个服从同一分布的随机数 ；,使用产生的,n,个随机数计算标准化随机变量值,设置循环变量,j,和循环的跳跃步长 ，当 时，重复步骤、，直至 ；,对,m,个 值进行正态性检验和描述性统计分析，包括直观的,QQ,图检验、正态性,W,检验以及偏度系数、峰度系数、均值和方差。,进入演示,五、应用,R,软件模拟验证中心极限定理,非常感谢！,