上海理工大学学报  2020, Vol. 42 Issue (2): 108-114, 121   PDF    
基于近似贝叶斯计算方法的排队模型参数估计
张岚, 钱夕元     
华东理工大学,理学院,上海 200237
摘要: 基于贝叶斯推断对单服务台的马尔科夫与非马尔科夫系统下的4个排队模型做参数估计与效果评价,并对银行排队叫号换汇数据做实证研究。采用近似贝叶斯计算的方法有效解决复杂排队模型似然函数难解析表达的问题。利用R中queue computer包给出M/M/1,M/G/1,G/M/1,GI/GI/1这4个排队模型的到达时间与服务时间的参数估计及后验分布图。其中M/M/1:估计值(真实值)分别约为1.02(1)、1.12(1/0.9),都与真实值接近;M/G/1:1.17(1)、1.21(1.2),对服务时间的估计效果优于到达时间;G/M/1:0.49(0.5)、0.96(1),从后验分布进一步得到对服务时间的估计更准确;GI/GI/1:1.16(1.1)、1.07(1)、0.23(0.251)、0.22(0.25),各分量的估计值与真实值相对接近。对于实际银行数据,估计得的参数所拟合的数据与Ausin的拟合数据分布接近。研究表明,近似贝叶斯的方法在排队模型的参数估计上有较大优势,在实际数据中取得了较好应用。
关键词: 排队论     近似贝叶斯计算     参数估计    
Approximate bayesian computation based parameter estimation for queuing models
ZHANG Lan, QIAN Xiyuan     
School of Science, East China University of Science and Technology, Shanghai 200237, China
Abstract: Based on Bayesian inference, parameters and their validity of four queuing models with a single server in Markov and non-Markov system were estimated, and face-to-face currency exchange data of a bank was tested. By using approximate Bayesian computation, it is effective to solve the difficulty in analytically expressing the likelihood function in complex queuing models. Based on R package, parameters’ estimations and posterior distributions of arrival time and service time for four queuing models: M/M/1, M/G/1, G/M/1, GI/GI/1 were given. For M/M/1: estimated value (true value) is 1.02(1)、1.12(1/0.9) respectively, which are close to the true ones; for M/G/1: they are 1.17(1)、1.21(1.2) and estimated parameter of service time is superior to that of arrival time; for G/M/1: they are 0.49(0.5)、0.96(1) and estimated parameter of service time is more accurate from the posterior distribution; for GI/GI/1: they are 1.16(1.1)、1.07(1)、0.23(0.251)、0.22(0.25) and estimated value of each component is relatively close to the true value. Furthermore, for real data, distribution of fitting data of estimated parameters is close to that of Ausin’s fitting data. The results show great efficiency of parameter estimation based on approximate Bayesian computation in queuing theory and great application to the real world.
Key words: queuing theory     approximate Bayesian computation     parameter estimation    

排队论在社会各个领域都应用广泛,包括银行排队叫号、机场排队取行李、电话排队服务、交通运输等。但在实际应用中,在衡量一个排队模型优劣时,常需基于已知的到达率和服务率计算排队模型的几个重要指标:系统平均队长、排队等候平均队长、平均逗留时间、平均排队时间,而通常排队数据无法直接观测到这两个参数。在估计到达率和服务率,且是否使得排队系统稳定时,贝叶斯推断和估计在排队论中有很大的优势。

贝叶斯早在1763年就提出贝叶斯公式,并给出贝叶斯估计的思想,根据每一类的样本估计每一类的类条件概率密度,把参数 $\theta $ 看成未知的随机变量,通过对样本观察再求贝叶斯估计。在贝叶斯估计时,首先需要选择先验分布,再确定似然函数,确定参数的后验分布,选择损失函数,最后估计参数。贝叶斯推断在排队系统中的优势: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系统的贝叶斯推断,同时指出他们的方法可以用于估计瞬时队伍大小以及等待时间的分布。

在描述一个给定排队系统的情况时,常见方法DES(discrete event simulation)算法中状态的变化是不连续的。Ebert等[8]提出QDC(queue daparture computation)算法中状态变化是连续的,因此比DES更具有计算效率,且可结合近似贝叶斯计算的方法估计一个非常复杂的队列模型中的参数。由于在贝叶斯推断时,对于似然函数比较复杂的排队模型,贝叶斯参数估计就变得复杂,因此考虑不依赖传统似然的ABC(approximate Bayesian computation)方法。ABC最早由Feldman等[9]提出,将贝叶斯的后验分布公式 $\pi \left( {\theta |X} \right) = \dfrac{{L\left( {\theta |X} \right)\pi \left( \theta \right)}}{{\int {L\left( {\theta |X} \right)\pi \left( \theta \right)d\theta } }}$ 中对似然函数 $L\left( {\theta {\rm{|}}X} \right)$ 计算工作替换为可以产生人工模拟数据 $Y$ 的模型,其中 $\theta $ 为待估参数, $\pi \left( \theta \right)$ 为数据 $X$ 的先验分布。Malmberg等[10]指出在构造模拟数据 $Y$ 时需要建立与观测数据 $X$ 之间的测度(距离),常通过使观测数据与实际数据概要统计量的平方误差的和最小而进一步估计参数。令模拟数据 $Y$ 中参数满足 $\theta = {\theta ^ * }$ ,当两组数据的某个距离 $\rho \left( {X,Y} \right) < {\varepsilon _0}$ ,这里 ${\varepsilon _0}$ 足够小时, $\pi \left( {\theta |\rho \left( {X,Y} \right) < {\varepsilon _0}} \right)$ 可以近似当做后验。Brandon等[11]认为概要统计量需要充分提供 $\theta $ 在整个数据中的信息,如果 $S\left( X \right)$ $\theta $ 的概要统计量,则后验分布可以写成 $\pi \left( {\theta {\rm{|}}X} \right) = \pi \left( {\theta {\rm{|}}S\left( X \right)} \right)$ ,只有 ${\varepsilon _0}$ $\rho \left( {S\left( X \right),S\left( Y \right)} \right)$ 选择恰当时,ABC估计较准确。且根据Fisher-Neyman因子分解理论,如果概率分布 $f\left( {x|\theta } \right)$ 可以分解成 $f\left( {x|\theta } \right) = g\left( {S\left( x \right)|\theta } \right)h\left( x \right)$ ,那么概要统计量能充分说明 $\theta $ ,常见的 $S\left( X \right)$ 可以为期望或者方差。

本文运用近似贝叶斯的方法,结合排队论QDC算法的思想对四个经典排队模型进行仿真模拟,给出参数估计与估计效果;并结合实际例子(银行的面对面排队换汇数据)估计参数。

1 模 型

一个排队系统中需要研究每个顾客 $i$ 的到达时间 ${a_i}$ (或者到达时间间隔 ${\delta _i} = {a_i} - {a_{i - 1}},{a_0} = 0$ )与服务台对他们的服务时间 ${s_i}$ 。排队系统主要有6个特征:a. 到达时间的分布 ${f_\delta }$ ;b. 服务时间间隔的分布 ${f_s}$ ;c. 服务台的个数K $K$ $ \in \mathbb {N}$ ;d. 系统的承载能力C $C$ $ \in \mathbb {N}$ ;e. 顾客数量n $n$ $ \in \mathbb {N}$ ;f. 服务规则 $R$

对于到达间隔时间和服务时间的分布:常用M表示指数或泊松且独立分布;GI表示广义独立分布;G表示广义无独立性假设的分布。C表示同一时刻系统承载的最多顾客数量。在系统中的顾客只有两种状态:等待或者正在被服务。顾客数量n包含系统中、即将到达或者已经离开的顾客。服务规则R表示在队列中顾客如何分配到对应的服务台,最普遍的规则是先到先服务(first come first serve,FCFS)。

本文介绍4个单服务台排队系统: $M/M/1$ $M/G/1$ $G/M/1$ $GI/GI/1$ ,其中最简单的为: ${\delta _i}\mathop \sim \limits^{iid} \exp \left( \lambda \right), $ $\forall i \in 1:n$ 、服务时间 ${s_i}\mathop {\rm{\sim }}\limits^{iid} \exp \left( \mu \right),\forall i \in 1:n$ $\lambda $ $\mu $ 是指数分布的参数。设定K为1;Cn为无穷;R为FCFS;表示为 $M/M/1/\infty /\infty /\rm {FCFS}$ ,简写为 $M/M/1$ ,其中到达间隔时间也可以服从参数为 $\lambda $ 的泊松分布。类似地 $M/G/1$ 则表示到达间隔时间服从指数(或泊松)分布,服务时间服从广义分布。记 $\rho = {\lambda }/{K\mu }$ 为系统的资源利用率,只有 $\lambda < K\mu $ 时排队系统会达到平衡状态。若假设到达间隔时间与服务过程随时间的变化改变,即参数为 $\lambda \left( t \right)$ $\mu \left( t \right)$ ,系统表示为 $M\left( t \right)/M\left( t \right)/1$

对于单服务台的排队系统,由于顾客和服务台只能处于两种状态,顾客等待服务台或者服务台等待顾客,因此第i个顾客的离开时间定义为: ${d_i} = \max \left( {{a_i},{d_{i - 1}}} \right) + {s_i}$

1.1 对于马尔科夫系统 $M/M/1$ 的贝叶斯推断

$M/M/1$ 排队系统中,假定到达率和服务率 $\lambda $ $\mu $ 是未知的。 ${t_a}$ 表示 ${n_a}$ 个顾客到达的总时间, ${t_s}$ 表示 ${n_s}$ 个顾客服务完成的时间,在排队系统中一般 ${n_a} = {n_s}$ 。由于 $n$ 个IID指数分布的和服从Erlang分布,因此似然函数为

$\begin{split}&l\left( {\lambda ,\mu |data} \right) = \frac{{{\lambda ^{{n_a}}}}}{{\Gamma \left( {{n_a}} \right)}}{t_a}^{{n_a} - 1}{\rm {e}^{ - \lambda {t_a}}}\frac{{{\mu ^{{n_s}}}}}{{\Gamma \left( {{n_s}} \right)}}{t_s}^{{n_s} - 1}{\rm {e}^{ - \mu {t_s}}}\\&\quad \propto {\lambda ^{{n_a}}}{\rm {e}^{ - \lambda {t_a}}}{\mu ^{{n_s}}}{\rm {e}^{ - \mu {t_s}}}\end{split}$ (1)

基于似然函数, $\lambda $ $\mu $ 的共轭先验分布为Gamma分布。

$\left\{ {\begin{array}{*{20}{l}} \lambda {\rm{\sim }}Ga\left( {{\alpha _a},{\beta _a}} \right) \\ \mu {\rm{\sim }}Ga\left( {{\alpha _s},{\beta _s}} \right) \\ \end{array}} \right.$ (2)

其中 ${\alpha _a},{\beta _a},{\alpha _s},{\beta _s} > 0$ ,后验分布易得

$\left\{ {\begin{array}{*{20}{l}} \lambda {\rm{|}}data\sim Ga\left( {{\alpha _a}^ * ,{\beta _a}^ * } \right) \\ \mu {\rm{|}}data\sim Ga\left( {{\alpha _s}^ * ,{\beta _s}^ * } \right) \\ \end{array}} \right.$ (3)

其中 $\alpha _a^ * = {\alpha _a} + n,\beta _a^ * = {\beta _a} + {t_a},\;$ $\alpha _s^ * = {\alpha _s} + n,\beta _s^ * = {\beta _s} + {t_s}$ 。当先验信息较少时,可以采用Jeffreys先验, $f\left( {\lambda ,\mu } \right) \propto {1}/{{\lambda \mu }}$ ,则对应后验为

$\left\{ {\begin{array}{*{20}{l}} \lambda {\rm{|}}data\sim Ga\left( {{n_a},{t_a}} \right) \\ \mu {\rm{|}}data\sim Ga\left( {{n_s},{t_s}} \right) \\ \end{array}} \right. $ (4)
1.2 对于非马尔科夫系统的贝叶斯推断 1.2.1 $ G/M/1 $ 模型

${E_r}/M/1$ 为例,即到达间隔时间服从Erlang分布。在这个排队系统中,假定到达时间是v个以 $\lambda /\nu $ 为参数的指数同分布,则到达间隔时间 $X{\rm{|}}\nu ,\lambda \sim Er\left( {\nu ,\lambda } \right)$ ,服务时间服从以 $\mu $ 为参数的指数分布,先验同 $M/M/1$ 中Gamma分布。 $\left( {\nu ,\lambda } \right)$ 的共轭先验分布为

$f\left( {\nu ,\lambda } \right) \propto \frac{{\theta _a^{\nu - 1}\nu {{\left( {\lambda \nu } \right)}^{{\alpha _a}\nu - 1}}{e^{ - {\beta _a}\nu \lambda }}}}{{{{\left( {\left( {\nu - 1} \right)!} \right)}^{{\alpha _a}}}}}$ (5)

式中: $\lambda > 0{\text{;}}\nu = 1,2, \cdots $ 。根据该先验分布可以得到基于 $\nu $ $\lambda $ 的条件分布、 $\nu $ 的边际分布分别为

$\left\{ {\begin{array}{*{20}{l}} \lambda {\rm{|}}\nu {\rm{\sim }}Ga\left( {\nu {\alpha _a},{\beta _a}\nu } \right) \\ P\left( \nu \right) \propto \dfrac{{\Gamma \left( {{\alpha _a}\nu } \right)}}{{{{\left( {\left( {\nu - 1} \right)!} \right)}^{{\alpha _a}}}}}{\left( {\dfrac{{{\theta _a}}}{{\beta _a^{{\alpha _a}}}}} \right)^{\nu - 1}} \\ \end{array}} \right. $ (6)

默认地设定 ${\theta _a} = 1,{\alpha _a} = {\beta _a} = 0$ ,则先验变为 $f\left( {\lambda ,\nu } \right) \propto {1}/{\lambda }$ 。基于共轭先验分布,联合后验分布为

$\left\{ {\begin{array}{*{20}{l}} \lambda {\rm{|}}\nu ,data\sim Ga\left( {\nu \alpha _a^ * ,\beta _a^ * \nu } \right) \\ P\left( {\nu {\rm{|}}data} \right) \propto \dfrac{{\Gamma \left( {\alpha _a^ * \nu } \right)}}{{{{\left( {\left( {\nu - 1} \right)!} \right)}^{\alpha _a^ * }}}}{\left( {\dfrac{{\theta _a^ * }}{{\beta _a^{ * \alpha _a^ * }}}} \right)^{\nu - 1}} \\ \end{array}} \right. $ (7)

式中: $\alpha _a^ * = {\alpha _a} + n{\text{;}}\; \beta _a^ * = {\beta _a} + {t_a}{\text{;}}\;\theta _a^ * = {\theta _a}{T_a}$ ${T_a}$ 是到达间隔时间的乘积。

1.2.2 $ M/G/1 $ 模型

$M/{E_r}/1$ 排队系统为例,Insua等[5]已经提出对该模型的贝叶斯推断过程。假定每个顾客到达的间隔时间是以 $\lambda $ 为参数的指数分布,服务时间服从Erlang分布。指数分布的共轭先验为Gamma分布,其后验分布也为Gamma分布;服务时间服从爱尔兰分布的先验不再赘述,且两者相互独立。

1.2.3 $GI/GI/1$ 模型

$GI/GI/1$ 模型为例,假定到达间隔时间与服务时间都服从Coxian分布或者混合广义爱尔兰(mixed generalized Erlang,MGE)分布,且两者独立。若X服从Coxian分布,则以 $L$ $P = \left( {{P_1},{P_2}, \cdots ,{P_L}} \right)$ $\lambda = \left( {{\lambda _1},{\lambda _2}, \cdots ,{\lambda _L}} \right)$ 为参数,满足

$X = \left\{ {\begin{array}{*{20}{l}} {{Y_1}}&{{\text{依概率}}{P_1}}\\ {{Y_1} + {Y_2}}&{{\text{依概率}}{P_2}}\\ {{Y_1} + {Y_2}+\cdots+{Y_r}}&{{\text{依概率}}{P_r}}\\ \vdots &{}\\ {{Y_1} + {Y_2} + \cdots+{Y_r}+\cdots+ {Y_L}}&{{\text{依概率}}{P_L}} \end{array}} \right.$ (8)

式中: ${Y_r}$ 是以 ${\lambda _r}$ r= 1,2,…,L)为参数的指数分布; ${P_r},{\lambda _r} > 0$ 且满足 $\displaystyle\sum\limits_{r = 1}^L {{P_r}} = 1$ 。即可以表示为

$f\left( {x|L,P,\lambda } \right) = \sum\limits_{r = 1}^L {{P_r}} {f_r}\left( {x|{\lambda _1},{\lambda _2}, \ldots, {\lambda _r}} \right),x \geqslant 0$ (9)

其中 ${f_r}\left( {x|{\lambda _1}, \ldots, {\lambda _r}} \right)$ 是MGE分布的密度函数。若 $\exists i \ne j,s.t.\; \; {\lambda _i} \ne {\lambda _j}$ ,则

${f_r}\left( {x|{\lambda _1}, \cdots, {\lambda _r}} \right) = \sum\limits_{j = 1}^r {\left( {\prod\limits_{\begin{array}{*{20}{c}} {i = 1}\\ {i \ne j} \end{array}}^r {\frac{{{\lambda _i}}}{{{\lambda _i} - {\lambda _j}}}} } \right)} {\lambda _j}\exp \left( { - {\lambda _j}x} \right),x \geqslant 0$ (10)

假定L是已知的,本文基于Ausín等[7]对Coxian分布中各参数的无信息先验分布的结论作修改

$\left\{ {\begin{array}{*{20}{l}} \forall r = 1, \cdots ,L;P{\rm{\sim }}Dirichlet\left( {{\phi _1}, \cdots ,{\phi _L}} \right) \\ f\left( {{\lambda _r}} \right) \propto \dfrac{1}{{{\lambda _r}}} \\ \end{array}} \right. $ (11)

特别地,令 ${\phi _r} = 1,\forall r = 1, \cdots ,L$

1.3 基于近似贝叶斯计算的排队系统

在贝叶斯推断中M/G/1,G/M/1,GI/GI/1模型的似然函数复杂,因此采用近似贝叶斯计算的方法,将重点从似然函数转为人工数据的模拟,则需要选取恰当的先验与概要统计量。在衡量排队系统好坏时,平均等待时间被指出为重要指标,因此将平均等待时间作为第一个概要统计量,其中顾客等待时间可以直接从数据获得,或定义为到达时间减去上一位顾客离开时间。并选取各排队模型下离开时间的各分位数作为第二个概要统计量,能充分反映到达率和服务率在各排队模型中的信息。

基于近似贝叶斯与QDC思想,给出单服务台的排队系统的参数估计步骤如下:

a. 给出顾客排队数据,并且假定顾客的到达时间 $a$ (或到达间隔时间)与接受服务的时间 $s$ 分别服从参数 $\lambda $ $\mu $ ,离开时间为 $d$ ,并给出阈值 ${\varepsilon _0}$ 以及先验分布 $\pi \left( \lambda \right)$ $\pi \left( \mu \right)$

b. 对于顾客i:首先,从先验分布 $\pi \left( \lambda \right)$ $\pi \left( \mu \right)$ 中模拟抽取参数 ${\lambda ^{\rm{*}}}$ ${\mu ^ * }$ ;其次,基于参数 ${\lambda ^{\rm{*}}}$ ${\mu ^ * }$ 模拟产生新的到达和服务数据集,并根据 ${d_i} \leftarrow {\rm{max}} $ $\left( {{{{a}}_i},{d_{i - 1}}} \right) + {s_i}$ 计算新的离开时间集 ${d^ * }$ ;最后,选取平均等待时间、离开时间的分位数为概要统计量,并计算 $\left| {\overline w - {{\overline w }^ * }} \right|$ $\left| {{q_j}\left( d \right) - {q_j}\left( {{d^ * }} \right)} \right|,j = 1,2,3$ ,如果概要统计量的距离小于 ${\varepsilon _0}$ ,接受 ${\lambda ^{\rm{*}}}$ ${\mu ^ * }$ ,否则拒绝。

2 仿真与实际数据实验

R软件中的queue computer包可以简单地将排队过程呈现出来,并且可得到队列的多个重要统计量,如平均等待时间、平均队列长度等。因此结合queue computer与近似贝叶斯方法对4个排队模型进行参数估计。

2.1 $M/M/1$ 排队系统

由于先验分布服从Gamma分布,因此给出选择分布中超参数的3种方法[12]

a. 设均值为 $\mu $ 、方差为 ${\sigma ^2}$ ,由关系 $\mu = \frac{\alpha }{\beta }$ ${\sigma ^2} = {\alpha }/{{{\beta ^2}}}$ 推得: $\alpha = {{{\mu ^2}}}/{{{\sigma ^2}}}$ $\beta = {\mu }/{{{\sigma ^2}}}$

b. 采用分位数,比如中位数和三分位数。检验最终结果是否合理,并与其他不同的选择情况进行比较;

c. 合理地实验,通过经验判断或者尝试后选取合适的 $\alpha $ $\beta $

这里选取第一种方法选取超参数,在该排队系统中,对顾客数 $n = 10\;000$ $\lambda = 1$ $\mu = 1/0.9$ 的系统进行模拟。在进行30次迭代后得到如表1中的结果,后验分布如图1



表 1 $\lambda $ $\mu $ 的估计值与标准差 Table 1 Estimation and standard deviation of $\lambda $ and $\mu $

图 1 $\lambda $ $\mu $ 的后验分布图 Fig. 1 Posterior distribution of $\lambda $ and $\mu $

表1可以得到 $\lambda $ $\mu $ 的估计值与真实值接近,且后验分布也表明估计效果较好。

2.2 $M/G/1$ 排队系统

在这个排队系统中,对顾客数 $n = 60$ $\lambda = 1$ $Er\left( {3,1.2} \right)$ 的系统进行模拟。经过实证比较,在选择超参数时,选取中位数和三分位数较合理;在进行30次迭代后得到如表2中的结果,后验分布如图2。从表2的结果可以观察到虽然估计值都靠近真实值,结合后验分布及标准差发现 $\mu $ 的估计效果比 $\lambda $ 好。


图 2 $\lambda $ $\mu $ 的后验分布图 Fig. 2 Posterior distribution of $\lambda $ and $\mu $

表 2 $\lambda $ $\mu $ 的估计值与标准差 Table 2 Estimation and standard deviation of $\lambda $ and $\mu $
2.3 $G/M/1$ 排队系统

在这个排队系统中,对顾客数 $n = 100$ $Er\left( {5,0.5} \right)$ $\mu = 1$ 的系统进行模拟。经过实证比较,在选择超参数时,选取中位数和三分位数较合理;在进行30次迭代后得到如表3中的结果,后验分布如图3。从表3的结果观察到 $\lambda $ $\mu $ 的估计值接近真实值,从后验分布与标准差的结果比较发现对 $\lambda $ 的估计更准确。


图 3 $\lambda $ $\mu $ 的后验分布图 Fig. 3 Posterior distribution of $\lambda $ and $\mu $

表 3 $\lambda $ $\mu $ 的估计值与标准差 Table 3 Estimation and standard deviation of $\lambda $ and $\mu $
2.4 $GI/GI/1$ 排队系统

在这个排队系统中,对顾客数 $n = 300$ $L = 4$ $P = \left( {0.09,0.7,0.01,0.2} \right)$ $\lambda = \left( {1.1,1,0.251,0.25} \right)$ 的系统进行模拟。基于修改的先验,在进行50次迭代后得到如表4中的结果,后验分布如图4。从表4中可以得到各分量的估计值与真实值相近,但 ${\lambda _3}$ ${\lambda _4}$ 估计值有偏离;结合后验分布发现 ${\lambda _2}$ ${\lambda _3}$ ${\lambda _4}$ 的估计效果相对较差,出现这种情况的原因可能在于Coxian分布中参数的约束受到概率P的影响。



表 4 $\lambda $ 各分量的估计值与标准差 Table 4 Estimation and standard deviation of each component of $\lambda $

图 4 $\lambda $ 各个分量的后验分布图 Fig. 4 Posterior distribution of each component of $\lambda $
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的研究,在这里设定 ${L_a} = 2$ ${L_s} = 1$ ,因此需要对3个参数 $\lambda $ 进行估计,在50次迭代后得到如表5中的结果,后验分布如图5。从表5中的标准差及后验分布得出, ${\lambda _1}$ ${\lambda _3}$ 的估计效果相对比 ${\lambda _2}$ 好。



表 5 ${\lambda _1},{\lambda _2},{\lambda _3}$ 的估计值与标准差 Table 5 Estimation and standard deviation of ${\lambda _1},{\lambda _2},{\lambda _3}$

图 5 ${\lambda _1},{\lambda _2},{\lambda _3}$ 的后验分布图 Fig. 5 Posterior distribution of each component of ${\lambda _1},{\lambda _2},{\lambda _3}$

为了进一步验证参数估计值的正确性,将Ausín等[7]研究得到的估计值 ${\lambda _1} = 0.1$ ${\lambda _2} = 0.1$ ${\lambda _3} = $ 0.14 与本文得到的值分别代入得的拟合数据、实际数据进行比较(如图6),比较发现本文估计参数拟合的数据与Ausin的拟合数据分布接近,但两者对于 λ1λ2与实际数据都有差距,因此可能需要进一步提高参数的精确度。


图 6 两组拟合数据的分布曲线图 Fig. 6 Distribution of two fitted data
3 结 论

本文将近似贝叶斯与排队论模型结合在一起,并在4个排队模型 $M/M/1$ $M/G/1$ $G/M/1$ $GI/GI/1$ 下给出对应的参数估计值以及估计效果。模拟研究表明,近似贝叶斯方法较排队论传统方法在参数估计上有优势,并且量化了不确定性,解决了复杂排队模型中似然函数难以解析表达的困难。在实际应用数据中也获得了较好的参数估计,但为了提高在复杂模型中参数的精确度需要进一步改进贝叶斯参数估计方法。

参考文献
[1]
ARMERO C, BAYARRI M J. Bayesian prediction in M/M/1 queues [J]. Queueing Systems, 1994, 15(1/4): 401-417. DOI:10.1007/BF01189248
[2]
ARMERO C. Bayesian inference in Markovian queues[J]. Queueing Systems, 1994, 15(1/4): 419-426. DOI:10.1007/BF01189249
[3]
CHOUDHURY A, BORTHAKUR A C. Bayesian inference and prediction in the single server Markovian queue[J]. Metrika, 2008, 67(3): 371-383. DOI:10.1007/s00184-007-0138-3
[4]
WIPER M P. Bayesian analysis ofEr/M/1 and Er/M/c queues [J]. Journal of Statistical Planning and Inference, 1998, 69(1): 65-79. DOI:10.1016/S0378-3758(97)00124-9
[5]
INSUA D R, WIPER M, RUGGERI F. Bayesian analysis of M/Er/1 and M/Hk/1 queues [J]. Queueing Systems, 1998, 30(3/4): 289-308. DOI:10.1023/A:1019173206509
[6]
Ruggeri F, Insua D R. Bayesian model selection for M/G/1 queues[EB/OL]. (2007-10-16) [2018-11-02]. https://www.researchgate.net/publication/2469168.
[7]
AUSÍN M C, WIPER M P, LILLO R E. Bayesian prediction of the transient behaviour and busy period in short- and long-tailed GI/G/1 queueing systems [J]. Computational Statistics and Data Analysis, 2008, 52(3): 1615-1635. DOI:10.1016/j.csda.2007.05.009
[8]
EBERT A, WU P, MENGERSEN K, et al. Computationally efficient simulation of queues: the R package queuecomputer[EB/OL]. (2008-10-16)[2018-11-02]. https://cran.r-project.org/web/packages/queuecomputer/index.html.
[9]
PRITCHARD J K, SEIELSTAD M T, PEREZ-LEZAYN A, et al. Population growth of human Y chromosomes: a study of Y chromosome microsatellites[J]. Molecular Biology and Evolution, 1999, 16(12): 1791-1798. DOI:10.1093/oxfordjournals.molbev.a026091
[10]
MALMBERG K J, ZEELENBERG R, SHIFFRIN R M. Turning up the noise or turning down the volume? on the nature of the impairment of episodic recognition memory by Midazolam[J]. Journal of Experimental Psychology: Learning, Memory,and Cognition, 2004, 30(2): 540-549. DOI:10.1037/0278-7393.30(2):540
[11]
TURNER B M, VAN ZANDT T. A tutorial on approximate Bayesian computation[J]. Journal of Mathematical Psychology, 2012, 56(2): 69-85. DOI:10.1016/j.jmp.2012.02.005
[12]
INSUA D R, RUGGERI F, WIPER M P. Bayesian analysis of stochastic process models[M]. Hoboken, N.J.: Wiley, 2012.