基于R软件的统计模拟

2019-04-10 10:53:47 浏览数 (1)
  • 统计模拟的基本概念

（一）统计模拟的定义

统计模拟即是计算机统计模拟，它实质上是计算机建模，而这里的计算机模型就是计算机方法、统计模型(如程序、流程图、算法等)，它是架于计算机理论和实际问题之间的桥梁。它与统计建模的关系如下图。

（二）统计模拟方法

一般地，统计模拟分类如下：

若按状态变量的变化性质分为连续随机模拟和离散随机模拟。

而按变量是否随时间变化又可分为动态随机模拟和静态随机模拟。

常用的统计模拟方法主要有以下几种：

  1. 1.蒙特卡罗法
  2. 2.系统模拟方法
  3. 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→∞这一趋势，在这种趋势下，样本的均值与理论分布期望的误差ε应该呈现出越来越小的趋势，同时，根据上述思想，分别对五种常用分布下的大数定律进行验证。

代码语言：javascript复制
> #n1为循环的初始值
代码语言：javascript复制
> #n2为循环的上限值，step为步长
代码语言：javascript复制
> #注意parameter是一个向量，其中第一个参数为均值
代码语言：javascript复制
> dashu<-function(n1,n2,steps,epesino,types,parameter){
代码语言：javascript复制
       #计算需模拟的数据集
代码语言：javascript复制
       datas<-seq(n1,n2,steps)
代码语言：javascript复制
       #通过switch语句选择理论分布的类型并调用相应类型的模拟子函数
代码语言：javascript复制
       switch(types,normal_d(datas,parameter,epesino),exponential_d(datas,parameter,epesino),
代码语言：javascript复制
       uniform_d(datas,parameter,epesino), poisson_d(datas,parameter,epesino),
代码语言：javascript复制
       binomial_d(datas,parameter,epesino))
代码语言：javascript复制
代码语言：javascript复制
  }
代码语言：javascript复制
>
代码语言：javascript复制
> #正态分布条件下的大数定律模拟函数
代码语言：javascript复制
> normal_d<-function(datas,parameter,epesino){
代码语言：javascript复制
        le<-length(datas)     #计算向量datas的长度
代码语言：javascript复制
        result<-1:le          #结果向量的初始化
代码语言：javascript复制
        for(i in 1:le)      #模拟次数，反映出样本空间越来越大这一趋势
代码语言：javascript复制
        {
代码语言：javascript复制
        #parameter[1]代表均值，parameter[2]代表标准差
代码语言：javascript复制
        temp<-rnorm(datas[i],parameter[1],parameter[2])  #产生datas[i]个独立正态分布的随机数
代码语言：javascript复制
        result[i]<-mean(temp)  #计算抽样的均值
代码语言：javascript复制
             }
代码语言：javascript复制
        result                 #显示抽样计算结果
代码语言：javascript复制
        #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图
代码语言：javascript复制
        plot(datas,result,type="l")
代码语言：javascript复制
        #作出平均值线，以反映出抽样平均值围绕总体均值波动的规律
代码语言：javascript复制
        abline(h=parameter[1])
代码语言：javascript复制
        abline(h=parameter[1] epesino)
代码语言：javascript复制
        abline(h=parameter[1]-epesino)
代码语言：javascript复制
        title(main="normal sitimulation")         #给所绘图加标题
代码语言：javascript复制
  }
代码语言：javascript复制
>
代码语言：javascript复制
> #指数分布条件下大数定律模拟函数
代码语言：javascript复制
> exponential_d<-function(datas,parameter,epesino){
代码语言：javascript复制
       le<-length(datas)     #计算向量datas的长度
代码语言：javascript复制
        result<-1:le          #结果向量的初始化
代码语言：javascript复制
        for(i in 1:le)      #模拟次数，反映出样本空间越来越大这一趋势
代码语言：javascript复制
        {
代码语言：javascript复制
        #parameter[1]代表lampta
代码语言：javascript复制
        temp<-rexp(datas[i],parameter[1])  #产生datas[i]个独立指数分布的随机数
代码语言：javascript复制
        result[i]<-mean(temp)  #计算抽样的均值
代码语言：javascript复制
             }
代码语言：javascript复制
        result                 #显示抽样计算结果
代码语言：javascript复制
        #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图
代码语言：javascript复制
        plot(datas,result,type="l")
代码语言：javascript复制
        #作出平均值线，以反映出抽样平均值围绕总体均值波动的规律
代码语言：javascript复制
        abline(h=1/parameter[1])
代码语言：javascript复制
        abline(h=1/parameter[1] epesino)
代码语言：javascript复制
        abline(h=1/parameter[1]-epesino)
代码语言：javascript复制
        title(main="exponential sitimulation")
代码语言：javascript复制
  }
代码语言：javascript复制
>
代码语言：javascript复制
>
代码语言：javascript复制
> #均匀分布条件下大数定律模拟函数
代码语言：javascript复制
> uniform_d<-function(datas,parameter,epesino){
代码语言：javascript复制
        le<-length(datas)     #计算向量datas的长度
代码语言：javascript复制
        result<-1:le          #结果向量的初始化
代码语言：javascript复制
        for(i in 1:le)      #模拟次数，反映出样本空间越来越大这一趋势
代码语言：javascript复制
        {
代码语言：javascript复制
        #parameter[1]代表下限，parameter[2]代表上限
代码语言：javascript复制
        temp<-runif(datas[i],parameter[1],parameter[2])  #产生datas[i]个独立均匀分布的随机数
代码语言：javascript复制
        result[i]<-mean(temp)  #计算抽样的均值
代码语言：javascript复制
             }
代码语言：javascript复制
        result                 #显示抽样计算结果
代码语言：javascript复制
         #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图
代码语言：javascript复制
        plot(datas,result,type="l")
代码语言：javascript复制
        #作出平均值线，以反映出抽样平均值围绕总体均值波动的规律
代码语言：javascript复制
        abline(h=(parameter[1] parameter[2])/2)
代码语言：javascript复制
        abline(h=(parameter[1] parameter[2])/2 epesino)
代码语言：javascript复制
        abline(h=(parameter[1] parameter[2])/2-epesino)
代码语言：javascript复制
        title(main="uniform sitimulation")
代码语言：javascript复制
  }
代码语言：javascript复制
>
代码语言：javascript复制
> #泊松分布条件下大数定律模拟函数
代码语言：javascript复制
> poisson_d<-function(datas,parameter,epesino){
代码语言：javascript复制
        le<-length(datas)     #计算向量datas的长度
代码语言：javascript复制
        result<-1:le          #结果向量的初始化
代码语言：javascript复制
        for(i in 1:le)      #模拟次数，反映出样本空间越来越大这一趋势
代码语言：javascript复制
        {
代码语言：javascript复制
        #parameter[1]代表均值lambda
代码语言：javascript复制
        temp<-rpois(datas[i],parameter[1])  #产生datas[i]个独立泊松分布的随机数
代码语言：javascript复制
        result[i]<-mean(temp)  #计算抽样的均值
代码语言：javascript复制
             }
代码语言：javascript复制
        result                 #显示抽样计算结果
代码语言：javascript复制
        #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图
代码语言：javascript复制
        plot(datas,result,type="l")
代码语言：javascript复制
        abline(h=parameter[1])
代码语言：javascript复制
        abline(h=parameter[1] epesino)
代码语言：javascript复制
        abline(h=parameter[1]-epesino)
代码语言：javascript复制
        title(main="Poisson sitimulation")
代码语言：javascript复制
  }
代码语言：javascript复制
>
代码语言：javascript复制
> #二项分布条件下大数定律模拟函数
代码语言：javascript复制
> binomial_d<-function(datas,parameter,epesino){
代码语言：javascript复制
        le<-length(datas)     #计算向量datas的长度
代码语言：javascript复制
        result<-1:le          #结果向量的初始化
代码语言：javascript复制
        for(i in 1:le)      #模拟次数，反映出样本空间越来越大这一趋势
代码语言：javascript复制
        {
代码语言：javascript复制
        #parameter[1]代表n即实验的次数，parameter[2]代表p,即每次成功的概率
代码语言：javascript复制
        temp<-rbinom(datas[i],parameter[1],parameter[2])  #产生datas[i]个独立二项分布的随机数
代码语言：javascript复制
        result[i]<-mean(temp)  #计算抽样的均值
代码语言：javascript复制
             }
代码语言：javascript复制
        result                 #显示抽样计算结果
代码语言：javascript复制
        #x轴代表样本容量，y轴代表每次抽样所得的样本的平均值，做序列图
代码语言：javascript复制
        plot(datas,result,type="l")
代码语言：javascript复制
        abline(h=parameter[1]*parameter[2])
代码语言：javascript复制
        abline(h=parameter[1]*parameter[2] epesino)
代码语言：javascript复制
        abline(h=parameter[1]*parameter[2]-epesino)
代码语言：javascript复制
        title(main="binominal sitimulation")
代码语言：javascript复制
  }
r语言

0 人点赞