1. 程式人生 > >R語言之隨機數與抽樣模擬篇

R語言之隨機數與抽樣模擬篇

3.1 隨機數的產生

3.1.1 均勻分佈隨機數

R語言生成均勻分佈隨機數的函式是runif()

句法是:runif(n,min=0,max=1)    n表示生成的隨機數數量,min表示均勻分佈的下限,max表示均勻分佈的上限;若省略引數min、max,則預設生成[0,1]上的均勻分佈隨機數。

例1:

> runif(5,0,1)     # 生成5個[0,1]的均勻分佈的隨機數
[1] 0.5993 0.7391 0.2617 0.5077 0.7199 

> runif(5)         # 預設生成5個[0,1]上的均勻分佈隨機數
[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)         #新增均勻分佈的密度函式線

3.1.2 正態分佈隨機數

正態分佈隨機數的生成函式是 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)

3.1.3 二項分佈隨機數

二項分佈是指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))

3.1.4  指數分佈隨機數

R生成指數分佈隨機數的函式是:rexp()  

其句法是:rexp(n,lamda=1) n表示生成的隨機數個數,lamda=1/mean

例:

>x=rexp(100,1/10)     # 生成100個均值為10的指數分佈隨機數
>hist(x,prob=T,col=gray(0.9),main=“均值為10的指數分佈隨機數”) 
>curve(dexp(x,1/10),add=T) #新增指數分佈密度線

3.1.5 常見的分佈函式

產生分佈的隨機數,只需要在相應的分佈前加r就行

表 3-1 常見分佈函式表 
分佈  中文名稱 R中的表達  引數
Beta  貝塔分佈 beta(a,b)  shape1,   shape2 
Binomial  二項分佈 binom(n,p)  size,       prob
Cauchy  柯西分佈 cauchy( )  location,   scale  Chi-square  卡方分佈 chisq(df)  df  Exponential  指數分佈 exp(lamda)  rate  F  F分佈 f(df1,df2)  df1         df2
Gamma  伽瑪分佈 gamma()  shape       rate
Geometric  幾何分佈 geom()  prob  Hypergeometric  超幾何分佈 hyper()  m,n,k 
Logistic  邏輯分佈 logis()  location    scale
Negative binomial  負二項分佈 nbinom()  size        prob
Normal  正態分佈 norm()  mean, sd  Multivariate normal  多元正態分佈 mvnorm()  mean,cov 
Poisson  泊松分佈 pois()  lambda  T  t 分佈 t()  df 
Uniform  均勻分佈 unif()  min,       max  Weibull  威布兒分佈 weibull()  shape,     scale 
Wilcoxon  威爾考可森分佈  wilcox()  m,           n 

表 3-2 與分佈相關的函式及代號

函式代號  函式作用
r-  生成相應分佈的隨機數
d-  生成相應分佈的密度函式
p-  生成相應分佈的累積概率密度函式
q-  生成相應分佈的分位數函式

例:

dnorm表示正態分佈密度函式
pnorm表示正態分佈累積概率密度函式
qnorm表示正態分佈分位數函式(即正態累積概率密度函式的逆函式)

3.2  隨機抽樣

3.2.1 放回與無放回抽樣

R可以進行有放回、無放回抽樣

sample()函式即可以實現 

句法為:sample(x,n,replace=F,prob=NULL)

3.3 統計模擬

3.3.1 幾種常見的模擬方法

1 中心極限定理:

2 二項分佈模擬中心極限定理


3 用函式進行模擬

指定模擬次數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 

4 正態概率模擬

能比直方圖更好判定隨機數是否近似服從正態分佈的是正態概率圖。

基本思想是:作實際資料的分位數與正態分佈資料的分位數的散點圖,也就是作樣本分位數與理論分位數的散點圖。

3.3.2 模擬函式的建立方法 

若每次模擬都要編寫一個迴圈,非常麻煩.

sim.fun()就是專門用來解決這類問題的

只需要編寫一個用來生成隨機數的函式,剩下的工作就交給sim.fun來完成

sim.fun <-function (m,f,...)   # 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)) } 

> x=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)-m
u)/(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)