R语言-统计分布和模拟

R语言中统计分布和模拟

前言

  很多应用都需要随机数。像interlink connection,密码系统、视频游戏、人工智能、优化、问题的初始条件,金融等都需要生成随机数。但实际上目前我们并没有“真正”的随机数生成器,尽管有一些伪随机数生成器也是非常有效的。

目录

 1. 概率统计分布概述

 2. 随机函数模拟介绍

 3. 密度函数模拟介绍

 4. 分布函数模拟介绍

 5. 分位数函数模拟介绍

 6. 函数模拟举例

1. 概率统计分布概述

  各种统计分布在R中的名称,这张表取自《An Introduction to R》中概率分布一章,基本涵盖了R中所有的概率函数。

  R给出了详尽的统计表。R 还提供了相关函数来 计算累计概率分布函数 X <= x), 概率密度函数和分位数函数(给定 q,符合 P(X <= x) > q的最小x就是对应的分位数), 和 基于概率分布的计算机模拟。

R中的各种概率统计分布
汉文名称英文名称R对应的名字附加参数
β分布betabetashape1, shape2, ncp
二项式分布binomialbinomsize, prob
柯西分布Cauchycauchylocation, scale
卡方分布chi-squaredchisqdf, ncp
指数分布exponentialexprate
F分布Ffdf1, df1, ncp
Gamma(γ)分布gammagammashape, scale
几何分布geometricgeomprob
超几何分布hypergeometrichyperm, n, k
对数正态分布log-normallnormmeanlog, sdlog
Logistic分布logisticlogislocation, scale
负二项式分布negative binomialnbinomsize, prob
正态分布normalnormmean, sd
泊松分布Poissonpoislambda
Wilcoxon分布signed ranksignrankn
t分布Student's ttdf, ncp
均匀分布uniformunifmin, max
韦伯分布Weibullweibullshape, scale
秩和分布Wilcoxonwilcoxm, n

  概率函数介绍

  在R中各种概率函数都有统一的形式,即一套统一的 前缀+分布函数名

   d 表示密度函数(density);

   p 表示分布函数(生成相应分布的累积概率密度函数);

   q 表示分位数函数,能够返回特定分布的分位数(quantile);

   r 表示随机函数,生成特定分布的随机数(random)。

  每一种分布有四个函数:d―density(密度函数),p―分布函数,q―分位数函数,r―随机数函数。比如,正态分布的这四个函数为dnorm,pnorm,qnorm,rnorm。dnorm 表示正态分布密度函数;pnorm 表示正态分布累积概率密度函数;qnorm 表示正态分布分位数函数(即正态累积概率密度函数的逆函数);rnorm 表示正态分布随机数。各分布后缀,前面加前缀d、p、q或r就构成函数名。

  不同的名字前缀表示不同的含义,d表示概率密度函数,p 表示 累积分布函数(cumulative distribution function,CDF),q 表 示分位函数以及 r 表示随机模拟(random deviates)或者随机数发生器。 dxxx 的第一个参数是xpxxx是qqxxx 是 p,和rxxx的是n(rhyperrwilcox例外,二者的参数是 nn)。偏态指数(non-centrality parameter) ncp 现在仅用于累积分布函数,大多数概率密度函数 和部分其他情况:更细节的内容可以参考帮助文档。

  pxxx 和 qxxx 函数都有逻辑 参数 lower.taillog.pdxxx 也有一个逻辑函数 log。 它们可以用来计算所要的函数值。 例如可以通过下式计算累计(“积分的”) 风险 (hazard)函数。

    - pxxx(t, ..., lower.tail = FALSE, log.p = TRUE)

  它们也可以直接用来计算更精确的对数似然值 (dxxx(..., log = TRUE))。

  此外还有函数 ptukeyqtukey 计算 来自正态分布的样本的标准化全距(studentized range) 的分布。

  这里是一些例子:

    > ## t分布的双侧p值
    > 2*pt(-2.43, df = 13)
    > ## F(2, 7)分布的上1%分位数
    > qf(0.99, 2, 7)

2. 随机函数模拟介绍

  各种分布的随机数生存函数

    rnorm(n, mean=0, sd=1)   #正态分布
    rexp(n, rate=1)   #指数
    rgamma(n, shape, rate=1, scale=1/rate)   #r 分布
    rpois(n, lambda)   #泊松
    rt(n, df, ncp)   #t 分布
    rf(n, df1, df2, ncp)   #f 分布
    rchisq(n, df, ncp=0)   #卡方分布
    rbinom(n, size, prob)   #二项分布
    rweibull(n, shape, scale=1)   #weibull 分布
    rbata(n, shape1, shape2)   #bata 分布

  均匀分布随机数

  R语言生成均匀分布随机数的函数是runif(),句法是:runif(n,min=0,max=1)。 n表示生成的随机数数量,min表示均匀分布的下限,max表示均匀分布的上限;若省略参数min、max,则默认生成[0,1]上的均匀分布随机数。

    # 例1:生成5个[0,1]的均匀分布的随机数
    > runif(5,0,1)     
    [1] 0.5993 0.7391 0.2617 0.5077 0.7199
    # 默认生成5个[0,1]上的均匀分布随机数
    > runif(5)         
    [1] 0.2784 0.7755 0.4107 0.8392 0.7455 
    
    # 例2:随机产生100个均匀分布随机数,作其概率直方图,再添加均匀分布的密度函数线,程序如下:
    > x=runif(100) 
    > hist(x,prob=T,col=gray(.9),main="uniform on [0,1]")
    # 添加均匀分布的密度函数线
    > curve(dunif(x,0,1),add=T)  

  正态分布随机数

  正态分布随机数的生成函数是 rnorm() 。句法是:rnorm(n,mean=0,sd=1)。其中n表示生成的随机数数量,mean是正态分布的均值,默认为0,sd是正态分布的标准差,默认时为1。

    # 例:随机产生100个正态分布随机数,作其概率直方图,再添加正态分布的密度函数线
    
    > x=rnorm(100) 
    > hist(x,prob=T,main="normal mu=0,sigma=1") 
    > curve(dnorm(x),add=T)

  二项分布随机数

  二项分布是指n次独立重复贝努力试验成功的次数的分布,每次贝努力试验的结果只有两个,成功和失败,记成功的概率为p。生成二项分布随机数的函数是:rbinom() 。句法是:rbinom(n,size,prob)。n表示生成的随机数数量,size表示进行贝努力试验的次数,prob表示一次贝努力试验成功的概率。

    # 例:产生100个n为10,15,50,概率p为0.25的二项分布随机数:
    
    > par(mfrow=c(1,3)) 
    > p=0.25 
    > for( n in c(10,20,50)) {  
        x=rbinom(100,n,p) 
        hist(x,prob=T,main=paste("n =",n)) 
        xvals=0:n 
        points(xvals,dbinom(xvals,n,p),type="h",lwd=3) 
      } 
    > par(mfrow=c(1,1))

  指数分布随机数

  R生成指数分布随机数的函数是:rexp()。其句法是:rexp(n,lamda=1)。 n表示生成的随机数个数,lamda=1/mean

    # 例:生成100个均值为10的指数分布随机数
    > x=rexp(100,1/10)     
    > hist(x,prob=T,col=gray(0.9),main="均值为10的指数分布随机数") 
    # 添加指数分布密度线
    > curve(dexp(x,1/10),add=T) 
    
    # 例:生成5个指数分布随机数(应和下面举例)
    > rexp(5, rate=1)
    [1] 0.6626410 1.4266883 0.2150661 1.5788140 0.4469142
    

3. 密度函数模拟介绍

  以指数分布(R中函数名为exp)为例进行示范

  密度函数调用形式:

    dexp(x,rate)

  参数解释:x随机变量,rate为指数概率密度函数的参数λ

    ## 例1:绘制0到4上,参数为1的指数分布的概率密度函数图像
    > x <- seq(0, 4, 0.5)
    > x
    [1] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
    > y <- dexp(x, rate=1)
    > y
    [1] 1.00000000 0.60653066 0.36787944 0.22313016 0.13533528
    [6] 0.08208500 0.04978707 0.03019738 0.01831564
    > plot(x,y)
    > plot(x,y,type='l')

4. 分布函数模拟介绍

  分布函数调用形式:

    pexp(x,rate, lower.tail =TRUE)

  参数解释:x随机变量,rate同上,参数lower.tail为一个逻辑值,TURE表示P(X ≤ x),也是默认值。

    ## 例:求取上图中x=2左侧的概率密度函数曲线下方面积
    > pexp(2, rate=1)
    [1] 0.8646647

5. 分位数函数模拟介绍

  分位数函数调用形式:

    qexp(p,rate, lower.tail =True )

  参数解释:p为概率值,其他同上

    ## 例:求取参数为1的指数分布函数的85%分位数
    > qexp(0.85, rate=1)
    [1] 1.89712

6. 函数模拟举例

  例如:指定模拟次数m=100,样本量n=10,概率=0.25,如果要改变这些参数来重新进行模拟将会很麻烦,下面将展示如何将上面的程序形成一个模拟函数再进行模拟。

    > sim.clt <- function (m=100,n=10,p=0.25) { 
         z = rbinom(m,n,p)                
         x = (z-n*p)/sqrt(n*p*(1-p))         
         hist(x,prob=T,breaks=20,main=paste("n =",n,"p =",p)) 
         curve(dnorm(x),add=T)              
      } 
    > sim.clt()               # 默认 m=100,n=10,p=0.25 
    > sim.clt(1000)           # 取 m=1000,n=10,p=0.25 
    > sim.clt(1000,30)        # 取 m=1000,n=30,p=0.25 
    > sim.clt(1000,30,0.5)    # 取 m=1000,n=30,p=0.5 

  模拟函数的建立方法

  若每次模拟都要编写一个循环,非常麻烦。sim.fun()就是专门用来解决这类问题的。只需要编写一个用来生成随机数的函数,剩下的工作就交给sim.fun来完成。

    # m 模拟样本次数,f需模拟的函数
    sim.fun <-function (m,f,...) { 
        sample <-1:m 
        for (i in 1:m) { 
            sample[i] <-f(...) 
         } 
        sample 
    } 

  正态概率模拟:

  能比直方图更好判定随机数是否近似服从正态分布的是正态概率图。基本思想:作实际数据的分位数与正态分布数据的分位数的散点图,也就是作样本分位数与理论分位数的散点图。

  二项分布模拟:

  先编写一个函数用来生成一个二项分布随机的标准化值。

    > f <- function(n=10,p=0.5){s=rbinom(1,n,p); (s-n*p)/sqrt(n*p*(1-p)) }
    > xf  <- sim.fun(1000,f)       # 模拟1000个二项随机数
    > hist(x,prob=T)

  均匀分布来模拟中心极限定理:

    > f <- function(n=10) { mean(runif(n)-1/2) / (1/sqrt(12*n)) }
    > x <- sim.fun(1000,f)        # 模拟1000个均匀随机数
    > hist(x,prob=T) 

  正态分布:

    > f <- function(n=10,mu=0,sigma=1){ r=rnorm(n,mu,sigma); (mean(r)-mu)/(sigma/sqrt(n)) }
    
    > x <- sim.fun(1000,f)          # 模拟1000个样本量为10的N(0,1)随机数
    > hist(x,breaks=10,prob=T)
    
    > x <- sim.fun(1000,f,30,5,2)   # 模拟1000个样本量为30的N(5,4)随机数
    > hist(x,breaks=10,prob=T)

转载于:https://www.cnblogs.com/cloudtj/articles/5676624.html

  • 3
    点赞
  • 33
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值