专栏首页机器学习与统计学基于R软件的统计模拟

基于R软件的统计模拟

  • 统计模拟的基本概念

(一)统计模拟的定义

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

(二)统计模拟方法

一般地,统计模拟分类如下:

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

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

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

  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→∞这一趋势,在这种趋势下,样本的均值与理论分布期望的误差ε应该呈现出越来越小的趋势,同时,根据上述思想,分别对五种常用分布下的大数定律进行验证。

> #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")
+ }

本文分享自微信公众号 - 机器学习与统计学(tjxj666)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2015-07-25

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 100天搞定机器学习|Day56 随机森林工作原理及调参实战(信用卡欺诈预测)

    前文对随机森林的概念、工作原理、使用方法做了简单介绍,并提供了分类和回归的实例。本期我们重点讲一下:

    统计学家
  • 蒙特卡罗方法入门

    蒙特卡罗方法是一种计算方法。原理是通过大量随机样本,去了解一个系统,进而得到所要计算的值。

    统计学家
  • 重温统计学⑨假设检验

    零假设(Null Hypothesis):零假设是指干预后的总体和当前总体参数之间没有显著性差别。零假设总是等式,通常如下表示:

    统计学家
  • 机器人钻颅器可帮助医生快速、安全地完成手术

    据美国犹他大学报道,在所有关于未来起机器起义的虚构描写中都有这样一种主题:因为我们已经教会机器人很多人类工作的方法,所以实际上就是创造了对我们的弱点无所不知的敌...

    人工智能快报
  • 早报:全球云服务营收到2021年将会达到5540亿美元

    1、全球云服务营收到2021年将会达到5540亿美元 腾讯科技讯 据外媒报道,市场研究公司IDC发布了它的第一份云服务营收预测报告。IDC公司预计,整体而言,...

    用户1335017
  • django+uwsgi+nginx部署

      说明:Linux系统内置了python2.7,如果你的Django项目依赖于Python3,请使用pip3 install django安装Python3环...

    py3study
  • 7天快速掌握MySQL-DAY4

    查询数据时,如果表名很长,使用起来不方便,此时,就可以为表取一个别名,用这个别名来代替表的名称 SELECT * FROM 表名 [AS] 别名; 注意,为表指...

    披头
  • Google将推中文信息流产品,今日头条们请注意~

    关于这款中文信息流产品,跟Google翻译、猜画小歌一脉相承,据称照样基于移动端、大概率是手机App,重点是为中文信息而生。

    量子位
  • 怎么在elementary OS中使用中文输入法

    能找到官网或官网论坛的帖子,那么就是解决问题的一半了。刚好我也搜到了这个帖子, http://elementary.io/answers/how-to-type...

    用户3579639
  • Go 脚本往InfluxDB插入数据

    简单、

扫码关注云+社区

领取腾讯云代金券