潘素娟，丁杰

(1.福建商學(xué)院 信息工程學(xué)院， 福建 福州 350012； 2.金融數(shù)學(xué)福建省高校重點(diǎn)實(shí)驗(yàn)室(莆田學(xué)院)，福建 莆田 351100； 3.廈門大學(xué) 經(jīng)濟(jì)學(xué)院， 福建 廈門 361005 )

0 引言

研究抽簽概率不僅能促進(jìn)統(tǒng)計和概率理論的發(fā)展，同時也能解決許多現(xiàn)實(shí)的問題，比如招標(biāo)或投標(biāo)項目、解決城市有限公共資源的供需矛盾、構(gòu)建和運(yùn)行馬拉松賽事信息系統(tǒng)、設(shè)計和優(yōu)化政府征收等.目前，國內(nèi)外已有許多學(xué)者對抽簽概率進(jìn)行了研究，并取得了較好的研究結(jié)果[1-6].在體育競技淘汰賽中，競賽對象的不同可能會直接影響到競賽的結(jié)果.為了解決競賽過程中出現(xiàn)的機(jī)遇性強(qiáng)和競賽結(jié)果的偶然性大等問題,通常采用抽簽來決定競賽對象，以此最大限度地保證競賽過程的公正性和競賽結(jié)果的合理性.目前，針對體育競技淘汰賽的抽簽方法主要有“1/4區(qū)公式控制抽簽法”“1/2區(qū)逐級分區(qū)抽簽法”“逐區(qū)雙分抽簽法”等方法.其中：“1/4區(qū)公式控制抽簽法”能夠較好地解決淘汰制抽簽中“機(jī)遇”與“控制”的矛盾，但該方法的抽簽過程較為繁瑣，且進(jìn)區(qū)概率存在失真性等問題[7].“1/2區(qū)逐級分區(qū)抽簽法”雖然能夠克服“1/4區(qū)公式控制抽簽法”所存在的缺陷，但該方法并未考慮受隨機(jī)因素影響的淘汰賽抽簽[8].“逐區(qū)雙分抽簽法”雖然能夠使得整個抽簽過程做到最大限度的隨機(jī)，但該方法并未考慮很多淘汰賽規(guī)則中所隱含的陷阱信息[9].目前為止，對淘汰賽中所隱含的陷阱還沒有較好的模擬方法可以實(shí)現(xiàn)其對陣概率的求解.

蒙特卡洛模擬算法[10]是以概率統(tǒng)計理論為基礎(chǔ)的一種計算方法，它可以隨機(jī)模擬各種變量間的動態(tài)關(guān)系，把每一種不確定性對結(jié)果的影響以概率分布的形式表示出來，因此該方法可以解決數(shù)值解的求解問題.目前，已有學(xué)者利用蒙特卡洛模擬方法研究了彩票的概率模型，并取得了較好的效果[11]，但利用該方法來研究抽簽概率的研究尚未見到報道.因此，基于上述難以求出對陣概率解析解的抽簽概率問題，本文以2017—2018賽季的歐冠淘汰賽抽簽概率為例，利用蒙特卡洛模擬方法求出進(jìn)入歐冠淘汰賽的各支球隊之間的對陣概率，并分別利用置信區(qū)間和分位點(diǎn)對模擬結(jié)果的可信度進(jìn)行了分析.

1 研究設(shè)計

1.1 研究對象

歐冠2017—2018賽季的淘汰賽一共有16支球隊.小組賽中以小組第1身份出線的球隊有：曼聯(lián)(A組，英超)、巴黎(B組，法甲)、羅馬(C組，意甲)、巴薩(D組，西甲)、利物浦(E組，英超)、曼城(F組，英超)、貝西克塔斯(G組，土超)和熱刺(H組，英超)；小組賽中以小組第2身份出線的球隊有：巴塞爾(A組，瑞超)、拜仁(B組，德甲)、切爾西(C組，英超)、尤文(D組，意甲)、塞維利亞(E組，西甲)、礦工(F組，烏超)、波爾圖(G組，葡超)和皇馬(H組，西甲).

1.2 抽簽規(guī)則

歐冠2017—2018賽季的淘汰賽采用抽簽的辦法確定競賽對象.若對所有球隊都不做任何約束條件，則每一支球隊都可能抽到任何一個競爭對象，即抽簽只是為了保證所有球隊都有一個同等的機(jī)遇條件.但比賽組織方對歐冠2017—2018賽季淘汰賽的抽簽制定了如下基本規(guī)則：

1)小組第1的球隊對陣小組第2的球隊.

2)同個小組的球隊之間回避，如曼聯(lián)不能和小組賽同組的巴塞爾對陣.

3)同國聯(lián)賽的球隊規(guī)避，如巴薩不能與同為西甲的皇馬對決.

1.3 模型的建立

建立模型的目標(biāo)是計算一個8×8的矩陣.矩陣內(nèi)的元素對應(yīng)的是兩支球隊之間的對陣概率，矩陣的行對應(yīng)的是小組第1出線的球隊(曼聯(lián)標(biāo)記為1，巴黎標(biāo)記為2，羅馬標(biāo)記為3，巴薩標(biāo)記為4，利物浦標(biāo)記為5，曼城標(biāo)記為6，貝西克塔斯標(biāo)記為7，熱刺標(biāo)記為8)，矩陣的列對應(yīng)的是小組第2出線的球隊(巴塞爾標(biāo)記為1，拜仁標(biāo)記為2，切爾西標(biāo)記為3，尤文標(biāo)記為4，塞維利亞標(biāo)記為5，礦工標(biāo)記為6，波爾圖標(biāo)記為7，皇馬標(biāo)記為8).

1.4 抽簽流程

抽簽的流程為：第1輪，先從小組第1的球隊中抽出1支球隊，然后從它可能對應(yīng)的(要考慮到規(guī)則2和規(guī)則3中的回避條款)小組第2中的對手中抽出1支球隊，然后這兩支球隊組成對陣雙方；第2輪，在剩下的球隊中再次進(jìn)行上一步操作.依據(jù)上述流程，只有完成第8輪抽簽才可完成對陣抽簽，但在實(shí)際抽簽過程中會存在如下陷阱：

1)雖然理論上要進(jìn)行8輪抽簽才能確定最后的對陣，但在實(shí)際中抽到第7輪就已經(jīng)能夠確定最后的對陣情況.其原因是， 7輪抽簽過后小組第1和小組第2的球隊都只剩下1支，所以肯定是他們之間進(jìn)行對決.

2)假設(shè)6輪過后，在剩下的球隊中小組第1的球隊是巴薩和貝西克塔斯，小組第2的球隊是礦工和皇馬.在這種情況下，經(jīng)過6輪抽簽即可確定對陣情況.因?yàn)樾〗M第1的巴薩隊不能對陣皇馬隊(同國聯(lián)賽回避)，所以巴薩隊只能對陣礦工隊，由此可以推出貝西克塔斯隊對陣皇馬隊.

3)在陷阱2)中，小組第2的皇馬隊不能對陣巴薩隊(同國聯(lián)賽回避)，所以皇馬隊只能對陣貝西克塔斯隊，由此可以推出巴薩隊對陣礦工隊.

陷阱2)和3)在規(guī)則方面較為接近，但在有些情形中兩者并不一致.例如：假設(shè)前2輪抽簽結(jié)束后，對陣結(jié)果是巴黎隊對陣皇馬隊，貝西克塔斯隊對陣拜仁隊.在這種情況下，由于存在同組規(guī)避和同國聯(lián)賽規(guī)避的條款，因此小組第2的切爾西隊對陣的球隊只有巴黎、巴薩和貝西克塔斯，然而在巴黎隊和貝西克塔斯隊都已經(jīng)被抽走的情況下，與切爾西隊對陣的只能是巴薩隊.上述情形所對應(yīng)的就是陷阱3)，而非陷阱2).根據(jù)歐冠規(guī)則，同一個聯(lián)賽原則上最多的歐冠名額是“3+1”個，其中3個球隊直接進(jìn)入歐冠小組賽，另1個球隊參加附加賽(如果附加賽獲勝就可以參加小組賽).另外，如果該聯(lián)賽有1支其他的球隊在上一年獲得了歐聯(lián)杯的冠軍，那么該支球隊也可以直接獲得參加歐冠小組賽的資格.所以從理論上說，某個聯(lián)賽參加歐冠小組賽的球隊最多可以是5支.如果這5支球隊都能小組出線，那么該聯(lián)賽有5支球隊能夠參加歐冠淘汰賽.如果其中4支球隊為小組第1，另1支球隊為小組第2；或者其中4支球隊為小組第2，另1支球隊小組第1：那么和其他4支球隊不在同一個組的那支球隊不僅需要回避4支同國聯(lián)賽的球隊，還要回避同小組的球隊.所以該球隊需要回避的球隊有5支，能夠選擇對陣的球隊只有3支.這種最極端的例子出現(xiàn)在歐冠2017—2018賽季的淘汰賽抽簽中，即切爾西隊就是那支需要回避5支球隊的球隊.在這種最極端的情形下，經(jīng)過兩輪抽簽后就可決定球隊的對陣情況，即至少經(jīng)過2輪抽簽過后才可能會出現(xiàn)陷阱2)或陷阱3)的現(xiàn)象.

2 模擬過程

2.1 利用蒙特卡洛模擬計算對陣概率數(shù)值解的步驟

第1步 生成1個8×8且元素全部為0的矩陣.

第2步 根據(jù)規(guī)則先抽出一組對陣結(jié)果(即16支球隊的8場對陣結(jié)果).如果2支球隊之間的抽簽結(jié)果為對陣，則在矩陣中對與該結(jié)果相對應(yīng)的元素加1； 如果2支球隊之間抽簽的結(jié)果不為對陣，則矩陣中相應(yīng)的元素不變.

第3步 重復(fù)步驟2，直到完成100 000次蒙特卡羅模擬.

第4步 將步驟3最終得到的矩陣中每個元素的最后累積值除以蒙特卡洛模擬的次數(shù)(100 000次)，即可得到對陣概率的數(shù)值解.

2.2 編程

2.1模擬步驟中的第2步的編程方法為：①根據(jù)規(guī)則任意抽取前兩輪(因在最極端的情形下，前兩輪也不會掉進(jìn)陷阱2)或者陷阱3)).②從第3輪開始，首先觀察每個小組第1的球隊會不會出現(xiàn)只剩下1個可以選擇的對手的情形.如果有，直接抽取出來；如果沒有，就進(jìn)入下一步.然后觀察每個小組第2的球隊會不會出現(xiàn)只剩下一個可以選擇的對手的情形.如果有，直接抽取出來；如果沒有，就可以和前兩輪一樣根據(jù)規(guī)則(即在同組回避和同國聯(lián)賽規(guī)避的條件下)進(jìn)行抽取.③第4輪到第7輪的抽簽思路與第3輪相同.在所有的情形下，實(shí)際上不需要進(jìn)行第8輪抽簽，因?yàn)榻?jīng)過第7輪抽簽后就可確定最終的對陣結(jié)果.經(jīng)編程(見附錄)得到的對陣概率如表1所示.

表1 各球隊之間的對陣概率

3 可信度分析

利用蒙特卡羅模擬方法研究數(shù)值解問題時，若輸入模式中的隨機(jī)數(shù)并不是真正的隨機(jī)數(shù)，則模擬及預(yù)測結(jié)果就會產(chǎn)生錯誤.為此，本文利用區(qū)間估計的方法(置信區(qū)間和分位點(diǎn))對上述模擬得到的對陣概率結(jié)果進(jìn)行可信度分析.

1)置信區(qū)間法.利用模型對對陣相關(guān)的64個參數(shù)進(jìn)行100 000次蒙特卡洛模擬后發(fā)現(xiàn)，每次模擬都服從I.I.D.假設(shè)[12].為了驗(yàn)證可信度，對上述對陣概率的置信區(qū)間進(jìn)行100次實(shí)驗(yàn)，每次實(shí)驗(yàn)均進(jìn)行100 000次蒙特卡洛模擬.根據(jù)實(shí)驗(yàn)結(jié)果計算出的對陣概率的95%水平下的置信區(qū)間和99%水平下的置信區(qū)間的結(jié)果如表2和表3所示.由表2和表3可以看出，表中的置信區(qū)間都比較狹窄，說明本文提出的模型具有較高的可信度.

表2 對陣概率在95%水平下的置信區(qū)間

表3 對陣概率在99%水平下的置信區(qū)間

2)分位點(diǎn)法.以下通過分析兩組分位點(diǎn)的結(jié)果來驗(yàn)證模擬結(jié)果的可信度.第1組為2.5百分位點(diǎn)到97.5百分位點(diǎn)，如表4所示；第2組為0.5百分位點(diǎn)到99.5百分位點(diǎn)，如表5所示.由表4和表5可以看出，表中的區(qū)間都比較狹窄，該結(jié)果再次說明本文提出的模擬方法具有較高的可信度.

表4 對陣概率在2.5到97.5的百分位點(diǎn)

表5 對陣概率在0.5到99.5的百分位點(diǎn)

4 結(jié)論

本文通過對歐冠2017—2018賽季淘汰賽的抽簽規(guī)則和抽簽流程進(jìn)行分析，利用蒙特卡洛模擬算法建立了一種新型的抽簽概率模型，并給出了對陣概率的數(shù)值解.利用置信區(qū)間和分位點(diǎn)對模型計算所得的對陣概率進(jìn)行可信度分析表明，本文的模擬方法具有較好的可信度.模擬結(jié)果顯示，巴薩與切爾西對陣的概率約為40%(該結(jié)果與實(shí)際比賽結(jié)果相符)，由此再次表明本文模擬方法具有可靠性.另外，由于本研究已經(jīng)考慮了最極端的情形，所以對于其他歐冠賽季，在不改比賽規(guī)則的情況下只要改變初始數(shù)據(jù)集就可以得到對陣概率.本文方法也可以為金融市場的風(fēng)險管理、項目的招標(biāo)與投標(biāo)等問題的抽簽提供借鑒.

附錄：

R代碼：

library(lubridate)

set.seed(1234)

sink(file="C:/Users/Ding/Desktop/UEFA Champions League.txt"，append=T)

n_1 <- 100000 # 100000 times Monte Carlos Simulations each experiment

n_2 <- 100 # 100 experiments

now()

result_all <- array(NA，dim=c(8，8，n_2))

for (x in 1:n_2) {

result_cum <- matrix(0，nrow=8，ncol=8)

FUN_1 <- function(a，b) {

result[a，1] <<- set_1[b]

result[a，2] <<- set_2_2[[set_1[b]]]

set_1 <<- setdiff(set_1，result[a，1])

for (j in 1∶8) {

set_1_1[[j]] <<- setdiff(set_1_1[[j]]，result[a，1])

}

set_2 <<- setdiff(set_2，result[a，2])

for (j in 1∶8) {

set_2_2[[j]] <<- setdiff(set_2_2[[j]]，result[a，2])

}

}

FUN_2 <- function(a，b) {

result[a，1] <<- set_1_1[[set_2[b]]]

result[a，2] <<- set_2[b]

set_1 <<- setdiff(set_1，result[a，1])

for (j in 1∶8) {

set_1_1[[j]] <<- setdiff(set_1_1[[j]]，result[a，1])

}

set_2 <<- setdiff(set_2，result[a，2])

for (j in 1∶8) {

set_2_2[[j]] <<- setdiff(set_2_2[[j]]，result[a，2])

}

}

for (k in 1∶n_1) {

result <- matrix(NA，nrow=8，ncol=2)

set_1 <- 1∶8

set_2 <- 1∶8

set_1_1 <- list(setdiff(1∶8，1)，setdiff(1∶8，2)，setdiff(1∶8，c(1，3，5，6，8))，setdiff(1∶8，c(3，4))，setdiff(1∶8，c(4，5))，setdiff(1∶8，6)，setdiff(1∶8，7)，setdiff(1∶8，c(4，8)))

set_2_2 <- list(setdiff(1∶8，c(1，3))，setdiff(1∶8，2)，setdiff(1∶8，c(3，4))，setdiff(1∶8，c(4，5，8))，setdiff(1∶8，c(3，5))，setdiff(1∶8，c(3，6))，setdiff(1∶8，7)，setdiff(1∶8，c(3，8)))

for (i in 1∶2) {

result[i，1] <- sample(set_1，1)

result[i，2] <- sample(set_2_2[[result[i，1]]]，1)

set_1 <- setdiff(set_1，result[i，1])

for (j in 1∶8) {

set_1_1[[j]] <- setdiff(set_1_1[[j]]，result[i，1])

}

set_2 <- setdiff(set_2，result[i，2])

for (j in 1∶8) {

set_2_2[[j]] <- setdiff(set_2_2[[j]]，result[i，2])

}

}

if (length(set_2_2[[set_1[1]]])==1) {

FUN_1(3，1)

} else if (length(set_2_2[[set_1[2]]])==1) {

FUN_1(3，2)

} else if (length(set_2_2[[set_1[3]]])==1) {

FUN_1(3，3)

} else if (length(set_2_2[[set_1[4]]])==1) {

FUN_1(3，4)

} else if (length(set_2_2[[set_1[5]]])==1) {

FUN_1(3，5)

} else if (length(set_2_2[[set_1[6]]])==1) {

FUN_1(3，6)

} else if (length(set_1_1[[set_2[1]]])==1) {

FUN_2(3，1)

} else if (length(set_1_1[[set_2[2]]])==1) {

FUN_2(3，2)

} else if (length(set_1_1[[set_2[3]]])==1) {

FUN_2(3，3)

} else if (length(set_1_1[[set_2[4]]])==1) {

FUN_2(3，4)

} else if (length(set_1_1[[set_2[5]]])==1) {

FUN_2(3，5)

} else if (length(set_1_1[[set_2[6]]])==1) {

FUN_2(3，6)

} else {

result[3，1] <- sample(set_1，1)

result[3，2] <- sample(set_2_2[[result[3，1]]]，1)

set_1 <- setdiff(set_1，result[3，1])

for (j in 1∶8) {

set_1_1[[j]] <- setdiff(set_1_1[[j]]，result[3，1])

}

set_2 <- setdiff(set_2，result[3，2])

for (j in 1∶8) {

set_2_2[[j]] <- setdiff(set_2_2[[j]]，result[3，2])

}

}

if (length(set_2_2[[set_1[1]]])==1) {

FUN_1(4，1)

} else if (length(set_2_2[[set_1[2]]])==1) {

FUN_1(4，2)

} else if (length(set_2_2[[set_1[3]]])==1) {

FUN_1(4，3)

} else if (length(set_2_2[[set_1[4]]])==1) {

FUN_1(4，4)

} else if (length(set_2_2[[set_1[5]]])==1) {

FUN_1(4，5)

} else if (length(set_1_1[[set_2[1]]])==1) {

FUN_2(4，1)

} else if (length(set_1_1[[set_2[2]]])==1) {

FUN_2(4，2)

} else if (length(set_1_1[[set_2[3]]])==1) {

FUN_2(4，3)

} else if (length(set_1_1[[set_2[4]]])==1) {

FUN_2(4，4)

} else if (length(set_1_1[[set_2[5]]])==1) {

FUN_2(4，5)

} else {

result[4，1] <- sample(set_1，1)

result[4，2] <- sample(set_2_2[[result[4，1]]]，1)

set_1 <- setdiff(set_1，result[4，1])

for (j in 1∶8) {

set_1_1[[j]] <- setdiff(set_1_1[[j]]，result[4，1])

}

set_2 <- setdiff(set_2，result[4，2])

for (j in 1∶8) {

set_2_2[[j]] <- setdiff(set_2_2[[j]]，result[4，2])

}

}

if (length(set_2_2[[set_1[1]]])==1) {

FUN_1(5，1)

} else if (length(set_2_2[[set_1[2]]])==1) {

FUN_1(5，2)

} else if (length(set_2_2[[set_1[3]]])==1) {

FUN_1(5，3)

} else if (length(set_2_2[[set_1[4]]])==1) {

FUN_1(5，4)

} else if (length(set_1_1[[set_2[1]]])==1) {

FUN_2(5，1)

} else if (length(set_1_1[[set_2[2]]])==1) {

FUN_2(5，2)

} else if (length(set_1_1[[set_2[3]]])==1) {

FUN_2(5，3)

} else if (length(set_1_1[[set_2[4]]])==1) {

FUN_2(5，4)

} else {

result[5，1] <- sample(set_1，1)

result[5，2] <- sample(set_2_2[[result[5，1]]]，1)

set_1 <- setdiff(set_1，result[5，1])

for (j in 1∶8) {

set_1_1[[j]] <- setdiff(set_1_1[[j]]，result[5，1])

}

set_2 <- setdiff(set_2，result[5，2])

for (j in 1∶8) {

set_2_2[[j]] <- setdiff(set_2_2[[j]]，result[5，2])

}

}

if (length(set_2_2[[set_1[1]]])==1) {

FUN_1(6，1)

} else if (length(set_2_2[[set_1[2]]])==1) {

FUN_1(6，2)

} else if (length(set_2_2[[set_1[3]]])==1) {

FUN_1(6，3)

} else if (length(set_1_1[[set_2[1]]])==1) {

FUN_2(6，1)

} else if (length(set_1_1[[set_2[2]]])==1) {

FUN_2(6，2)

} else if (length(set_1_1[[set_2[3]]])==1) {

FUN_2(6，3)

} else {

result[6，1] <- sample(set_1，1)

result[6，2] <- sample(set_2_2[[result[6，1]]]，1)

set_1 <- setdiff(set_1，result[6，1])

for (j in 1∶8) {

set_1_1[[j]] <- setdiff(set_1_1[[j]]，result[6，1])

}

set_2 <- setdiff(set_2，result[6，2])

for (j in 1∶8) {

set_2_2[[j]] <- setdiff(set_2_2[[j]]，result[6，2])

}

}

if (length(set_2_2[[set_1[1]]])==1) {

FUN_1(7，1)

} else if (length(set_2_2[[set_1[2]]])==1) {

FUN_1(7，2)

} else if (length(set_2_2[[set_1[3]]])==1) {

FUN_1(7，3)

} else if (length(set_1_1[[set_2[1]]])==1) {

FUN_2(7，1)

} else if (length(set_1_1[[set_2[2]]])==1) {

FUN_2(7，2)

} else if (length(set_1_1[[set_2[3]]])==1) {

FUN_2(7，3)

} else {

result[7，1] <- sample(set_1，1)

result[7，2] <- sample(set_2_2[[result[7，1]]]，1)

set_1 <- setdiff(set_1，result[7，1])

for (j in 1∶8) {

set_1_1[[j]] <- setdiff(set_1_1[[j]]，result[7，1])

}

set_2 <- setdiff(set_2，result[7，2])

for (j in 1∶8) {

set_2_2[[j]] <- setdiff(set_2_2[[j]]，result[7，2])

}

}

result[8，1] <- setdiff(1∶8，result[1∶7，1])

result[8，2] <- setdiff(1∶8，result[1∶7，2])

for (i in 1∶8) {

result_cum[result[i，1]，result[i，2]] <- result_cum[result[i，1]，result[i，2]]+1

}

#print(k)

}

#print(x)

print(now())

result_all[，，x] <- result_cum/n_1

}

result_mean <- matrix(NA，nrow=8，ncol=8)

for (i in 1∶8) {

for (j in 1∶8) {

result_mean[i，j] <- mean(result_all[i，j，])

}

}

result_sd <- matrix(NA，nrow=8，ncol=8)

for (i in 1∶8) {

for (j in 1∶8) {

result_sd[i，j] <- sd(result_all[i，j，])

}

}

result_all[，，1]

#### Confidence Interval

ci_95 <- data.frame(cbind(rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)))

#ci denotes confidence interval

colnames(ci_95) <- as.character(1∶8)

# 95% confidence interval

for (i in 1∶8) {

for (j in 1∶8) {

ci_95[i，j] <- paste0("("，round((result_mean+qnorm(.025)*result_sd/sqrt(n_2))[i，j]*100，4)，"，"，round((result_mean+qnorm(.975)*result_sd/sqrt(n_2))[i，j]*100，4)，")")

}

}

ci_95 # in percentage

ci_99 <- data.frame(cbind(rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)))

colnames(ci_99) <- as.character(1∶8)

# 99% confidence interval

for (i in 1∶8) {

for (j in 1∶8) {

ci_99[i，j] <- paste0("("，round((result_mean+qnorm(.005)*result_sd/sqrt(n_2))[i，j]*100，4)，"，"，round((result_mean+qnorm(.995)*result_sd/sqrt(n_2))[i，j]*100，4)，")")

}

}

ci_99 # in percentage

#### Quantile

# From .025th quantile to .975th quantile

quantile_1 <- data.frame(cbind(rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)))

colnames(quantile_1) <- as.character(1∶8)

for (i in 1∶8) {

for (j in 1∶8) {

quantile_1[i，j] <- paste0("("，round(quantile(result_all[i，j，]，.025)*100，4)，"，"，round(quantile(result_all[i，j，]，.975)*100，4)，")")

}

}

quantile_1 # in percentage

#### From .005th quantile to .995th quantile

quantile_2 <- data.frame(cbind(rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)，rep(NA，8)))

colnames(quantile_2) <- as.character(1∶8)

for (i in 1∶8) {

for (j in 1∶8) {

quantile_2[i，j] <- paste0("("，round(quantile(result_all[i，j，]，.005)*100，4)，"，"，round(quantile(result_all[i，j，]，.995)*100，4)，")")

}

}

quantile_2 # in percentage

sink()

rm(list=ls())

1/4+1/5.5+1/6+1/8.5+1/9+1/15+1/15+1/17+1/28+1/34

注：

1) set_1表示在第1小組中可抽的球隊，set_2表示在第2小組中可抽的球隊.

2) set_1_1表示一個長度為8的列表，其中的每一個元素都是一個向量.元素i表示第2小組中的第i支球隊目前能夠?qū)﹃嚨?小組的球隊編號.

3) set_2_2表示一個長度為8的列表，其中的每一個元素都是一個向量.元素i表示第1小組中的第i支球隊目前能夠?qū)﹃嚨?小組的球隊編號.

4)為了提高代碼的效率，可以去掉顯式循環(huán)，使用并行或者RCPP來提高R語言的循環(huán)速度.若需進(jìn)一步提高模擬的精確度，可增加n_1和n_2的數(shù)值，由此則必然存在某個正數(shù)，只要n_1和n_2大于該數(shù)字就可滿足設(shè)定的精確性要求.