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中的各种概率统计分布
汉⽂名称
β分布⼆项式分布柯西分布卡⽅分布指数分布F分布⼏何分布超⼏何分布对数正态分布Logistic分布负⼆项式分布正态分布泊松分布Wilcoxon分布t分布均匀分布韦伯分布秩和分布英⽂名称
betabinomialCauchychi-squaredexponentialFgeometrichypergeometriclog-normallogisticnormalPoissonsigned rankStudent's tuniformWeibullWilcoxon
R对应的名字
betabinomcauchychisqexpfgammageomhyperlnormlogisnormpoissignranktunifweibullwilcox
附加参数
shape1, shape2, ncpsize, problocation, scaledf, ncpratedf1, df1, ncpshape, scaleprobm, n, kmeanlog, sdloglocation, scalesize, probmean, sdlambdandf, ncpmin, maxshape, scalem, n
Gamma(γ)分布gamma
negative binomialnbinom
概率函数介绍
在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 的第⼀个参数是x，pxxx是q， qxxx 是 p，和rxxx的是n(rhyper 和 rwilcox例外，⼆者的参数是 nn)。偏态指数（non-centrality parameter） ncp 现在仅⽤于累积分布函数，⼤多数概率密度函数 和部分其他情况：更细节的内容可以参考帮助⽂档。
pxxx 和 qxxx 函数都有逻辑 参数 lower.tail 和 log.p。dxxx 也有⼀个逻辑函数 log。 它们可以⽤来计算所要的函数值。 例如可以通过下式计算累计(“积分的”) 风险 （hazard）函数。
- pxxx(t, ..., lower.tail = FALSE, log.p = TRUE)
它们也可以直接⽤来计算更精确的对数似然值 (dxxx(..., log = TRUE))。
此外还有函数 ptukey 和 qtukey 计算 来⾃正态分布的样本的标准化全距（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 =\ xvals=0:n
points(xvals,dbinom(xvals,n,p),type=\"h\ }
> 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.66210 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.018315 > 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.867
5. 分位数函数模拟介绍
分位数函数调⽤形式：
qexp（p，rate, lower.tail =True ）
参数解释：p为概率值，其他同上
## 例：求取参数为1的指数分布函数的85%分位数 > qexp(0.85, rate=1) [1] 1.712
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 =\ 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)
