基于R軟件實現隨機抽樣及其應用

張效嘉，胡良平,2*

(1.軍事醫學科學院生物醫學統計學咨詢中心，北京 100850；2.世界中醫藥學會聯合會臨床科研統計學專業委員會，北京 100029

基于R軟件實現隨機抽樣及其應用

張效嘉1，胡良平1,2*

(1.軍事醫學科學院生物醫學統計學咨詢中心，北京 100850；2.世界中醫藥學會聯合會臨床科研統計學專業委員會，北京 100029

本文的目的是使讀者快速掌握用R軟件產生服從多種分布規律的隨機數和多種隨機抽樣的方法。通過介紹隨機數生成器的概念、對應的函數、實際生成服從正態分布和均勻分布的隨機數，讀者很容易獲得基于R軟件產生所需要的隨機數；通過介紹有放回與無放回隨機抽樣、系統與分層隨機抽樣、整群與bootstrap隨機抽樣的概念和具體實現，使讀者輕松地基于R軟件實現多種隨機抽樣。事實表明：R軟件易于獲取、易學易用；R軟件功能強大、適用面寬，可方便快捷地解決試驗設計中的隨機抽樣問題。

R軟件；隨機數生成器；正態分布；均勻分布；隨機抽樣；bootstrap隨機抽樣

1 產生隨機數[1-3]

1.1 隨機數及其種類

1.1.1 概述

所謂隨機數，就是它們在數量大小和先后順序等方面，出現的規律是不確定的。然而，人們還可以在一定的前提條件限制下，去生成一系列的隨機數，它們出現的先后順序是隨機的，但它們在整體上卻是服從某種特定規律的。

由此可知，當人們將“前提條件”設定為“均勻分布”，就可生成均勻分布隨機數；當人們將“前提條件”設定為“正態分布”，就可生成正態分布隨機數；……；以此類推，還可生成很多服從其他特定分布規律的隨機數。

1.1.2 隨機數生成器種類

在計算機上，生成隨機數需要基于一定的算法。不同的算法會生成不同的隨機數。R提供了多種隨機數生成器(random number generators，RNG)，默認的RNG是由Makoto Matsumoto與Takuji Nishimura于1997年-1998年提出的RNG，即Mersenne twister隨機數生成器。該RNG的循環周期為219937-1。R還提供了其他幾種RNG，通過在R的控制臺發送命令“help(RNGkind)”，可以獲得有關的詳細信息。

若用戶想把默認的RNG更改為自己指定的某種其他的RNG，可通過在R的控制臺發送命令“RNGkind(kind=“new_rng”)”。注意：這里的“new_rng”只能是下列幾種備選內容：Mersenne(默認項)、Wich、Mars、Super、Knuth-TAOCP-2002、Knuth-TAOCP、L’Ecuyer-CMRG。

1.1.3 確保每批生成的隨機數相同

一般來說，生成隨機數的函數每批生成的隨機數是不同的。但有時，人們需要生成相同批次的隨機數(如在臨床試驗設計和模擬試驗時)，此時，可以在運行某種隨機數函數時，先運行一下set.seed(n)函數(這里n為一個確定的正整數，如1或2等)。

現以生成均勻分布隨機數的runif( )函數為例，在運行此隨機數函數之前先運行set.seed( )函數，可使runif( )函數每批生成的隨機數相同(生成的隨機數的長度應相同)。

> set.seed(1)

> runif(10)

以上為第一次運行，要求生成在[0,1]區間內且服從均勻分布的10個隨機數。

[1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968

[7] 0.94467527 0.66079779 0.62911404 0.06178627

以上是第一次運行輸出的10個隨機數。

> set.seed(1)

> runif(10)

以上為第二次運行，仍要求生成在[0,1]區間內且服從均勻分布的10個隨機數(與前面的10個隨機數完全相同，此處從略)。

顯然，以上兩次運行，生成的在[0,1]區間內且服從均勻分布的10個隨機數完全相同。

1.2 用R生成正態分布隨機數

【例1】試生成100個服從均值為0、標準差為1的標準正態分布的隨機數。

解答：在R中，使用rnorm( )函數可以生成服從正態分布的隨機數。其語法規則是：rnorm(n,mean=k1,sd=k2)，其中，參數n代表需要生成的隨機數的個數、k1為用戶指定的均值、k2為用戶指定的標準差。當k1=0、k2=1時，就是希望生成服從標準正態分布的n(n必須為一個具體的正整數)個隨機數。rnorm(n,mean=0,sd=1)=rnorm(n,0,1)=rnorm(n)。

> r_number1<- rnorm(100,0,1);r_number1

以上語句的目的是生成100個服從均值=0、標準差=1的標準正態分布的隨機數并將其賦值給變量(本質上是一個向量)r_number，然后，將其輸出(輸出結果從略)。

1.3 用R生成均勻分布隨機數

在R中，使用runif( )函數可以生成服從均勻分布的隨機數。其語法規則是：runif(n,min=k1,max=k2)，其中，參數n代表需要生成的隨機數的個數、k1為用戶指定的服從均勻分布的隨機數的下限、k2為用戶指定的服從均勻分布的隨機數的上限。若省略參數min、max，則默認生成[0,1]上的均勻分布隨機數。

【例2】試生成1000個服從下限為2、上限為12的均勻分布的隨機數，并用直方圖展示它們。

解答：為實現題中的目標，可使用下面的語句。

> x<- runif(1000,min=2,max=12)

> hist(x,prob=T,main="uniform distribution(min=2,max=12)")

第一句的目的是生成1000個服從下限為2、上限為12的均勻分布的隨機數。

第二句的目的是繪制這1000個隨機數的直方圖(無法繪制非標準均勻分布密度函數曲線)。產生的直方圖從略。

1.4 用R生成其他分布隨機數

1.4.1 生成指數分布隨機數

生成指數分布隨機數的函數為rexp()，其語法規則是：rexp(n,lamda=1/mean)，這里n代表生成的隨機數的個數，而lamda代表隨機數的平均值的倒數，mean代表隨機數的平均數。

> x<- rexp(1000,1/110);x

以上語句的目的是生成1000個服從均值為110的指數分布的隨機數并將其輸出(輸出結果此處從略)。

1.4.2 生成貝塔分布隨機數

生成貝塔分布隨機數的函數為rbeta()，其語法規則是：rbeta(n,a,b)，這里，n代表生成的隨機數的個數，而a和b代表該分布的兩個非負的實數(a>0，b>0)。

> x<- rbeta(100,2,8);x

以上語句的目的是生成100個服從參數a=2、參數b=8的貝塔分布的隨機數并將其輸出(輸出結果此處從略)。

1.4.3 生成伽瑪分布隨機數

生成伽瑪分布隨機數的函數為rgamma()，其語法規則是：rgamma(n,shape,rate=1,scale=1/rate)，這里，n代表生成的隨機數的個數，而shape代表該分布的形狀參數、scale代表尺度參數，這兩個參數都應該為非負實數(shape>0、scale>0)。注意：rate與scale兩個中只能給定一個的具體值。

> x<- rgamma(100,1,4);x

以上語句的目的是生成100個服從形狀參數shape=1、尺度參數scale=1/4的伽瑪分布的隨機數并將其輸出(輸出結果此處從略)。

說明：在R中，還可生成許多服從其他分布的隨機數，因篇幅所限，此處從略。

2 隨機抽樣[1-4]

2.1 概述

從給定的一個總體中抽取一定數目的樣品，稱為抽樣。若總體中的所有樣品被放在一個黑色的袋子中，樣品的形狀、大小、顏色和重量是相同的，只是上面的編號不同。當人們從中抽取若干個樣品時，就稱為隨機抽樣。這里的“隨機”是指被抽取的樣品是哪些，事先是不知道的。若從n個樣品(假定它們構成一個總體)中一次性抽取m(m

在R中，可以使用sample( )函數來實現無放回與有放回隨機抽樣。其語法規則如下：

Sample(x,n,replace=F或T，prob=NULL)

x代表總體向量，可以是數據、字符、邏輯向量；n代表樣本含量；replace=F，代表無放回隨機抽樣(這是默認的)，replace=T，代表有放回隨機抽樣；prob可以設置各個抽樣單元不同的入樣概率，進行不等概率抽樣。

2.2 有放回隨機抽樣舉例

【例3】假定擲一枚質地均勻的硬幣，出現正面用“H”表示、出現反面用“B”表示，現重復拋硬幣20次，顯示有放回隨機抽樣的結果。

解答：在R中使用下面的語句就可實現前述的目的。

> sample(c("B","H"),20,rep=T)

[1] "B""B""H""B""H""B""H""B""H""H""B""H""B""H""H""H""H""H"

[19] "B""B"

以上是20次試驗的結果，每次試驗的結果要么是“B”、要么是“H”。其中，“B”出現了9次、“H”出現了11次。

【例4】假定擲一枚質地均勻的骰子(6個面上分別有1、2、3、4、5、6個點)，現重復拋骰子20次，顯示有放回隨機抽樣的結果。

解答：在R中使用下面的語句就可實現前述的目的。

> sample(c(1:6),20,rep=T)

[1] 3 2 4 1 5 1 6 1 6 4 3 6 3 6 5 6 4 4 2 2

以上是20次試驗的結果，每次試驗的結果是1到6六個數字中的一個出現。這批試驗的結果表明，1出現了3次、2出現了3次、3出現了3次、4出現了4次、5出現了2次、6出現了4次。

2.3 無放回隨機抽樣舉例

【例5】現有編號為1～200的200位受試對象，希望從中隨機抽取20位。試顯示被隨機抽取的20位受試者的編號。

解答：在R中使用下面的語句就可實現前述的目的。

> sample(200,20,rep=F)

[1] 180 47 91 124 22 141 36 186 193 8 53 67 113 68 154 143 108 59

[19] 123 71

以上是從1～200個編號中無放回地隨機抽取的20個編號。

2.4 系統隨機抽樣舉例

在R中，要想實現分層隨機抽樣或整群隨機抽樣，需要從R軟件中下載sampling包，分別調用其中的strata( )函數與cluster( )函數；與此同時，這個包中還有getdata( )函數，可以讀取隨機抽樣的全部樣品對應的全部信息(否則，只有與抽樣有關的部分信息)。

但是在R中，筆者目前尚未找到專門用于系統隨機抽樣的函數。雖然在下面的分層隨機抽樣的strata( )函數與整群隨機抽樣的cluster( )函數中，都涉及一個參數，即“method=”， 在等號后面可以輸入四種抽樣方法之一，即srswor(無放回隨機抽樣)、srswr(有放回隨機抽樣)、poisson(泊松隨機抽樣)和systematic(系統隨機抽樣)，但這四種具體的抽樣方法都是在前面的分層隨機抽樣或整群隨機抽樣基礎之上的，而不是直接實施這四種具體隨機抽樣的。故此處暫不能進行系統隨機抽樣，若用戶今后從R軟件中發現了某個包中的某些函數可直接進行系統隨機抽樣或泊松隨機抽樣，可通過R中的幫助功能，很容易學會使用。

2.5 分層隨機抽樣舉例

【例6】試以文獻[5]第208-209頁“表9-1”的資料作為“抽樣框”，并以“sex”為分層因素，實施分層隨機抽樣。

解答：在R中，實現分層隨機抽樣的strata( )函數被放置在sampling包中，使用前必須先下載并導入到R軟件中。操作的步驟如下：①將安裝了R軟件的計算機聯網并啟動R軟件；②選擇“程序包(Packages)”并選擇一個鏡像(Set CRAN mirror…)；③選擇“安裝程序包(Install package(s)…)”；④在彈出的長條形“程序包(packages)”窗口內(約有數千行，其中少部分是“程序包”，絕大多數是“函數”)瀏覽sampling；⑤在R中使用下面的語句就可實現前述的目的。

在strata( )函數中未明確寫出參數method="srswor"，系統默認為“分層(無放回)隨機抽樣”。

> data1<- read.table("F:/studyr/sex_age.txt",header=T)

> sub1<- strata(data1,stratanames="sex",size=c(4,5),description=TRUE)

第一句表明，從F盤上文件夾為studyr內讀取文本文件sex_age.txt賦值給變量(本質上是數據框)data1。

第二句表明，調用strata( )函數，采用默認的隨機抽樣方法(即分層隨機抽樣)；設置分層因素為sex；分別從F層與M層抽取4例與5例受試對象或樣品；需要對隨機抽樣出來的樣品進行詳細描述“description=TRUE”。

Stratum 1

Population total and number of selected units: 12 4

Stratum 2

Population total and number of selected units: 8 5

以上結果表明，從第1層(即F層)12個樣品中隨機抽取了4個樣品；從第2層(即M層)8個樣品中隨機抽取了5個樣品。

Number of strata 2

Total number of selected units 9

以上結果表明，總體中總共有兩層，共隨機抽取9個樣品。

Warning message:

In strata(data1, stratanames = "sex", size = c(4, 5), description = TRUE) :

the method is not specified; by default, the method is srswor

以上結果表明，系統給出了“警告信息”，在調用strata( )函數時，關于抽樣方法的參數未指定，作為默認方法，采取的是分層(無放回)隨機抽樣(srswor)。

> sub1

該語句表明，要求系統呈現出分層隨機抽樣的結果sub1中的與抽樣有關的樣品信息。

SexIDunitProbStratum1F10.333333314F40.3333333116F160.3333333119F190.333333317M70.625000028M80.6250000210M100.6250000211M110.6250000218M180.62500002

以上結果表明，第1列為被抽取的樣品的編號、第2列為分層因素sex的水平代碼、第3列為被抽取的樣品在總體中的編號、第4列為各行上的樣品被抽取的概率(F層中總共12個樣品，抽取4個，每一個被抽取的概率為4/12=1/3=0.333333；M層中總共8個樣品，抽取5個，每一個被抽取的概率為5/8=0.625)、第5列為層的編號。

>getdata(data1,sub1)

這一句的目的是調用getdata( )函數，顯示被抽取的樣本中每個樣品的全部信息(即各變量及其取值)。

idagesexIDunitProbStratum1160F10.333333314457F40.33333331161633F160.33333331191949F190.333333317760M70.625000028864M80.62500002101016M100.62500002111158M110.62500002181840M180.62500002

以上結果表明，與前面的輸出結果相比較，增加了“年齡(age)”及其取值。如果原數據集(或抽樣總體)中除分層因素外，還有50個變量及其取值，此時，都將全部顯示出來。

在strata( )函數中明確寫出參數method="srswor"，指明要進行分層(無放回)隨機抽樣。

> data1<- read.table("F:/studyr/sex_age.txt",header=T)

> sub2<- strata(data1,stratanames="sex",size=c(4,5),method="srswor",description=TRUE)

與前面的程序相比較，多了參數method="srswor"。

> sub2

SexIDunitProbStratum4F40.333333316F60.3333333117F170.3333333119F190.333333313M30.625000027M70.625000028M80.6250000211M110.6250000212M120.62500002

以上結果表明，被抽取的樣品編號發生了改變。

>getdata(data1,sub2)

idagesexIDunitProbStratum4457F40.333333316631F60.33333331171739F170.33333331191949F190.333333313337M30.625000027760M70.625000028864M80.62500002111158M110.62500002121263M120.62500002

以上結果表明，分層隨機抽樣的結果發生了變化。

2.6 整群隨機抽樣舉例

2.6.1 數據集介紹

在R軟件的MASS包中，有一個名為Insurance的數據集。其內記錄了某保險公司1973年第三季度車險投保人的相關信息。其各觀測(即行)代表投保人，5列的列名分別為：第1列為District(投保人家庭住址所在區域，代號為1至4)；第2列為Group(所投保汽車的發動機排放量，其各水平分別為：<1升、1.0～1.5升、>1.5～2升、>2升共四個等級)；第3列為Age[投保人的年齡，其各年齡檔為：<25歲、25～30歲(不含30歲)、30～35歲、>35歲]；第4列為Holders(投保人數量)；第5列為Claims(要求索賠的投保人數量)。

2.6.2 使用sampling包中的cluster( )函數進行整群隨機抽樣

> sub3<- cluster(Insurance,clustername="District",size=2,method="srswor",description=TRUE)

以上語句的目的是調用cluster( )函數對Insurance數據集進行整群隨機抽樣，群的標志為“District”(即區域)、隨機抽取的群數為2、將采取簡單(無放回)隨機抽樣(即srswor)方法從4個群中隨機地抽取兩個群、希望對被抽取的樣品信息進行描述(即description=TRUE)。

Number of selected clusters: 2

Population total and number of selected units 64 32

以上結果表明，隨機抽取了兩個群，共有64個樣品，而被抽取的樣品數為32個。

> getdata(Insurance,sub3)

以上語句表明，調用getdata( )函數顯示從Insurance數據集中抽取的兩個群中的全部樣品。

GroupAgeHoldersClaimsDistrictIDunitProb1<1L<2519738110.52<1L25～2926435120.551～1.5L<2528463150.561～1.5L25～2953684160.5……18<1L25～29139192180.519<1L30～35151222190.517<1L<2585222170.5261.5～2L25～29175462260.5……

以上是被抽取的32個樣品，為節省篇幅，兩個群內都只保留了4個樣品，其他樣品被省略了。由于共有4個群，因此，抽取兩個群時，對于每個群而言，被抽中的概率都為0.5。

2.7 bootstrap隨機重抽樣舉例

2.7.1 概述

bootstrap隨機重抽樣是以原始數據為基礎的模擬抽樣統計推斷法，其基本思想是：在原始數據范圍(假定總樣本含量為n)內做有放回的再抽樣，每次抽取的樣本含量都是m，原始數據中每個觀測(樣品或個體)每次被抽中的概率相等，即1/n，所抽得的樣本被稱為bootstrap樣本。

bootstrap樣本有何用途呢？當人們無法駕馭總體但可以反復獲得bootstrap樣本時，由它們所呈現的規律就非常接近無法駕馭的那個總體。下面用一個可以駕馭的總體為例，先呈現出總體中一個定量變量的取值的頻數分布情況，再采取有放回的隨機抽樣，獲得一個樣本含量較大(如1000或10000)的bootstrap樣本，再呈現該bootstrap樣本的頻數分布。比較這兩個頻數分布的形狀是否接近。

2.7.2 用于抽取bootstrap樣本的數據集

在R中，有一個內置的faithful數據集(Old Faithful Geyser Data)，它記錄了美國懷俄明州黃石國家公園內的古老而又忠實的噴泉的噴射持續時間eruptions(min)和等待時間waiting(min)(即相鄰兩次噴射之間的時間間隔)。

> faithful

此語句用于顯示faithful數據集中的全部數據(n=272)。

eruptionswaiting13.6007921.8005433.33374……2704.417902711.817462724.46774……

以上顯示出faithful數據集中的開始和結尾的各3個觀測結果。

2.7.3 顯示總體中eruptions變量的頻數分布

> attach(faithful)

> hist(eruptions,breaks=25)

第一句的目的是綁定數據；第二句的目的是繪制eruptions變量的頻數直方圖，因篇幅所限，圖形從略。

2.7.4 從總體中進行2000次有放回的隨機抽樣，獲得樣本含量n=2000的bootstrap樣本

> boot.sample<- sample(eruptions,2000,rep=T)

以上語句的目的就是為了實現從總體中進行2000次有放回的隨機抽樣，獲得樣本含量n=2000的bootstrap樣本。

2.7.5 將總體與bootstrap樣本中的eruptions變量的頻數分布直方圖平行地繪制出來

> par(mfrow=c(1,2))

> hist(eruptions,breaks=25)

> hist(boot.sample,breaks=25)

> par(mfrow=c(1,1))

以上四個語句的目的分別為：

第一句：設置作圖窗口為一行兩列(就是希望把兩幅圖左右放置)。

第二句：在左邊放置總體中的eruptions變量的頻數分布直方圖。

第三句：在右邊放置bootstrap樣本中的eruptions變量的頻數分布直方圖。

第四句：設置作圖窗口為一行一列(就是把前面的一行上繪制兩幅圖的設置改回默認狀態，即一行上只繪制一幅圖)。

由圖形繪制的結果(因篇幅所限，此處從略)可以清楚地看出：bootstrap樣本呈現的圖形幾乎與總體一模一樣。若總體真是不可駕馭的，那就可以通過bootstrap樣本呈現的形狀去推斷它。

[1] 黃文, 王正林. 數據挖掘: R語言實戰[M]. 北京: 電子工業出版社, 2015: 34-39.

[2] 李詩羽, 張飛, 王正林. 數據分析: R語言實戰[M]. 北京: 電子工業出版社, 2015: 88-156.

[3] 方匡南, 朱建平, 姜葉飛. R數據分析: 方法與案例詳解[M]. 北京: 電子工業出版社, 2015: 54-168.

[4] Joseph Adler. R語言核心技術手冊[M]. 2版. 劉思喆, 李艦, 陳鋼, 等譯. 北京: 電子工業出版社, 2015: 417-421.

[5] 胡良平. 科研設計與統計分析[M]. 北京: 軍事醫學科學出版社, 2012:129-274.

(本文編輯：陳 霞)

The realization of the random sampling and its application based on R software

ZhangXiaojia1,HuLiangping1,2*

(1.ConsultingCenterofBiomedicalStatistics,AcademyofMilitaryMedicalSciences,Beijing100850,China; 2.SpecialtyCommitteeofClinicalScientificResearchStatisticsofWorldFederationofChineseMedicineSocieties,Beijing100029,China

*Correspondingauthor:HuLiangping,E-mail:lphu812@sina.com)

The paper aims to make it convenient for readers to utilize R software to generate random number from various distributions and implement varieties of random sampling methods. The paper introduced the conception and corresponding functions of random number generators. We presented cases of generating random number from both normal distribution and uniform distribution. Consequently, readers may easily generate random number by R software. The paper also introduced the conception and realization of random sampling both with and without replacement, systematic random sampling, stratified random sampling, cluster random sampling and bootstrap random sampling. Therefore, readers may use R software to realize varieties of random sampling methods. The fact indicated that R software may be simply obtained and utilized. R software may be used to solve random sampling problem in experimental designs conveniently, due to its powerful function and wide application.

R software; Random number generator; Normal distribution; Uniform distribution; Random sampling; Bootstrap random sampling

*通信作者：胡良平，E-mail：lphu812@sina.com)

R195.1

A

10.11886/j.issn.1007-3256.2016.06.002

國家高技術研究發展計劃課題資助(2015AA020102)

2016-12-03)