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對應的名字 | 附加參數 |
| β分布 | beta | beta | shape1, shape2, ncp |
| 二項式分布 | binomial | binom | size, prob |
| 柯西分布 | Cauchy | cauchy | location, scale |
| 卡方分布 | chi-squared | chisq | df, ncp |
| 指數分布 | exponential | exp | rate |
| F分布 | F | f | df1, df1, ncp |
| Gamma(γ)分布 | gamma | gamma | shape, scale |
| 幾何分布 | geometric | geom | prob |
| 超幾何分布 | hypergeometric | hyper | m, n, k |
| 對數正態分布 | log-normal | lnorm | meanlog, sdlog |
| Logistic分布 | logistic | logis | location, scale |
| 負二項式分布 | negative binomial | nbinom | size, prob |
| 正態分布 | normal | norm | mean, sd |
| 泊松分布 | Poisson | pois | lambda |
| Wilcoxon分布 | signed rank | signrank | n |
| t分布 | Student's t | t | df, ncp |
| 均勻分布 | uniform | unif | min, max |
| 韋伯分布 | Weibull | weibull | shape, scale |
| 秩和分布 | Wilcoxon | wilcox | m, 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 的第一個參數是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 =",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)