國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
3
第二章
"rBeta2009"套件之介紹
本章將介紹"rBeta2009"套件中，用來生成貝他分配所使用的方針，以及生成 狄式分配亂數所使用的方針與演算法。
第一節 生成貝他分配之演算法
對於貝他分配 𝐵𝑒𝑡𝑎(𝛼, 𝛽)，Hung et al. (2009)提出三點方針，根據貝他分配的 參數大小來選擇生成貝他分配亂數最快速的演算法。以下是其方針：
方針一：當 𝛼, 𝛽 < 1 ，若 𝛼 + 𝛽 > 1.2 ，則使用演算法 Kennedy’s MK [20]；
其餘則使用演算法 B00 [30]；
方針二：當 α < 1 < 𝛽 或 α > 1 > β ，則使用演算法 B01 [30]；
方針三：當 𝛼, 𝛽 > 1 ，若其中一個參數很接近 1 ，且另一個參數很大(例 如 > 4 )，則使用演算法 B4PE [32]；其餘則使用演算法 BPRS [39]。
"rBeta2009"套件使用以上三個方針來產生貝他分配亂數。根據以上方針，就 能利用演算法 B00、B01、B4PE 及 BPRS 的特性，以最快的時間生成各種參數 的貝他分配亂數。值得一提的是，Kennedy’s MK 演算法為一"近似"貝他分配的 隨機搜尋方法，其產生的亂數有較大的不確定性(或誤差)，所以在"rBeta2009"套 件中，此演算法被 B00 所取代。
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
4
第二節 生成狄氏分配之演算法
狄式分配可以由貝他分配經過轉換而得到。若 𝑋1, … , 𝑋𝑘 為互相獨立的隨機 變 數 服 從 貝 他 分 配 𝑋𝑖~𝐵𝑒𝑡𝑎(𝛼𝑖, ∑𝑘𝑗=𝑖+1𝑎𝑗), 𝑖 = 1, … , 𝑘 ， 令 𝑌1 = 𝑋1, 𝑌𝑖 = 𝑋𝑖∏𝑖−1𝑗=1(1 − 𝑋𝑗), 𝑖 = 2, … , 𝑘 ， 則 (𝑌1, … , 𝑌𝑘) 服 從 所 謂 的 狄 式 分 配 ， 並 記 為 𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼1, … , 𝛼𝑘; 𝛼𝑘+1)。為了使此轉換過程速度更快，Hung et al. (2011) 提出 了以下三個方針：
方針一：根據參數 𝛼𝑖 選擇最快的演算法來生成貝他分配亂數。(例如，使用
第一節介紹的方針來選擇演算法。)
方針二：將狄氏分配的參數重新排序。其想法是將原本的參數 𝛼1, … , 𝛼𝑘+1 由 小到大排序，稱為 𝛼(1), … , 𝛼(𝑘+1) ，並令 𝑚 (≤ 𝑘 + 1)為參數中小於 1 的個 數。
若 𝛼(𝑘+1) < 1 且 𝑘 + 1 是奇數，則生成
𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼(1), 𝛼(3), … , 𝛼(𝑘−1), 𝛼(𝑘+1), 𝛼(𝑘), 𝛼(𝑘−2), 𝛼(𝑘−4), … , 𝛼(4), 𝛼(2)) 隨機向量；
若𝛼(𝑘+1) < 1 且 𝑘 + 1 是偶數，則生成
𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼(1), 𝛼(2), 𝛼(4), … , 𝛼(𝑘−1), 𝛼(𝑘+1), 𝛼(𝑘−2), 𝛼(𝑘−4), … , 𝛼(5), 𝛼(3)) 隨機向量；
若 1 ≤ 𝛼(𝑘+1) ≤ 3 且 𝛼(1) > 1 ，則生成
𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼(𝑘+1), 𝛼(𝑘−1), 𝛼(𝑘−2), … , 𝛼(1), 𝛼(𝑘)) 隨機向量；
若 𝛼(𝑘+1) > 3 且 𝛼(1)≥ 1 ，則生成
𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼(𝑘+1), 𝛼(𝑘), 𝛼(𝑘−2), 𝛼(𝑘−3)… , 𝛼(1), 𝛼(𝑘−1)) 隨機向量；
若 𝛼(𝑘+1) ≥ 1, 𝛼(1) < 1, 𝑚 >𝑘+22 且 𝛼(𝑚)≤ 0.5 ，則生成
𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼(1), 𝛼(2), … , 𝛼(𝑚−1), 𝛼(𝑚+1), 𝛼(𝑚+2), … , 𝛼(𝑘+1), 𝛼(𝑚))
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
5
隨機向量；
其餘則生成
𝐷𝑖𝑟𝑖𝑐ℎ𝑙𝑒𝑡(𝛼(𝑘+1), 𝛼(𝑘−1), 𝛼(𝑘−2), … , 𝛼(1), 𝛼(𝑘)) 隨機向量。
方 針 三 ： 減 少 乘 法 的 運 算 次 數 。 以 𝑌𝑖 = 𝑋𝑖(1 − ∑𝑖−1𝑗=1𝑌𝑗) 代 替 原 本 的 𝑌𝑖 = 𝑋𝑖∏𝑖−1𝑗=1(1 − 𝑋𝑗)，減少乘法的運算次數，加快運算速度。
套件"rBeta2009"使用以上三個方針來生成狄式分配亂數；並且由 Hung et al.
(2011) 將以上整理成演算法，稱為 BETA-M。
BETA-M 演算法
步 驟 一 ： 根 據 上 述 的 方 針 二 重 新 排 序 狄 式 分 配 的 參 數 𝛼1, … , 𝛼𝑘+1， 定 為 𝛼𝑠1, … , 𝛼𝑠𝑘+1 ；
步驟二：以方針一產生獨立的貝他分配 𝑋𝑠𝑖~𝐵𝑒𝑡𝑎 (𝛼𝑠𝑖, ∑𝑘+1𝑗=𝑖+1𝛼𝑠𝑗) , 𝑖 = 1, … , 𝑘 ； 步驟三：令 𝑌𝑠1 = 𝑋𝑠1 及 𝑌𝑠𝑖 = 𝑋𝑠𝑖(1 − ∑𝑖−1𝑗=1𝑌𝑠𝑗) , 𝑖 = 2, … , 𝑘 ；
步驟四：根據步驟一的排序將 𝑌𝑠1, … , 𝑌𝑠𝑘 排回 𝑌1, … , 𝑌𝑘 並輸出。
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
6
第三節 "rBeta2009"之指令執行
本節將介紹如何使用 R 軟體套件"rBeta2009"中的指令 rbeta()及 rdirichlet()。
程式 rbeta(n, shape1, shape2)是以本章第一節所述之演算法，生成 𝑛 個貝他分配 𝐵𝑒𝑡𝑎(shape1, shape2) 亂 數 。 其 中 n 為 生 成 的 亂 數 個 數 ， 若 長 度 大 於 1 (length(n)>1)，則會以 n 的維度個數當作生成的亂數個數；而 shape1 及 shape2 為 正形狀參數(shape parameter)。要注意的是，R 軟體的預設套件"stats"中有同名的 程式，但只需讀取"rBeta2009"套件之後，在套件"stats"中的同名程式就能被取代。
以下會舉例介紹如何使用"rBeta2009"套件中的 rbeta()。首先，將"rBeta2009"套件 讀入 R 軟體並隨機設一個種子(set a random seed)：
R> library("rBeta2009") R> set.seed(63522)
接著我們生成 10 個貝塔分配 𝐵𝑒𝑡𝑎(1.5, 0.3) 亂數：
R>rbeta(10, 1.5, 0.3)
[1] 0.7461113 0.8843041 0.9960109 0.9968034 0.2236435 0.5429963 0.5130998 [8] 0.9999988 0.9817266 0.7213977
生成 10 個貝塔分配 𝐵𝑒𝑡𝑎(4.2, 13) 亂數：
R > rbeta(10, 4.3, 13)
[1] 0.1456586 0.2540797 0.4101915 0.2163145 0.1844294 0.4010572 0.1043971 [8] 0.2315380 0.2456694 0.2543450
而程式 rdirichlet(n, shape)是以本章第二節之演算法生成 𝑛 個 𝑘 維度的狄 氏分配亂數，其中 shape 為一個長度為 𝑘 + 1 正的形狀參數向量(𝑘 ≥ 2)，而 n 為生成的亂數個數，若長度大於 1 (length(n)>1)，則會以 n 的長度當作生成的亂
‧
R>rdirichlet(10, c(1.5, 0.5, 3, 1.7))
[,1] [,2] [,3] [,4]
R>rdirichlet(10, c(-1.5, 0.5, 3, 1.7))
錯誤在 rdirichlet(10, c(-1.5, 0.5, 3, 1.7)) : Shape parameters should be all positive
而若是 shape 向量長度為 2 時，會出現警告，但程式仍可以執行：
R> rdirichlet(10, c(3, 0.5))
[1] 0.9998994 0.9994591 0.8965159 0.9985557 0.8156766 0.9676691 0.9319676 [8] 0.9593898 0.7713237 0.7098055
警告訊息：
In rdirichlet(10, c(3, 0.5)) :
'shape' is of length 2: reduced to beta variates
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
8