基于R软件的统计模拟
- 统计模拟的基本概念
（一）统计模拟的定义
统计模拟即是计算机统计模拟，它实质上是计算机建模，而这里的计算机模型就是计算机方法、统计模型(如程序、流程图、算法等)，它是架于计算机理论和实际问题之间的桥梁。它与统计建模的关系如下图。
（二）统计模拟方法
一般地，统计模拟分类如下：
若按状态变量的变化性质分为连续随机模拟和离散随机模拟。
而按变量是否随时间变化又可分为动态随机模拟和静态随机模拟。
常用的统计模拟方法主要有以下几种：
- 1.蒙特卡罗法
- 2.系统模拟方法
- 3.其它方法：包括Bootstrap(自助法)、MCMC（马氏链蒙特卡罗法）等。
（三）统计模拟的一般步骤
- 赶火车问题
一列列车从A站开往B站，某人每天赶往B站上车。他已经了解到火车从A站到B站的运行时间是服从均值为30min，标准差为2min的正态随机变量。火车大约下午13：00离开A站，此人大约13:30到达B站。火车离开A站的时刻及概率如表1所示，此人到达B站的时刻及概率如表2所示。问此人能赶上火车的概率有多大？
——问题的分析——
这个问题用概率论的方法求解十分困难，它涉及此人到达时刻、火车离开站的时刻、火车运行时间几个随机变量，而且火车运行时间是服从正态分布的随机变量，没有有效的解析方法来进行概率计算。在这种情况下可以用计算机模拟的方法来解决。
à为了便于建模，对模型中使用的变量作出如下假定：
à为了分析简化，假定13时为时刻t=0，则变量 、 的分布律为：
- 此人能及时赶上火车的充分必要条件为：
，所以此人能赶上火车的概率模型为：
。
> windows(7, 3) > prb = replicate(5,{ + 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 <= r) { + i = i + 1 + segments(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...missed the train!", "Bingo! + catched thetrain!")) + Sys.sleep(1) + x + r > y + }) > mean(prb) [1] 0.4
三、R软件的统计模拟功能
1、R软件优秀的随机数模拟功能
生产某概率分布的随机数是实现统计模拟的前提条件，而使用R命令可以生成以下常用分布的随机数
2、优良的编程环境和编程语言
R所拥有的好的兼容性、拓展性和强大的内置函数有利于统计模拟的实现。
3、高效率的向量运算功能
使用R拥有的向量运算功能可以大大减少程序运行的时间，提高程序运行的效率。
- 应用R软件模拟验证大数定律
2、在R软件实现的算法思想：
由大数定律可知，当n→∞，样本的均值趋向与理论分布的期望，因此利用样本容量 逐渐增大这一趋势来模拟n→∞这一趋势，在这种趋势下，样本的均值与理论分布期望的误差ε应该呈现出越来越小的趋势，同时，根据上述思想，分别对五种常用分布下的大数定律进行验证。
> #n1为循环的初始值> #n2为循环的上限值，step为步长> #注意parameter是一个向量，其中第一个参数为均值> dashu<-function(n1,n2,steps,epesino,types,parameter){ + #计算需模拟的数据集+ datas<-seq(n1,n2,steps) + #通过switch语句选择理论分布的类型并调用相应类型的模拟子函数+ switch(types,normal_d(datas,parameter,epesino),exponential_d(datas,parameter,epesino),+ uniform_d(datas,parameter,epesino), poisson_d(datas,parameter,epesino),+ binomial_d(datas,parameter,epesino)) + + }> > #正态分布条件下的大数定律模拟函数 > normal_d<-function(datas,parameter,epesino){+ le<-length(datas) #计算向量datas的长度+ result<-1:le #结果向量的初始化+ for(i in 1:le) #模拟次数，反映出样本空间越来越大这一趋势+ {+ #parameter[1]代表均值，parameter[2]代表标准差+ temp<-rnorm(datas[i],parameter[1],parameter[2]) #产生datas[i]个独立正态分布的随机数+ result[i]<-mean(temp) #计算抽样的均值+ }+ result #显示抽样计算结果+ #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图+ plot(datas,result,type="l")+ #作出平均值线，以反映出抽样平均值围绕总体均值波动的规律+ abline(h=parameter[1])+ abline(h=parameter[1]+epesino)+ abline(h=parameter[1]-epesino)+ title(main="normal sitimulation") #给所绘图加标题+ }> > #指数分布条件下大数定律模拟函数> exponential_d<-function(datas,parameter,epesino){+ le<-length(datas) #计算向量datas的长度+ result<-1:le #结果向量的初始化+ for(i in 1:le) #模拟次数，反映出样本空间越来越大这一趋势+ {+ #parameter[1]代表lampta+ temp<-rexp(datas[i],parameter[1]) #产生datas[i]个独立指数分布的随机数+ result[i]<-mean(temp) #计算抽样的均值+ }+ result #显示抽样计算结果+ #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图+ plot(datas,result,type="l")+ #作出平均值线，以反映出抽样平均值围绕总体均值波动的规律+ abline(h=1/parameter[1])+ abline(h=1/parameter[1]+epesino)+ abline(h=1/parameter[1]-epesino)+ title(main="exponential sitimulation")+ }> > > #均匀分布条件下大数定律模拟函数> uniform_d<-function(datas,parameter,epesino){+ le<-length(datas) #计算向量datas的长度+ result<-1:le #结果向量的初始化+ for(i in 1:le) #模拟次数，反映出样本空间越来越大这一趋势+ {+ #parameter[1]代表下限，parameter[2]代表上限+ temp<-runif(datas[i],parameter[1],parameter[2]) #产生datas[i]个独立均匀分布的随机数+ result[i]<-mean(temp) #计算抽样的均值+ }+ result #显示抽样计算结果+ #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图+ plot(datas,result,type="l")+ #作出平均值线，以反映出抽样平均值围绕总体均值波动的规律+ abline(h=(parameter[1]+parameter[2])/2)+ abline(h=(parameter[1]+parameter[2])/2+epesino)+ abline(h=(parameter[1]+parameter[2])/2-epesino)+ title(main="uniform sitimulation")+ }> > #泊松分布条件下大数定律模拟函数> poisson_d<-function(datas,parameter,epesino){+ le<-length(datas) #计算向量datas的长度+ result<-1:le #结果向量的初始化+ for(i in 1:le) #模拟次数，反映出样本空间越来越大这一趋势+ {+ #parameter[1]代表均值lambda + temp<-rpois(datas[i],parameter[1]) #产生datas[i]个独立泊松分布的随机数+ result[i]<-mean(temp) #计算抽样的均值+ }+ result #显示抽样计算结果+ #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图+ plot(datas,result,type="l")+ abline(h=parameter[1])+ abline(h=parameter[1]+epesino)+ abline(h=parameter[1]-epesino)+ title(main="Poisson sitimulation")+ }> > #二项分布条件下大数定律模拟函数> binomial_d<-function(datas,parameter,epesino){+ le<-length(datas) #计算向量datas的长度+ result<-1:le #结果向量的初始化+ for(i in 1:le) #模拟次数，反映出样本空间越来越大这一趋势+ {+ #parameter[1]代表n即实验的次数，parameter[2]代表p,即每次成功的概率+ temp<-rbinom(datas[i],parameter[1],parameter[2]) #产生datas[i]个独立二项分布的随机数+ result[i]<-mean(temp) #计算抽样的均值+ }+ result #显示抽样计算结果+ #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图+ plot(datas,result,type="l")+ abline(h=parameter[1]*parameter[2])+ abline(h=parameter[1]*parameter[2]+epesino)+ abline(h=parameter[1]*parameter[2]-epesino)+ title(main="binominal sitimulation")+ }腾讯云开发者
