張 嵐， 錢夕元

（華東理工大學(xué)，理學(xué)院，上海 200237）

排隊(duì)論在社會(huì)各個(gè)領(lǐng)域都應(yīng)用廣泛，包括銀行排隊(duì)叫號(hào)、機(jī)場(chǎng)排隊(duì)取行李、電話排隊(duì)服務(wù)、交通運(yùn)輸?shù)?。但在?shí)際應(yīng)用中，在衡量一個(gè)排隊(duì)模型優(yōu)劣時(shí)，常需基于已知的到達(dá)率和服務(wù)率計(jì)算排隊(duì)模型的幾個(gè)重要指標(biāo)：系統(tǒng)平均隊(duì)長(zhǎng)、排隊(duì)等候平均隊(duì)長(zhǎng)、平均逗留時(shí)間、平均排隊(duì)時(shí)間，而通常排隊(duì)數(shù)據(jù)無法直接觀測(cè)到這兩個(gè)參數(shù)。在估計(jì)到達(dá)率和服務(wù)率，且是否使得排隊(duì)系統(tǒng)穩(wěn)定時(shí)，貝葉斯推斷和估計(jì)在排隊(duì)論中有很大的優(yōu)勢(shì)。

貝葉斯早在1763 年就提出貝葉斯公式，并給出貝葉斯估計(jì)的思想，根據(jù)每一類的樣本估計(jì)每一類的類條件概率密度，把參數(shù)θ看成未知的隨機(jī)變量，通過對(duì)樣本觀察再求貝葉斯估計(jì)。在貝葉斯估計(jì)時(shí)，首先需要選擇先驗(yàn)分布，再確定似然函數(shù)，確定參數(shù)的后驗(yàn)分布，選擇損失函數(shù)，最后估計(jì)參數(shù)。貝葉斯推斷在排隊(duì)系統(tǒng)中的優(yōu)勢(shì)：a. 系統(tǒng)穩(wěn)定性的不確定性可以容易被量化；b. 在參數(shù)空間中的限制容易被處理；c. 預(yù)期的參數(shù)可以直接得出；d. 貝葉斯決策技巧可以直接用于計(jì)算排隊(duì)論中的最優(yōu)決策。對(duì)于M/M/1，除了上述Bagchi 等，Armero 等[1-2]和Choudhury 等[3]研究了這個(gè)排隊(duì)系統(tǒng)下常用衡量效果測(cè)度的分布，如系統(tǒng)和排隊(duì)的顧客數(shù)、等待時(shí)間以及忙碌期、空閑期的時(shí)間長(zhǎng)度；并結(jié)合泊松間隔到達(dá)時(shí)間和指數(shù)分布服務(wù)時(shí)間進(jìn)行貝葉斯推斷，并給出聯(lián)合后驗(yàn)分布。對(duì)于非馬爾科夫系統(tǒng)：關(guān)于G/M/1，G/M/c系統(tǒng)研究較少，其中Wiper[4]分析了Er/M/1 和Er/M/c，研究了間隔到達(dá)時(shí)間服從愛爾蘭分布的排隊(duì)系統(tǒng)，他用蒙特卡洛模擬方法估計(jì)隊(duì)列長(zhǎng)度和等待時(shí)間的分布，且闡述了采用獨(dú)立先驗(yàn)時(shí)許多隊(duì)列的測(cè)度不存在。而對(duì)于非馬爾科夫服務(wù)時(shí)間，Insua 等[5-6]對(duì)M/Er/1 系統(tǒng)強(qiáng)度估計(jì)的基礎(chǔ)上進(jìn)一步用蒙特卡洛方法給出了M/G/1 的貝葉斯推斷和估計(jì)，并使得隊(duì)列達(dá)到均衡狀態(tài)。Ausín 等[7]通過Coxian 分布模型給出GI/G/1 系統(tǒng)的貝葉斯推斷，同時(shí)指出他們的方法可以用于估計(jì)瞬時(shí)隊(duì)伍大小以及等待時(shí)間的分布。

本文運(yùn)用近似貝葉斯的方法，結(jié)合排隊(duì)論QDC算法的思想對(duì)四個(gè)經(jīng)典排隊(duì)模型進(jìn)行仿真模擬，給出參數(shù)估計(jì)與估計(jì)效果；并結(jié)合實(shí)際例子（銀行的面對(duì)面排隊(duì)換匯數(shù)據(jù)）估計(jì)參數(shù)。

1 模 型

一個(gè)排隊(duì)系統(tǒng)中需要研究每個(gè)顧客i的到達(dá)時(shí)間 ai（或者到達(dá)時(shí)間間隔δi=ai?ai?1,a0=0）與服務(wù)臺(tái)對(duì)他們的服務(wù)時(shí)間 si。排隊(duì)系統(tǒng)主要有6 個(gè)特征：a. 到達(dá)時(shí)間的分布 fδ；b. 服務(wù)時(shí)間間隔的分布fs；c. 服務(wù)臺(tái)的個(gè)數(shù)K， K∈N；d. 系統(tǒng)的承載能力C，C∈N；e. 顧客數(shù)量n，n∈N；f. 服務(wù)規(guī)則 R。

對(duì)于到達(dá)間隔時(shí)間和服務(wù)時(shí)間的分布：常用M 表示指數(shù)或泊松且獨(dú)立分布；GI 表示廣義獨(dú)立分布；G 表示廣義無獨(dú)立性假設(shè)的分布。C 表示同一時(shí)刻系統(tǒng)承載的最多顧客數(shù)量。在系統(tǒng)中的顧客只有兩種狀態(tài)：等待或者正在被服務(wù)。顧客數(shù)量n 包含系統(tǒng)中、即將到達(dá)或者已經(jīng)離開的顧客。服務(wù)規(guī)則R 表示在隊(duì)列中顧客如何分配到對(duì)應(yīng)的服務(wù)臺(tái)，最普遍的規(guī)則是先到先服務(wù)（first come first serve，F(xiàn)CFS）。

對(duì)于單服務(wù)臺(tái)的排隊(duì)系統(tǒng)，由于顧客和服務(wù)臺(tái)只能處于兩種狀態(tài)，顧客等待服務(wù)臺(tái)或者服務(wù)臺(tái)等待顧客，因此第i 個(gè)顧客的離開時(shí)間定義為：di=max(ai,di?1)+si。

1.1 對(duì)于馬爾科夫系統(tǒng)M/M/1的貝葉斯推斷

在M/M/1排隊(duì)系統(tǒng)中，假定到達(dá)率和服務(wù)率λ 和 μ是未知的。 ta表示 na個(gè)顧客到達(dá)的總時(shí)間，ts表示 ns個(gè)顧客服務(wù)完成的時(shí)間，在排隊(duì)系統(tǒng)中一般na=ns。由于n個(gè)IID 指數(shù)分布的和服從Erlang分布，因此似然函數(shù)為

基于似然函數(shù)， λ 和μ的共軛先驗(yàn)分布為Gamma 分布。

其中αa,βa,αs,βs＞0，后驗(yàn)分布易得

1.2 對(duì)于非馬爾科夫系統(tǒng)的貝葉斯推斷

1.2.1 G/M/1模型

以Er/M/1為例，即到達(dá)間隔時(shí)間服從Erlang分布。在這個(gè)排隊(duì)系統(tǒng)中，假定到達(dá)時(shí)間是v 個(gè)以λ/ν為參數(shù)的指數(shù)同分布，則到達(dá)間隔時(shí)間X|ν,λ ～Er(ν,λ)，服務(wù)時(shí)間服從以μ為參數(shù)的指數(shù)分布，先驗(yàn)同M/M/1中Gamma 分布。(ν,λ)的共軛先驗(yàn)分布為

式中：λ ＞0；ν=1,2,···。根據(jù)該先驗(yàn)分布可以得到基于ν的 λ的條件分布、 ν的邊際分布分別為

默認(rèn)地設(shè)定θa=1,αa=βa=0，則先驗(yàn)變?yōu)閒(λ,ν)∝1/λ?；诠曹椣闰?yàn)分布，聯(lián)合后驗(yàn)分布為

1.2.2 M/G/1模型

以M/Er/1排隊(duì)系統(tǒng)為例，Insua 等[5]已經(jīng)提出對(duì)該模型的貝葉斯推斷過程。假定每個(gè)顧客到達(dá)的間隔時(shí)間是以 λ為參數(shù)的指數(shù)分布，服務(wù)時(shí)間服從Erlang 分布。指數(shù)分布的共軛先驗(yàn)為Gamma 分布，其后驗(yàn)分布也為Gamma 分布；服務(wù)時(shí)間服從愛爾蘭分布的先驗(yàn)不再贅述，且兩者相互獨(dú)立。

1.2.3 GI/GI/1模型

以GI/GI/1模型為例，假定到達(dá)間隔時(shí)間與服務(wù)時(shí)間都服從Coxian 分布或者混合廣義愛爾蘭（mixed generalized Erlang，MGE）分布，且兩者獨(dú)立。若X 服從Coxian 分布，則以L ，P=(P1,P2,···,PL)、λ=(λ1,λ2,···,λL)為參數(shù)，滿足

其中fr(x|λ1,...,λr)是MGE 分布的密度函數(shù)。若?i ≠j,s.t. λi≠λj，則

假定L 是已知的，本文基于Ausín 等[7]對(duì)Coxian 分布中各參數(shù)的無信息先驗(yàn)分布的結(jié)論作修改

特別地，令φr=1,?r=1,···,L。

1.3 基于近似貝葉斯計(jì)算的排隊(duì)系統(tǒng)

在貝葉斯推斷中M/G/1，G/M/1，GI/GI/1 模型的似然函數(shù)復(fù)雜，因此采用近似貝葉斯計(jì)算的方法，將重點(diǎn)從似然函數(shù)轉(zhuǎn)為人工數(shù)據(jù)的模擬，則需要選取恰當(dāng)?shù)南闰?yàn)與概要統(tǒng)計(jì)量。在衡量排隊(duì)系統(tǒng)好壞時(shí)，平均等待時(shí)間被指出為重要指標(biāo)，因此將平均等待時(shí)間作為第一個(gè)概要統(tǒng)計(jì)量，其中顧客等待時(shí)間可以直接從數(shù)據(jù)獲得，或定義為到達(dá)時(shí)間減去上一位顧客離開時(shí)間。并選取各排隊(duì)模型下離開時(shí)間的各分位數(shù)作為第二個(gè)概要統(tǒng)計(jì)量，能充分反映到達(dá)率和服務(wù)率在各排隊(duì)模型中的信息。

基于近似貝葉斯與QDC 思想，給出單服務(wù)臺(tái)的排隊(duì)系統(tǒng)的參數(shù)估計(jì)步驟如下：

a. 給出顧客排隊(duì)數(shù)據(jù)，并且假定顧客的到達(dá)時(shí)間a（ 或到達(dá)間隔時(shí)間）與接受服務(wù)的時(shí)間s分別服從參數(shù) λ 和μ，離開時(shí)間為 d，并給出閾值 ε0以及先驗(yàn)分布π(λ)和 π(μ)。

2 仿真與實(shí)際數(shù)據(jù)實(shí)驗(yàn)

R 軟件中的queue computer 包可以簡(jiǎn)單地將排隊(duì)過程呈現(xiàn)出來，并且可得到隊(duì)列的多個(gè)重要統(tǒng)計(jì)量，如平均等待時(shí)間、平均隊(duì)列長(zhǎng)度等。因此結(jié)合queue computer 與近似貝葉斯方法對(duì)4 個(gè)排隊(duì)模型進(jìn)行參數(shù)估計(jì)。

2.1 M/M/1排隊(duì)系統(tǒng)

由于先驗(yàn)分布服從Gamma 分布，因此給出選擇分布中超參數(shù)的3 種方法[12]：

b. 采用分位數(shù)，比如中位數(shù)和三分位數(shù)。檢驗(yàn)最終結(jié)果是否合理，并與其他不同的選擇情況進(jìn)行比較；

c. 合理地實(shí)驗(yàn)，通過經(jīng)驗(yàn)判斷或者嘗試后選取合適的 α 與 β。

這里選取第一種方法選取超參數(shù)，在該排隊(duì)系統(tǒng)中，對(duì)顧客數(shù)n=10 000，λ=1， μ=1/0.9的系統(tǒng)進(jìn)行模擬。在進(jìn)行30 次迭代后得到如表1 中的結(jié)果，后驗(yàn)分布如圖1。

表1 λ 和μ的估計(jì)值與標(biāo)準(zhǔn)差Tab.1 Estimation and standard deviation of λ and μ

從表1 可以得到 λ 和 μ的估計(jì)值與真實(shí)值接近，且后驗(yàn)分布也表明估計(jì)效果較好。

2.2 M/G/1排隊(duì)系統(tǒng)

在這個(gè)排隊(duì)系統(tǒng)中，對(duì)顧客數(shù)n=60，λ=1，Er(3,1.2)的系統(tǒng)進(jìn)行模擬。經(jīng)過實(shí)證比較，在選擇超參數(shù)時(shí)，選取中位數(shù)和三分位數(shù)較合理；在進(jìn)行30 次迭代后得到如表2 中的結(jié)果，后驗(yàn)分布如圖2。從表2 的結(jié)果可以觀察到雖然估計(jì)值都靠近真實(shí)值，結(jié)合后驗(yàn)分布及標(biāo)準(zhǔn)差發(fā)現(xiàn)μ的估計(jì)效果比 λ好。

圖1 λ 和μ的后驗(yàn)分布圖Fig.1 Posterior distribution of λ and μ

表2 λ 和μ的估計(jì)值與標(biāo)準(zhǔn)差Tab.2 Estimation and standard deviation of λ and μ

圖2 λ 和μ的后驗(yàn)分布圖Fig.2 Posterior distribution of λ and μ

2.3 G/M/1排隊(duì)系統(tǒng)

在這個(gè)排隊(duì)系統(tǒng)中，對(duì)顧客數(shù)n=100 ，Er(5,0.5)，μ=1的系統(tǒng)進(jìn)行模擬。經(jīng)過實(shí)證比較，在選擇超參數(shù)時(shí)，選取中位數(shù)和三分位數(shù)較合理；在進(jìn)行30 次迭代后得到如表3 中的結(jié)果，后驗(yàn)分布如圖3。從表3 的結(jié)果觀察到 λ和μ的估計(jì)值接近真實(shí)值，從后驗(yàn)分布與標(biāo)準(zhǔn)差的結(jié)果比較發(fā)現(xiàn)對(duì) λ的估計(jì)更準(zhǔn)確。

圖3 λ 和μ的后驗(yàn)分布圖Fig.3 Posterior distribution of λ and μ

表3 λ 和μ的估計(jì)值與標(biāo)準(zhǔn)差Tab.3 Estimation and standard deviation of λ and μ

2.4 GI/GI/1排隊(duì)系統(tǒng)

在這個(gè)排隊(duì)系統(tǒng)中，對(duì)顧客數(shù)n=300 ，L=4，P=(0.09,0.7,0.01,0.2)，λ=(1.1,1,0.251,0.25)的系統(tǒng)進(jìn)行模擬?；谛薷牡南闰?yàn)，在進(jìn)行50 次迭代后得到如表4 中的結(jié)果，后驗(yàn)分布如圖4。從表4中可以得到各分量的估計(jì)值與真實(shí)值相近，但 λ3和λ4估計(jì)值有偏離；結(jié)合后驗(yàn)分布發(fā)現(xiàn) λ2， λ3， λ4的估計(jì)效果相對(duì)較差，出現(xiàn)這種情況的原因可能在于Coxian 分布中參數(shù)的約束受到概率P 的影響。

表4 λ各分量的估計(jì)值與標(biāo)準(zhǔn)差Tab.4 Estimation and standard deviation of each component of λ

圖4 λ各個(gè)分量的后驗(yàn)分布圖Fig.4 Posterior distribution of each component of λ

2.5 實(shí)際數(shù)據(jù)應(yīng)用：銀行面對(duì)面排隊(duì)外匯交換數(shù)據(jù)

從http：//iew3.technion.ac.il/serveng 中可以直接得到銀行的面對(duì)面排隊(duì)外匯交換的數(shù)據(jù)，數(shù)據(jù)記錄了14 d 內(nèi)270 位顧客的排隊(duì)信息。這類服務(wù)需要一個(gè)服務(wù)臺(tái)在一周內(nèi)3 d 從8：30 工作到12：00；兩天從8：30 到12：30，16：00 到18：00。為了滿足QDC算法中連續(xù)性的要求，本文將所有信息根據(jù)到達(dá)時(shí)間升序排序，并做標(biāo)準(zhǔn)化處理。該數(shù)據(jù)中，平均到達(dá)間隔時(shí)間為11.99 min；平均服務(wù)時(shí)間為7.16 min。根據(jù)Ausín 等[7]對(duì)于L 的研究，在這里設(shè)定La=2 ，Ls=1，因此需要對(duì)3 個(gè)參數(shù) λ進(jìn)行估計(jì)，在50 次迭代后得到如表5 中的結(jié)果，后驗(yàn)分布如圖5。從表5 中的標(biāo)準(zhǔn)差及后驗(yàn)分布得出，λ1與 λ3的估計(jì)效果相對(duì)比 λ2好。

表5 λ1,λ2,λ3的估計(jì)值與標(biāo)準(zhǔn)差Tab.5 Estimation and standard deviation of λ1,λ2,λ3

為了進(jìn)一步驗(yàn)證參數(shù)估計(jì)值的正確性，將Ausín 等[7]研究得到的估計(jì)值λ1=0.1， λ2=0.1，λ3=0.14 與本文得到的值分別代入得的擬合數(shù)據(jù)、實(shí)際數(shù)據(jù)進(jìn)行比較（如圖6），比較發(fā)現(xiàn)本文估計(jì)參數(shù)擬合的數(shù)據(jù)與Ausin 的擬合數(shù)據(jù)分布接近，但兩者對(duì)于 λ1，λ2與實(shí)際數(shù)據(jù)都有差距，因此可能需要進(jìn)一步提高參數(shù)的精確度。

圖5 λ1,λ2,λ3的后驗(yàn)分布圖Fig.5 Posterior distribution of each component of λ1,λ2,λ3

圖6 兩組擬合數(shù)據(jù)的分布曲線圖Fig.6 Distribution of two fitted data

3 結(jié) 論

本文將近似貝葉斯與排隊(duì)論模型結(jié)合在一起，并在4 個(gè)排隊(duì)模型M/M/1， M/G/1，G/M/1，GI/GI/1下給出對(duì)應(yīng)的參數(shù)估計(jì)值以及估計(jì)效果。模擬研究表明，近似貝葉斯方法較排隊(duì)論傳統(tǒng)方法在參數(shù)估計(jì)上有優(yōu)勢(shì)，并且量化了不確定性，解決了復(fù)雜排隊(duì)模型中似然函數(shù)難以解析表達(dá)的困難。在實(shí)際應(yīng)用數(shù)據(jù)中也獲得了較好的參數(shù)估計(jì)，但為了提高在復(fù)雜模型中參數(shù)的精確度需要進(jìn)一步改進(jìn)貝葉斯參數(shù)估計(jì)方法。