前言
在Matlab、R或者S-PLUS等軟件中做隨機數模擬時,經過會遇到set.seed()這個函數。隨機數的產生需要有一個隨機的種子,因為用計算機產生的隨機數是通過遞推的方法得來的,必須有一個初始值。用同一台電腦,且在初始值和遞推方法相同的情況下,可以產生相同的隨機序列。
用計算機產生的是“偽隨機數”。用投色子計數的方法產生真正的隨機數 , 但電腦若也這樣做 , 將會占用大量內存;用噪聲發生器或放射性物質也可產生真正的隨機數,但不可重復。而用數學方法產生最適合計算機,這就是周期有限,易重復的 “偽隨機數”。
set.seed函數括號內的參數是一個數,可能會經常看到該數取不同的值。這個函數的目的主要是為了產生相同的隨機數,這樣一來,每次運行的結果就是一樣的,這使得結果可重復被驗證。參數的選擇很隨意,如2、23、124等等。
為了驗證上面的結論,我分別取 set.seed(2)、 set.seed(2)、 set.seed(3),生成一組隨機變量,對生成的隨機變量作圖,由此驗證了前兩組隨機變量在種子數相同的情況下,得到的值是一致的(黑、紅),而第三組情況(藍)則與前兩者不同。
目錄
1. 隨機數生成
2. 抽樣模擬
3. 設置隨機數種子
1. 隨機數生成
作為一種統計分析語言,R是一個生成各種統計分布功能隨機數的綜合性圖書館。在這篇文章中,我想專注於這個簡單的問題:我如何生成一個隨機數?
答案取決於你想要什么樣的隨機數生成?讓我們通過例子說明。
a) 生成 5.0 和 7.5 之間的隨機數
如果你想生成一個十進制數規定的最低和最高之間的任何值(包括分數值)同樣是可能的,使用runif()功能。 這個函數生成均勻分布的值。 這里是如何生成一個 5.0 和 7.5 之間的隨機數的方法:
> x1 <- runif(1, 5.0, 7.5)
> x1
[1] 5.573882
當然,當你運行這個,你會得到一個不同的數字,但它一定會在 5.0 和 7.5 之間。 你不會得到准確值 5.0 或 7.5 。
如果你想生成多個隨機的值,不要使用一個循環。您可以生成多個值一次通過指定您要作為第一個參數runif()值的數目。這里是如何產生10個 5.0 和 7.5 之間的值:
> x2 <- runif(10, 5.0, 7.5)
> x2
[1] 6.871893 6.518833 5.372662 6.808643 5.490282 7.381824
[7] 7.476985 5.569176 6.182591 5.285149
b) 生成 1 到 10 之間的隨機整數
這看起來像最后一個相同的運動,但現在我們只希望完整的數字,而不是分數值。 為此,我們使用的示例函數:
> x3 <- sample(1:10, 1)
> x3
[1] 3
c) 各種分布的隨機數生存函數
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 分布
2. 抽樣模擬
a) 在 1 到 10 之間隨機抽取 1 個數
> x1 <- sample(1:10, 1)
> x1
[1] 3
第一個參數是一個有效的數字向量生成(這里的數字1到10),第二個參數表示應返回一個數字。
b) 在 “一組數” 之間隨機抽取 多 個數
如果我們要生成多個隨機數,我們必須增加一個額外的參數,表示允許重復:
# 有放回抽取
> x2 <- sample(1:10, 5, replace=T)
> x2
[1] 1 5 9 9 8
如果我們要生成多個隨機數,我們必須增加一個額外的參數,表示不允許重復:
# 無放回抽取
> x3 <- sample(1:40, 6, replace=F)
> x3
[1] 3 40 20 19 28 11
c) 對一組數據進行亂序排序(隨機抽取,順序隨機)
你可以使用同樣的想法產生的任何載體的隨機子集,甚至不包含數字。 例如,選擇10個不同的美國各州隨機:
> sample(state.name, 10)
[1] "Delaware" "Vermont" "Rhode Island" "Tennessee"
[5] "Arizona" "Mississippi" "Virginia" "Alaska"
[9] "Georgia" "Louisiana"
> sample(state.name, 52)
Error in sample.int(length(x), size, replace, prob) :
cannot take a sample larger than the population when 'replace = FALSE'
> length(state.name)
[1] 50
> sample(state.name, 50)
[1] "Oregon" "Georgia" "Maine"
[4] "Idaho" "Alaska" "Tennessee"
[7] "Indiana" "Wyoming" "Montana"
[10] "Utah" "Florida" "North Carolina"
[13] "Nevada" "Virginia" "Pennsylvania"
[16] "North Dakota" "Wisconsin" "Alabama"
[19] "Mississippi" "New Hampshire" "Delaware"
[22] "Arizona" "Massachusetts" "Vermont"
[25] "Maryland" "Missouri" "Michigan"
[28] "Connecticut" "Colorado" "New York"
[31] "Oklahoma" "California" "Washington"
[34] "South Carolina" "West Virginia" "Hawaii"
[37] "Illinois" "Arkansas" "Kentucky"
[40] "Iowa" "Kansas" "Rhode Island"
[43] "New Mexico" "South Dakota" "New Jersey"
[46] "Nebraska" "Ohio" "Texas"
[49] "Louisiana" "Minnesota"
> sample(state.name, 50)
[1] "Arkansas" "California" "Texas"
[4] "Tennessee" "Montana" "Massachusetts"
[7] "North Dakota" "Oregon" "Delaware"
[10] "Hawaii" "South Dakota" "Connecticut"
[13] "South Carolina" "Kansas" "Washington"
[16] "West Virginia" "Georgia" "Maine"
[19] "Wyoming" "Illinois" "Nebraska"
[22] "Idaho" "Maryland" "Ohio"
[25] "Indiana" "Louisiana" "Mississippi"
[28] "Iowa" "New York" "Minnesota"
[31] "Rhode Island" "New Mexico" "New Jersey"
[34] "Pennsylvania" "Kentucky" "Utah"
[37] "Nevada" "Arizona" "Missouri"
[40] "North Carolina" "New Hampshire" "Wisconsin"
[43] "Alaska" "Vermont" "Alabama"
[46] "Florida" "Oklahoma" "Michigan"
[49] "Colorado" "Virginia"
>
d) 隨機從矩陣(數據框)中選取一部分對象
-----------------------------------------------------
(col.name=colnames(mtcars))
(row.name=rownames(mtcars))
#列名向量不返回抽樣
(sam.col.name=sample(col.name,10,replace=FALSE))
#行名向量不返回抽樣
(sam.row.name=sample(row.name,10,replace=FALSE))
B=mtcars[sam.row.name,sam.col.name]
B
-----------------------------------------------------
> (col.name=colnames(mtcars))
[1] "mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs" "am"
[10] "gear" "carb"
> (row.name=rownames(mtcars))
[1] "Mazda RX4" "Mazda RX4 Wag"
[3] "Datsun 710" "Hornet 4 Drive"
[5] "Hornet Sportabout" "Valiant"
[7] "Duster 360" "Merc 240D"
[9] "Merc 230" "Merc 280"
[11] "Merc 280C" "Merc 450SE"
[13] "Merc 450SL" "Merc 450SLC"
[15] "Cadillac Fleetwood" "Lincoln Continental"
[17] "Chrysler Imperial" "Fiat 128"
[19] "Honda Civic" "Toyota Corolla"
[21] "Toyota Corona" "Dodge Challenger"
[23] "AMC Javelin" "Camaro Z28"
[25] "Pontiac Firebird" "Fiat X1-9"
[27] "Porsche 914-2" "Lotus Europa"
[29] "Ford Pantera L" "Ferrari Dino"
[31] "Maserati Bora" "Volvo 142E"
> #列名向量不返回抽樣
> (sam.col.name=sample(col.name,10,replace=FALSE))
[1] "mpg" "qsec" "wt" "disp" "carb" "cyl" "vs" "am" "hp"
[10] "gear"
> #行名向量不返回抽樣
> (sam.row.name=sample(row.name,10,replace=FALSE))
[1] "Merc 230" "Pontiac Firebird" "Maserati Bora"
[4] "Hornet Sportabout" "Merc 450SE" "Merc 280"
[7] "Duster 360" "Merc 280C" "Dodge Challenger"
[10] "Camaro Z28"
> B=mtcars[sam.row.name,sam.col.name]
> B
mpg qsec wt disp carb cyl vs am hp gear
Merc 230 22.8 22.90 3.150 140.8 2 4 1 0 95 4
Pontiac Firebird 19.2 17.05 3.845 400.0 2 8 0 0 175 3
Maserati Bora 15.0 14.60 3.570 301.0 8 8 0 1 335 5
Hornet Sportabout 18.7 17.02 3.440 360.0 2 8 0 0 175 3
Merc 450SE 16.4 17.40 4.070 275.8 3 8 0 0 180 3
Merc 280 19.2 18.30 3.440 167.6 4 6 1 0 123 4
Duster 360 14.3 15.84 3.570 360.0 4 8 0 0 245 3
Merc 280C 17.8 18.90 3.440 167.6 4 6 1 0 123 4
Dodge Challenger 15.5 16.87 3.520 318.0 2 8 0 0 150 3
Camaro Z28 13.3 15.41 3.840 350.0 4 8 0 0 245 3
>
-----------------------------------------------------
e) 有放回的隨機抽樣示例
在e里隨機有放回地取N數,形成一個向量
e <- c(1:6)
N <- 100
s <- sample(x=e, size=N, replace=TRUE)
table(s)
> e <- c(1:6)
> N <- 100
> s <- sample(x=e, size=N, replace=TRUE)
> table(s)
s
1 2 3 4 5 6
16 19 14 21 19 11
f) 用sample函數實現隨機排序(抽樣)
sample函數第一個參數x表示被抽樣的向量,第二個參數size表示抽取的樣本個數,第三個參數replace表示是不是有放回的抽樣。
-----------------------------------------------------
Vec <- 1:10
# 無放回抽樣
sample(x = Vec, size = 5)
sample(x = Vec, size = 5)
# 只要第二個參數size的數值和第一個參數x的長度相同就可以實現隨機排序。
sample(x = Vec, size = 10)
sample(x = Vec, size = 10)
# 有放回抽樣
sample(x = Vec, size = 10, replace = TRUE)
sample(x = Vec, size = 10, replace = TRUE)
-----------------------------------------------------
> Vec <- 1:10
> # 無放回抽樣
> sample(x = Vec, size = 5)
[1] 9 3 1 5 4
> sample(x = Vec, size = 5)
[1] 3 1 4 9 10
> # 只要第二個參數size的數值和第一個參數x的長度相同就可以實現隨機排序。
> sample(x = Vec, size = 10)
[1] 1 6 2 10 7 8 3 9 5 4
> sample(x = Vec, size = 10)
[1] 6 9 10 7 1 5 8 3 4 2
> # 有放回抽樣
> sample(x = Vec, size = 10, replace = TRUE)
[1] 3 6 2 2 8 8 3 6 10 8
> sample(x = Vec, size = 10, replace = TRUE)
[1] 7 5 10 5 7 3 1 6 4 7
-----------------------------------------------------
g) rep()和replicate()批量取隨機數
問題:假設我想從符合正態分布的數據集中隨機抽取2個數據,然后排序,而這樣的數據我需要20對,你會怎么做?
# 方法1:rep()函數
rep(sort(sample(rnorm(n=100,mean = 0,sd = 1),2)),20)
> rep(sort(sample(rnorm(n=100,mean = 0,sd = 1),2)),20)
[1] -0.7433281 0.1923704 -0.7433281 0.1923704
[5] -0.7433281 0.1923704 -0.7433281 0.1923704
[9] -0.7433281 0.1923704 -0.7433281 0.1923704
[13] -0.7433281 0.1923704 -0.7433281 0.1923704
[17] -0.7433281 0.1923704 -0.7433281 0.1923704
[21] -0.7433281 0.1923704 -0.7433281 0.1923704
[25] -0.7433281 0.1923704 -0.7433281 0.1923704
[29] -0.7433281 0.1923704 -0.7433281 0.1923704
[33] -0.7433281 0.1923704 -0.7433281 0.1923704
[37] -0.7433281 0.1923704 -0.7433281 0.1923704
# 方法2:replicate()函數
replicate(n=20,expr=sort(sample(rnorm(n=100,mean = 0,sd = 1),2)))
> replicate(n=20,expr=sort(sample(rnorm(n=100,mean = 0,sd = 1),2)))
[,1] [,2] [,3] [,4]
[1,] 0.1214492 0.1682693 -1.4858967 -0.2435413
[2,] 0.5112295 0.6430365 0.1062271 0.1608918
[,5] [,6] [,7] [,8]
[1,] -0.7942942 -0.9092859 -0.38338123 -1.3953220
[2,] 0.6556288 1.0844620 -0.06132436 -0.3722372
[,9] [,10] [,11] [,12]
[1,] -0.6628561 -0.9390822 -1.772896 -1.886010
[2,] 0.6303991 1.2778667 -1.588206 -1.099241
[,13] [,14] [,15] [,16]
[1,] -1.1779682 -0.6320272 -0.7623491 -0.03874195
[2,] -0.5253614 0.1441117 0.4199327 0.70510390
[,17] [,18] [,19] [,20]
[1,] 0.1193869 -1.774997 -1.4304266 -0.7024574
[2,] 0.3463737 -1.338911 0.1309915 0.3383929
>
小結:rep()返回的是向量,replicate()返回的是矩陣。
下面列出兩個函數的用法:
rep():
rep(x, ...)
rep.int(x, times) # 每個元素重復次數
rep_len(x, length.out) # 生成向量長度
# 隨機數生成器
replicate(),replicate(n, expr, simplify = "array")
3. 設置隨機數種子
> set.seed(1234)
> x <- rnorm(10)
> x
[1] -1.2070657 0.2774292 1.0844412 -2.3456977 0.4291247
[6] 0.5060559 -0.5747400 -0.5466319 -0.5644520 -0.8900378
>
> # 設置了隨機數種子后,每次產生的分布數都是相同的
> set.seed(1234)
> x <- rnorm(10)
> x
[1] -1.2070657 0.2774292 1.0844412 -2.3456977 0.4291247
[6] 0.5060559 -0.5747400 -0.5466319 -0.5644520 -0.8900378
>
> # 移除了隨機數種子后,產生的隨機分布改變了
> rm(.Random.seed)
>
> x <- rnorm(10)
> x
[1] 0.6977298 0.6580453 0.5422976 0.1453601 -0.4882267
[6] 0.6037503 -0.8052218 0.3463612 -2.3368599 0.6643264
>