基于近似貝葉斯計算方法的排隊模型參數估計

張 嵐， 錢夕元

（華東理工大學，理學院，上海 200237）

排隊論在社會各個領域都應用廣泛，包括銀行排隊叫號、機場排隊取行李、電話排隊服務、交通運輸等。但在實際應用中，在衡量一個排隊模型優劣時，常需基于已知的到達率和服務率計算排隊模型的幾個重要指標：系統平均隊長、排隊等候平均隊長、平均逗留時間、平均排隊時間，而通常排隊數據無法直接觀測到這兩個參數。在估計到達率和服務率，且是否使得排隊系統穩定時，貝葉斯推斷和估計在排隊論中有很大的優勢。

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

本文運用近似貝葉斯的方法，結合排隊論QDC算法的思想對四個經典排隊模型進行仿真模擬，給出參數估計與估計效果；并結合實際例子（銀行的面對面排隊換匯數據）估計參數。

1 模 型

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

對于到達間隔時間和服務時間的分布：常用M 表示指數或泊松且獨立分布；GI 表示廣義獨立分布；G 表示廣義無獨立性假設的分布。C 表示同一時刻系統承載的最多顧客數量。在系統中的顧客只有兩種狀態：等待或者正在被服務。顧客數量n 包含系統中、即將到達或者已經離開的顧客。服務規則R 表示在隊列中顧客如何分配到對應的服務臺，最普遍的規則是先到先服務（first come first serve，FCFS）。

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

1.1 對于馬爾科夫系統M/M/1的貝葉斯推斷

在M/M/1排隊系統中，假定到達率和服務率λ 和 μ是未知的。 ta表示 na個顧客到達的總時間，ts表示 ns個顧客服務完成的時間，在排隊系統中一般na=ns。由于n個IID 指數分布的和服從Erlang分布，因此似然函數為

基于似然函數， λ 和μ的共軛先驗分布為Gamma 分布。

其中αa,βa,αs,βs＞0，后驗分布易得

1.2 對于非馬爾科夫系統的貝葉斯推斷

1.2.1 G/M/1模型

以Er/M/1為例，即到達間隔時間服從Erlang分布。在這個排隊系統中，假定到達時間是v 個以λ/ν為參數的指數同分布，則到達間隔時間X|ν,λ ～Er(ν,λ)，服務時間服從以μ為參數的指數分布，先驗同M/M/1中Gamma 分布。(ν,λ)的共軛先驗分布為

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

默認地設定θa=1,αa=βa=0，則先驗變為f(λ,ν)∝1/λ。基于共軛先驗分布，聯合后驗分布為

1.2.2 M/G/1模型

以M/Er/1排隊系統為例，Insua 等[5]已經提出對該模型的貝葉斯推斷過程。假定每個顧客到達的間隔時間是以 λ為參數的指數分布，服務時間服從Erlang 分布。指數分布的共軛先驗為Gamma 分布，其后驗分布也為Gamma 分布；服務時間服從愛爾蘭分布的先驗不再贅述，且兩者相互獨立。

1.2.3 GI/GI/1模型

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

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

假定L 是已知的，本文基于Ausín 等[7]對Coxian 分布中各參數的無信息先驗分布的結論作修改

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

1.3 基于近似貝葉斯計算的排隊系統

在貝葉斯推斷中M/G/1，G/M/1，GI/GI/1 模型的似然函數復雜，因此采用近似貝葉斯計算的方法，將重點從似然函數轉為人工數據的模擬，則需要選取恰當的先驗與概要統計量。在衡量排隊系統好壞時，平均等待時間被指出為重要指標，因此將平均等待時間作為第一個概要統計量，其中顧客等待時間可以直接從數據獲得，或定義為到達時間減去上一位顧客離開時間。并選取各排隊模型下離開時間的各分位數作為第二個概要統計量，能充分反映到達率和服務率在各排隊模型中的信息。

基于近似貝葉斯與QDC 思想，給出單服務臺的排隊系統的參數估計步驟如下：

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

2 仿真與實際數據實驗

R 軟件中的queue computer 包可以簡單地將排隊過程呈現出來，并且可得到隊列的多個重要統計量，如平均等待時間、平均隊列長度等。因此結合queue computer 與近似貝葉斯方法對4 個排隊模型進行參數估計。

2.1 M/M/1排隊系統

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

b. 采用分位數，比如中位數和三分位數。檢驗最終結果是否合理，并與其他不同的選擇情況進行比較；

c. 合理地實驗，通過經驗判斷或者嘗試后選取合適的 α 與 β。

這里選取第一種方法選取超參數，在該排隊系統中，對顧客數n=10 000，λ=1， μ=1/0.9的系統進行模擬。在進行30 次迭代后得到如表1 中的結果，后驗分布如圖1。

表1 λ 和μ的估計值與標準差Tab.1 Estimation and standard deviation of λ and μ

從表1 可以得到 λ 和 μ的估計值與真實值接近，且后驗分布也表明估計效果較好。

2.2 M/G/1排隊系統

在這個排隊系統中，對顧客數n=60，λ=1，Er(3,1.2)的系統進行模擬。經過實證比較，在選擇超參數時，選取中位數和三分位數較合理；在進行30 次迭代后得到如表2 中的結果，后驗分布如圖2。從表2 的結果可以觀察到雖然估計值都靠近真實值，結合后驗分布及標準差發現μ的估計效果比 λ好。

圖1 λ 和μ的后驗分布圖Fig.1 Posterior distribution of λ and μ

表2 λ 和μ的估計值與標準差Tab.2 Estimation and standard deviation of λ and μ

圖2 λ 和μ的后驗分布圖Fig.2 Posterior distribution of λ and μ

2.3 G/M/1排隊系統

在這個排隊系統中，對顧客數n=100 ，Er(5,0.5)，μ=1的系統進行模擬。經過實證比較，在選擇超參數時，選取中位數和三分位數較合理；在進行30 次迭代后得到如表3 中的結果，后驗分布如圖3。從表3 的結果觀察到 λ和μ的估計值接近真實值，從后驗分布與標準差的結果比較發現對 λ的估計更準確。

圖3 λ 和μ的后驗分布圖Fig.3 Posterior distribution of λ and μ

表3 λ 和μ的估計值與標準差Tab.3 Estimation and standard deviation of λ and μ

2.4 GI/GI/1排隊系統

在這個排隊系統中，對顧客數n=300 ，L=4，P=(0.09,0.7,0.01,0.2)，λ=(1.1,1,0.251,0.25)的系統進行模擬。基于修改的先驗，在進行50 次迭代后得到如表4 中的結果，后驗分布如圖4。從表4中可以得到各分量的估計值與真實值相近，但 λ3和λ4估計值有偏離；結合后驗分布發現 λ2， λ3， λ4的估計效果相對較差，出現這種情況的原因可能在于Coxian 分布中參數的約束受到概率P 的影響。

表4 λ各分量的估計值與標準差Tab.4 Estimation and standard deviation of each component of λ

圖4 λ各個分量的后驗分布圖Fig.4 Posterior distribution of each component of λ

2.5 實際數據應用：銀行面對面排隊外匯交換數據

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

表5 λ1,λ2,λ3的估計值與標準差Tab.5 Estimation and standard deviation of λ1,λ2,λ3

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

圖5 λ1,λ2,λ3的后驗分布圖Fig.5 Posterior distribution of each component of λ1,λ2,λ3

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

3 結 論

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