基于蒙特卡洛模擬算法的歐冠淘汰賽抽簽概率的研究

潘素娟，丁杰

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

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

1 研究設(shè)計(jì)

1.1 研究對(duì)象

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

1.2 抽簽規(guī)則

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

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

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

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

1.3 模型的建立

建立模型的目標(biāo)是計(jì)算一個(gè)8×8的矩陣.矩陣內(nèi)的元素對(duì)應(yīng)的是兩支球隊(duì)之間的對(duì)陣概率，矩陣的行對(duì)應(yīng)的是小組第1出線的球隊(duì)(曼聯(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)，矩陣的列對(duì)應(yīng)的是小組第2出線的球隊(duì)(巴塞爾標(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的球隊(duì)中抽出1支球隊(duì)，然后從它可能對(duì)應(yīng)的(要考慮到規(guī)則2和規(guī)則3中的回避條款)小組第2中的對(duì)手中抽出1支球隊(duì)，然后這兩支球隊(duì)組成對(duì)陣雙方；第2輪，在剩下的球隊(duì)中再次進(jìn)行上一步操作.依據(jù)上述流程，只有完成第8輪抽簽才可完成對(duì)陣抽簽，但在實(shí)際抽簽過(guò)程中會(huì)存在如下陷阱：

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

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

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

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

2 模擬過(guò)程

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

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

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

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

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

2.2 編程

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

表1 各球隊(duì)之間的對(duì)陣概率

3 可信度分析

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

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

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

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

2)分位點(diǎn)法.以下通過(guò)分析兩組分位點(diǎn)的結(jié)果來(lái)驗(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é)果再次說(shuō)明本文提出的模擬方法具有較高的可信度.

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

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

4 結(jié)論

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

附錄：

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小組中可抽的球隊(duì)，set_2表示在第2小組中可抽的球隊(duì).

2) set_1_1表示一個(gè)長(zhǎng)度為8的列表，其中的每一個(gè)元素都是一個(gè)向量.元素i表示第2小組中的第i支球隊(duì)目前能夠?qū)﹃嚨?小組的球隊(duì)編號(hào).

3) set_2_2表示一個(gè)長(zhǎng)度為8的列表，其中的每一個(gè)元素都是一個(gè)向量.元素i表示第1小組中的第i支球隊(duì)目前能夠?qū)﹃嚨?小組的球隊(duì)編號(hào).

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