# Addressing Deep Learning Model Calibration Using Evidential Neural Networks
and Uncertainty-Aware Training

###### Abstract

In terms of accuracy, deep learning (DL) models have had considerable success
in classification problems for medical imaging applications. However, it is
well-known that the outputs of such models, which typically utilise the
SoftMax function in the final classification layer can be over-confident, i.e.
they are poorly _calibrated_. Two competing solutions to this problem have
been proposed: uncertainty-aware training and evidential neural networks
(ENNs). In this paper we perform an investigation into the improvements to
model calibration that can be achieved by each of these approaches
individually, and their combination. We perform experiments on two
classification tasks: a simpler MNIST digit classification task and a more
complex and realistic medical imaging artefact detection task using Phase
Contrast Cardiac Magnetic Resonance images. The experimental results
demonstrate that model calibration can suffer when the task becomes
challenging enough to require a higher capacity model. However, in our complex
artefact detection task we saw an improvement in calibration for both a low
and higher capacity model when implementing both the ENN and uncertainty-aware
training together, indicating that this approach can offer a promising way to
improve calibration in such settings. The findings highlight the potential use
of these approaches to improve model calibration in a complex application,
which would in turn improve clinician trust in DL models.

Index Terms—  Calibration, Evidential Neural Networks, Medical Imaging

##  1 Introduction

In recent years, artificial intelligence (AI), and more specifically deep
learning (DL) models have dominated medical imaging research for various
classification problems [1]. In many such applications, DL will likely be
utilised as a decision support tool to enable clinical users to make informed
decisions with AI assistance that comply with evidence-based medical practice
[2]. In this context, as well as the model prediction, the confidence of this
prediction becomes a critical component of the AI model.

The concept of _model calibration_ refers to the relationship between the
accuracy of predictions and their confidence: a well-calibrated model will be
less confident when making wrong predictions and more confident when making
correct predictions. With this in mind, measures of model calibration have
been proposed that provide a more complete understanding of the performance of
a model by estimating how closely the confidence matches the accuracy [3].

Good model calibration is of particular importance in medical imaging-based
AI. Patient prognosis may be adversely affected if poorly calibrated models
are deployed in a clinical setting, e.g. if models produce confident
predictions for incorrect diagnoses. Therefore it is important that AI models
for medical imaging are evaluated with appropriate metrics that can quantify
calibration to more suitably assess a model’s performance. Well-calibrated
models for medical AI applications can provide an opportunity to develop safe
and more reliable clinical decision-support tools [4].

It is well-known that DL models which utilise the SoftMax function in their
final layer can be poorly calibrated [5, 6]. In the literature, two types of
approaches have been proposed to attempt to improve the calibration
performance of DL models: evidential neural networks (ENN) and uncertainty-
aware training. ENNs [7] are deterministic models that remove the dependency
on the point estimates of the predictive distribution for confidence or
uncertainty when using the SoftMax function. ENNs utilise probability
distributions to obtain uncertainty in predictions for each category/class by
gathering evidence during training using a Dirichlet distribution and have
shown promising results towards improving model calibration in computer vision
tasks [8, 9, 10, 11, 12].

Uncertainty-aware training involves modifying the training of a standard DL
model (e.g. by changing the loss function) to explicitly encourage weight
updates that improve calibration. For example, [13] proposed an additional
loss term based on the relationship between accuracy and uncertainty. [9]
utilised a similar approach but tailored their loss to be more suitable for
their particular computer vision application. In the medical imaging domain
interest in uncertainty-aware training has also grown with recent work [14]
proposing a technique for improving calibration of a model for predicting
response to treatment from cardiac magnetic resonance images. In medical image
segmentation a novel uncertainty-aware framework was developed with a student
and teacher model, where the student model learns from the teacher model by
minimizing a segmentation loss and a consistency loss [15].

Contributions: In this paper, for the first time, we investigate and compare
the use of ENNs and uncertainty-aware training (as well as their combination)
to improve model calibration, and furthermore evaluate performance using both
a low capacity model (LeNet5) and a higher capacity model (ResNet18). We
perform two experiments to test our approaches: (1) a simple application in
which we aim to classify handwritten digits using the publicly available MNIST
database; and (2) a more complex and realistic medical application, in which
we aim to detect the presence of artefacts in 2D flow Phase Contrast Cardiac
Magnetic Resonance (PC-CMR) cine images.

##  2 Materials and Methods

Two different experiments were performed to evaluate the impact of different
methods to address the problem of model calibration. For the first experiment,
we used the publicly available MNIST digit dataset, which contains 50,000
training. 10,000 validation and 10,000 test images. For the second experiment,
we utilised two retrospective 2D flow PC-CMR databases, in which all images
were manually annotated for the presence of artefacts. The first database
contained images from 80 patients from the UK Biobank [16] and the second
database was retrieved from the clinical imaging database of Guy’s and St
Thomas’ NHS Foundation Trust (GSTFT) and contained images from 830 patients
with a range of different cardiovascular diseases acquired under routine
clinical PC-CMR practice. For the second experiment, we combined the two
databases and used 80% of the total data for training, 10% for validation, and
10% for testing. Further details about the data and protocol can be found in
[17].

We evaluated the performance of two different DL models: a low-capacity model
(LeNet5 [18]) and a higher capacity model (ResNet18 [19]). To address model
calibration, we utilised three approaches. The first approach was an ENN [7].
The second approach was the _Uncertainty versus Accuracy_ (UvAC) uncertainty-
aware training approach, which uses a differentiable loss function that
quantifies the relationship between accuracy and uncertainty [13]. This loss
function aims to improve uncertainty calibration and can be used as an
additional penalty term in model training. The final approach was a
combination of the ENN and UvAC. The aim here was to determine if the ENN’s
more robust estimate of uncertainty would improve the performance of the
uncertainty-aware training method.

##  3 Experiments and Results

We used the test set accuracy to measure model performance. To quantify model
calibration we utilised the Adaptive Expected Calibration Error (AECE) metric
[20], again computed on the test set.

Experiment 1: MNIST
In this experiment our motivation was to systematically modify model
performance and calibration by adding noise of varying levels, and then to
observe the impact of different methods aimed at improving calibration. After
normalising the image intensities to a mean of zero and standard deviation of
one, we added Gaussian noise (with mean of zero and different standard
deviations) to both the train and test sets, and for each noise level the test
accuracy and AECE were computed for each evaluated model. The values of the
learning rate, number of epochs and batch size were taken from typical values
in the literature: $\mathrm{10^{-3}}$, $10$ and $128$ (research has shown that
the MNIST dataset is easily trainable with low error rates seen at 3-6 epochs
[21]). The AvUC loss weight and ENN scaling factor were optimised using a grid
search strategy (ranges 0-4 and 0-40 respectively). See Table 1 (left) for
final hyperparameter values, which were selected to maximise validation set
accuracy. Table 1 (left) shows the results from this experiment, which are
averaged over 10 train/test runs for each approach.

Table 1: Summary of hyperparameters for Experiment 1 (MNIST) and Experiment 2
(PC-CMR) for all approaches using the LeNet and ResNet18 architectures.

| MNIST | PC-CMR
---|---|---
Method | AvUC | ENN | AvUC | ENN
LeNet5+UvAC | 2 | - | 1.5 | -
LeNet5+ENN | - | 30 | - | 10
LeNet5+ENN+UvAC | 1 | 40 | 2 | 10
ResNet18+UvAC | 2 | - | 2 | -
ResNet18+ENN | - | 30 | - | 10
ResNet18+ENN+UvAC | 2 | 40 | 1 | 10

Experiment 2: Phase Contrast Cardiac Magnetic Resonance (PC-CMR) Images
Here we utilised a realistic clinical data set to address a more complex
problem and observe the outcomes from our approaches. The problem of artefact
detection from PC-CMR was chosen because these images are utilised by
clinicians to quantify and assess cardiac function, but they can contain
artefacts which make their analysis unreliable. Identifying these artefacts
using DL models is a challenging problem. Similar to Experiment 1 we
quantified performance using the test accuracy and AECE. The values of the
learning rate, number of epochs and batch size ($\mathrm{10^{-3}}$, $200$ and
$16$ respectively) were taken from from values used in [17], whilst the AvUC
loss weight and ENN scaling factor were optimised using a grid search strategy
(ranges 0-4 and 0-40 respectively). See Table 1 (right) for final
hyperparameters, which were selected to maximise validation set accuracy.
Table 2 (right) shows the results from this experiment, which were again
averaged over 10 train/test runs for each approach. Sample reliability
diagrams using AECE are shown in Figure 1, which illustrate the relationship
between confidence and accuracy for all four approaches for the artefact
classification problem.

(a) ResNet18

(b) ResNet18+UvAC

(c) ResNet18+ENN

(d) ResNet18+ENN+UvAC

Fig. 1: AECE reliability diagrams for PC-CMR data with observed improvement
wrt AECE when training with (a) baseline ResNet18 model versus (b) ResNet18
with the UvAC loss only (c) ResNet18 with ENN only and lastly (d) ResNet18
with ENN and UvAC loss.

##  4 Discussion and Conclusion

We observe that there is inconsistency of performance across the two
experiments. On the MNIST experiment, both the LeNet5 and ResNet18 models show
high accuracy, but the ResNet models were more robust to noise. Both the
LeNet5 and ResNet18 baseline models were well calibrated, and applying the ENN
and/or uncertainty-aware training did not improve calibration and sometimes
degraded performance, both in terms of accuracy and calibration.

For the PC-CMR experiment, the higher capacity ResNet18 model performed better
than the low capacity LeNet5 model in terms of accuracy. For the ResNet18, the
ENN alone and the ENN combined with UvAC uncertainty aware training provided a
significant improvement to model calibration but the UvAC alone did not. In
contrast, for the LeNet5, only the combination of the ENN and UvAC
uncertainty-aware training improved calibration, whilst the ENN or UvAC alone
did not. The only approach that improved calibration for both the LeNet5 and
ResNet18 models was the combination of the ENN and UvAC uncertainty-aware
training.

We speculate that model calibration can suffer when the task becomes more
challenging, requiring a higher capacity model (as in the PC-CMR experiment,
compared to the easier MNIST task). In such situations it can be beneficial to
apply a method to improve model calibration. For example, noticeable
improvements are evident in the reliability diagrams in Figure 1 c-d over the
baseline ResNet18 in Figure 1a. Based on our results it seems as if the best
way to improve model calibration may be dependent on the capacity of the
model, although the combination of ENN and uncertainty-aware training seems to
be robust to differences in model capacity.

To the best of our knowledge the combination of an ENN with uncertainty-aware
training has not previously been investigated in a medical imaging context but
such a combination has been proposed in a computer vision context with some
preliminary success [9]. Our results indicate that a robust and deterministic
evidence-based and uncertainty-aware learning approach may improve calibration
of DL models by improving the modelling of uncertainty. In future work we will
investigate other uncertainty-aware training approaches in order to develop
trustworthy and reliable DL models for healthcare applications.

Table 2: Accuracy and calibration results of all models applied to the MNIST
and PC-CMR datasets. Evaluation metrics are: test accuracy (Acc) and Adaptive
Expected Calibration Error (AECE) in blue. Noise level indicated is the
standard deviation of the Gaussian distribution. All results averaged over 10
train/test runs.

MNIST (Acc, AECE) PC-CMR (Acc, AECE) Method/Noise levels: 0 0.2 0.3 0.5 6
LeNet5 98.9 $\pm$ $0.1$, 0.003 $\pm$ $0.001$ 98.7 $\pm$ $0.05$, 0.003 $\pm$
$0.001$ 98.7 $\pm$ $0.07$, 0.003 $\pm$ $0.0005$ 98.3 $\pm$ $0.12$, 0.003 $\pm$
$0.001$ 62.5 $\pm$ $0.13$ , 0.01 $\pm$ $0.005$ 52.0 $\pm$ $4.6$, 0.12 $\pm$
$0.02$ LeNet5+UvAC 99.0 $\pm$ $0.048$, 0.004 $\pm$ $0.001$ 98.9 $\pm$ $0.011$,
0.005 $\pm$ $0.001$ 98.8 $\pm$ $0.11$ , 0.004 $\pm$ $0.001$ 99.0 $\pm$ $0.10$,
0.004 $\pm$ $0.001$ 63.1 $\pm$ $0.01$, 0.01 $\pm$ $0.002$ 50.0 $\pm$ $7.6$,
0.11 $\pm$ $0.05$ LeNet5+ENN 98.0 $\pm$ $0.06$, 0.05 $\pm$ $0.04$ 86.70 $\pm$
$12.0$, 0.12 $\pm$ $0.10$ 85.30 $\pm$ $9.2$, 0.14 $\pm$ $0.10$ 89.30 $\pm$
$8.4$ 0.1 $\pm$ $0.08$ 57.40 $\pm$ $2.8$, 0.21 $\pm$ $0.03$ 54.60 $\pm$ $7.0$,
0.12 $\pm$ $0.05$ LeNet5+ENN+UvAC 97.1 $\pm$ $2.90$, 0.03 $\pm$ $0.03$ 93.20
$\pm$ $4.37$ 0.06 $\pm$ $0.04$ 82.04 $\pm$ $7.42$, 0.17 $\pm$ $0.06$ 84.91
$\pm$ $9.47$, 0.14 $\pm$ $0.10$ 50.89 $\pm$ $6.55$, 0.33 $\pm$ $0.04$ 55.1
$\pm$ $8.9$, 0.06 $\pm$ $0.03$ ResNet18 99.0 $\pm$ $0.04$, 0.003 $\pm$ $0.001$
98.9 $\pm$ $0.09$, 0.003 $\pm$ $0.001$ 98.9 $\pm$ $0.09$, 0.003 $\pm$ $0.001$
98.8 $\pm$ $0.14$, 0.002 $\pm$ $0.001$ 61.6, $\pm$ $0.37$, 0.01 $\pm$ $0.002$
64.8 $\pm$ $5$, 0.17 $\pm$ $0.05$ ResNet18+UvAC 99.0 $\pm$ $0.005$, 0.003
$\pm$ $0.001$ 99.0 $\pm$ $0.004$, 0.003 $\pm$ $0.001$ 99.0 $\pm$ $0.11$, 0.003
$\pm$ $0.001$ 99.0 $\pm$ $0.1$, 0.002 $\pm$ $0.001$ 62.0 $\pm$ $0.37$, 0.01
$\pm$ $0.002$ 64.5$\pm$ $3.2$, 0.17 $\pm$ $0.04$ ResNet18+ENN 99.0 $\pm$
$0.35$, 0.01 $\pm$ $0.01$ 99.2 $\pm$ $0.08$, 0.008 $\pm$ $0.002$ 99.0 $\pm$
$0.09$, 0.01 $\pm$ $0.01$ 98.4 $\pm$ $1.58$, 0.01 $\pm$ $0.01$ 61.0 $\pm$
$3.13$, 0.13 $\pm$ $0.01$ 66.0 $\pm$ $4.90$, 0.06 $\pm$ $0.02$
ResNet18+ENN+UvAC 98.7 $\pm$ $0.35$, 0.05 $\pm$ $0.03$ 98.8 $\pm$ $0.13$, 0.05
$\pm$ $0.01$ 99.0 $\pm$ $0.17$, 0.06 $\pm$ $0.01$ 98.7 $\pm$ $0.23$, 0.07
$\pm$ $0.03$ 61.4 $\pm$ $0.53$, 0.31 $\pm$ $0.05$ 66.0 $\pm$ $4.70$, 0.06
$\pm$ $0.02$

##  5 Compliance with ethical standards

This research study was conducted retrospectively using human subject data
made available using the UK Biobank Resource under Application Number 17806.

##  6 Acknowledgments

This work was supported by the Kings DRIVE Health CDT for Data-Driven Health
and further funded/supported by the National Institute for Health Research
(NIHR) Biomedical Research Centre at Guy’s and St Thomas’ NHS Foundation Trust
and King’s College London. Additionally this research was funded in whole, or
in part, by the Wellcome Trust WT203148/Z/16/Z. For the purpose of open
access, the author has applied a CC BY public copyright license to any author
accepted manuscript version arising from this submission. This research has
been conducted using the UK Biobank Resource under Application Number 17806.

The views expressed in this paper are those of the authors and not necessarily
those of the NHS, EPSRC, the NIHR or the Department of Health and Social Care.

## References

  * [1] Xuxin Chen, Ximin Wang, Ke Zhang, et al.,  “Recent advances and clinical applications of deep learning in medical image analysis,”  Medical Image Analysis, p. 102444, 2022.
  * [2] Weina Jin, Xiaoxiao Li, Mostafa Fatehi, and Ghassan Hamarneh,  “Guidelines and evaluation for clinical explainable ai on medical image analysis,”  arXiv preprint arXiv:2202.10553, 2022.
  * [3] Jeremy Nixon, Michael W Dusenberry, Linchuan Zhang, et al.,  “Measuring calibration in deep learning.,”  in CVPR Workshops, 2019, vol. 2.
  * [4] Annika Reinke, Minu D Tizabi, Michael Baumgartner, Matthias Eisenmann, Doreen Heckmann-Nötzel, A Emre Kavu, Tim Rädsch, Carole H Sudre, Laura Acion, Michela Antonelli, et al.,  “Understanding metric-related pitfalls in image analysis validation,”  arXiv preprint arXiv:2302.01790, 2023.
  * [5] Benjamin Kompa, Jasper Snoek, and Andrew L Beam,  “Second opinion needed: communicating uncertainty in medical machine learning,”  NPJ Digital Medicine, vol. 4, no. 1, pp. 1–6, 2021.
  * [6] Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, et al.,  “A survey of uncertainty in deep neural networks,”  arXiv preprint arXiv:2107.03342, 2021.
  * [7] Murat Sensoy, Lance Kaplan, and Melih Kandemir,  “Evidential deep learning to quantify classification uncertainty,”  Advances in neural information processing systems, vol. 31, 2018\.
  * [8] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger,  “On calibration of modern neural networks,”  in International conference on machine learning. PMLR, 2017, pp. 1321–1330.
  * [9] Wentao Bao, Qi Yu, and Yu Kong,  “Evidential deep learning for open set action recognition,”  in Proceedings of the IEEE/CVF International Conference on Computer Vision, 2021, pp. 13349–13358.
  * [10] Pauline Minary, Frédéric Pichon, et al.,  “Face pixel detection using evidential calibration and fusion,”  International Journal of Approximate Reasoning, vol. 91, pp. 202–215, 2017.
  * [11] Christian Tomani and Florian Buettner,  “Towards trustworthy predictions from deep neural networks with fast adversarial calibration,”  in Proceedings of the AAAI Conference on Artificial Intelligence, 2021, vol. 35, pp. 9886–9896.
  * [12] Archit Karandikar, Nicholas Cain, et al.,  “Soft calibration objectives for neural networks,”  Advances in Neural Information Processing Systems, vol. 34, 2021\.
  * [13] Ranganath Krishnan and Omesh Tickoo,  “Improving model calibration with accuracy versus uncertainty optimization,”  Advances in Neural Information Processing Systems, 2020.
  * [14] Tareen Dawood, Chen Chen, Robin Andlauer, et al.,  “Uncertainty-aware training for cardiac resynchronisation therapy response prediction,”  in Statistical Atlases and Computational Models of the Heart. 2022, pp. 189–198, Springer International Publishing.
  * [15] Lequan Yu, Shujun Wang, Xiaomeng Li, Chi-Wing Fu, and Pheng-Ann Heng,  “Uncertainty-aware self-ensembling model for semi-supervised 3D left atrium segmentation,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2019, pp. 605–613.
  * [16] Cathie Sudlow, John Gallacher, Naomi Allen, et al.,  “Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age,”  PLoS medicine, vol. 12, no. 3, pp. e1001779, 2015.
  * [17] Emily Chan, Ciaran O’Hanlon, Carlota Asegurado Marquez, et al.,  “Automated quality controlled analysis of 2d phase contrast cardiovascular magnetic resonance imaging,”  in Statistical Atlases and Computational Models of the Heart (In Press), 2022.
  * [18] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner,  “Gradient-based learning applied to document recognition,”  Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.
  * [19] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [20] Yukun Ding, Jinglan Liu, Jinjun Xiong, and Yiyu Shi,  “Revisiting the evaluation of uncertainty estimation and its application to explore model complexity-uncertainty trade-off,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, 2020, pp. 4–5.
  * [21] Dan Claudiu Ciresan, Ueli Meier, Jonathan Masci, Luca Maria Gambardella, and Jürgen Schmidhuber,  “Flexible, high performance convolutional neural networks for image classification,”  in Twenty-second international joint conference on artificial intelligence, 2011.


0

Zuk and Kirszenblat

Non-Preemptive Two-Level Priority Queue

Analytic Approach to the Non-Preemptive Markovian Priority Queue

Josef Zuk Defence Science and Technology Group, Melbourne, Australia,
<EMAIL_ADDRESS>David Kirszenblat Defence Science and Technology
Group, Melbourne, Australia<EMAIL_ADDRESS>
A new approach is developed for the joint queue-length distribution of the
two-level non-preemptive Markovian priority queue that allows explicit and
exact results to be obtained. Marginal distributions are derived for the
general multi-level problem. The results are based on a representation of the
joint queue-length probability mass function as a single-variable complex
contour integral, that reduces to a real integral on a finite interval arising
from a cut on the real axis. Both numerical quadrature rules and exact finite
sums, involving Legendre polynomials and their generalization, are presented
for the joint and marginal distributions. A high level of accuracy is
demonstrated across the entire ergodic region. Relationships are established
with the waiting-time distributions. Asymptotic behaviour in the large queue-
length regime is extracted.

queueing theory; non-preemptive priority; queue length distribution Primary:
90B22; secondary: 60K25, 60J74 Primary: Queues: Priority; secondary: Queues:
Markovian

Date created: September 05, 2023. Last update: December 07, 2023.

##  1 Introduction

This work is concerned with both algebraically exact and highly accurate
numerical evaluations of the joint and marginal queue length distributions for
the non-preemptive Markovian priority queue, based on a representation of the
joint queue-length probability mass function (PMF) as a single-variable
complex contour integral, making use of its pole and cut singularities on the
real axis.

There exist numerous real-world applications of the non-preemptive priority
queue, including telecommunications [5], health care [1, 15, 28], beam
scheduling in phased-array radars [21, 22], and air traffic control [23], to
name just a few. Here, we revisit the two-level problem from a different
perspective to that considered by us in [31] and the previous literature cited
therein. The two-level case is amenable to a number of exact analytical
approaches that are not easily accessible to the general multi-level problem.
Furthermore, the two-level problem occupies a distinguished position, because
knowledge of the low-priority marginal queue-length distribution for the two-
level case is sufficient for the evaluation of all the marginal queue-length
distributions in the general problem. For this reason, the two-level problem
justifies some special attention.

The two-level problem was recently discussed in [31], where a detailed survey
of previous work was also given. The most recent work on the problem, prior to
our own, appeared in [27]. The multi-level problem for an arbitrary number of
priority levels was solved in [32]. In that work, an explicit exact expression
was derived for the multi-variate joint queue-length probability generating
function (PGF), but the extraction of the joint and marginal queue-length PMFs
from the PGF was numerical, rather than analytical. Here, we describe an
analytical approach to the two-level problem. It is, at present, not clear
whether an analogous analytic method could be easily extended to the general
multi-level problem.

In the present approach, both the marginal and joint PMFs are represented as
one-dimensional real-valued integrals over a finite interval. The integral
representation is amenable to efficient numerical evaluation via Gaussian
quadrature in an iterative form that allows a desired level of accuracy to be
set a priori. It is equally amenable to an exact evaluation as a finite sum of
easily computable functions that generalize the familiar associated Legendre
functions in a novel way. Both methods produce highly accurate results for
parameter values spanning the entire ergodic region, including values of the
total traffic intensity close to unity.

The analytical structure of the integrals arising here is very similar to that
encountered for the partial busy period of an M/M/$c$ queue [17, 30], and for
the waiting-time marginal distributions in [6]. In fact, we are able to
provide alternative derivations of the waiting-time marginal distributions,
either from the queue-length marginal PGFs, or from the queue-length marginal
PMFs. In the latter case, we find a new representation as a convergent
Laguerre series, where the Laguerre polynomials carry the time dependence,
while the sequence of coefficients is a binomial transform of the queue-length
PMF.

For multi-class non-priority queueing models, Bertsimas and Mourtzinou [2]
derive a representation of the multi-variate joint PGF for the queue-length
distribution as a multi-dimensional time integral over all the waiting-time
marginals, based on the distributional form of Little’s law. However, this
approach does not extend to the priority queue because the underlying no-
overtaking assumption is violated. Consequently, in the present problem, the
joint queue-length distribution cannot be derived from the waiting-time
marginals. Here, and in [32], a direct derivation of the joint queue-length
distribution is provided.

The integral representations of the marginal and joint PMFs facilitate the
extraction of the asymptotic behaviour in the large queue-length limits;
something that appears to be difficult to achieve from other starting points.
The integral representations are also robust to the inclusion of various
modelling complications. An example arises in the application to the ambulance
ramping problem [1]. Here, the priority levels are determined by patient
acuity levels, but patients also fall into separate arrival classes depending
on whether they are entering the hospital emergency department directly from
an ambulance, from an ambulance offload zone or as walk-ins. Each priority
level may span a number of arrival classes, and each arrival class has a
distinct mean arrival rate. This results in a multi-level, multi-class
queueing model with priorities and classes entwined. The queue-length PMFs
acquire a somewhat more complex integral representation, but the simple
quadrature rules derived for the present model remain applicable and
efficient. We shall make use of this in a forthcoming paper.

The remainder of the paper is organized as follows: In Section 2, we introduce
an approach to the problem based on a PGF of a single complex variable that
sums the joint distribution only over the low-priority indices (i.e. queue
lengths), and we elucidate the analytic structure of the PGF in the complex
plane. In Section 3, we derive an integral representation for the low-priority
marginal queue-length distribution, noting that the high-priority marginal is
trivially geometric, and proceed to establish iterative quadrature rules. We
show that the results are related to one of Appell’s bivariate hypergeometric
functions. The asymptotic behaviour in the limit of large queue length is also
derived directly from the integral representation. In Section 4, we discuss
the marginal waiting-time distributions for the general multi-level non-
preemptive priority problem. As a by-product of our results on the queue-
length marginal distributions, we are able to invoke the distributional form
of Little’s law to reproduce a previously established integral form for the
waiting-time marginals, which may be related to a bivariate confluent
hypergeometric function of Humbert, as well as deriving a new representation
as a Laguerre series. In Section 5, we derive an integral representation for
the joint queue-length distribution, and establish the relevant iterative
quadrature rules. In Section 6, we first show that the low-priority marginal
PMF can be expressed as a finite sum of Legendre polynomials. For the joint
PMF, an exact representation as a finite sum from a family of functions
$Q_{m}^{\ell}(\chi)$ is derived. These functions reduce to the Legendre
polynomials in the special case $Q_{m}^{0}(\chi)=P_{m}(\chi)$, and are closely
related to the associated Legendre functions when $\absolutevalue{\ell}\leq
m$, but are defined for all integral values of $\ell$. Various properties of
these functions are elucidated. The section concludes with an efficient and
robust numerical implementation of the exact expressions. In Section 7, we
introduce a number of measures of performance (MOP) to study the numerical
behaviour of the two approaches developed in the present paper, testing them
against one other, against known exact relationships, and comparing with two
other algorithms previously derived for this problem [31]. Some concluding
remarks follow, and various technical details appear in the appendices.

##  2 Probability Generating Function

As in [31], we consider a non-preemptive priority queue with $K$ priority
levels, each with a distinct Poisson arrival rate $\lambda_{k}$,
$k=1,2,\ldots,K$ and corresponding level traffic intensity111Consistent with
[8, 16], the notation $\rho_{k}$ reserved for $\rho_{k}\equiv\lambda_{k}/\mu$,
so that $r_{k}=\rho_{k}/c$. $r_{k}=\lambda_{k}/(c\mu)$, leading to a total
traffic intensity for the aggregation of all arrivals of
$r=\sum_{k=1}^{K}r_{k}$. We adopt the usual convention that smaller priority-
level indices $k$ represent higher priorities. Thus, $r_{1}$ denotes the
traffic intensity associated with the highest priority level. We assume a
common exponential service rate $\mu$ among all priority levels. The number of
servers is arbitrary and is denoted by $c$. For the discussion of the joint
queue-distribution, we limit the number of priority levels to $K=2$, and adopt
the notation $r_{\text{hi}}\equiv r_{1}$, $r_{\text{lo}}\equiv r_{2}$, so that
$r=r_{\text{hi}}+r_{\text{lo}}$. We also denote the fraction of high-priority
arrivals (abbreviated ‘hifrac’) by $\nu\equiv r_{\text{hi}}/r$, so that
$0\leq\nu\leq 1$. This allows us to study all possibilities for a given total
traffic intensity $r$ by independently varying $\nu$ over its full range.

Let $P_{\text{full}}(\ell,m)$ denote the unconditional queue-length PMF with
$\ell$ high-priority clients and $m$ low-priority clients in the queue. Let
$P(\ell,m)$ denote the same queue-length PMF conditioned on all servers being
busy, to which we shall refer as the wait-conditional PMF. These two
quantities are related by

$P_{\text{full}}(\ell,m)=P_{\text{NW}}\delta_{\ell
0}\delta_{m0}+(1-P_{\text{NW}})P(\ell,m)\;,$ (1)

where $P_{\text{NW}}$ is the no-wait (or delay) probability. An expression for
its value is derived in Appendix 9. The wait-conditional PMF does not directly
depend on the number of servers, but only indirectly through its dependence on
the total traffic intensity. Therefore, the model for $P(\ell,m)$ depends on
just two parameters, namely, $r$ and $\nu$, that may be varied independently
in the intervals $0\leq r<1$, $0\leq\nu\leq 1$. Cohen [5] introduced a PGF
with a discrete index for the unconditional PMF given by

$G_{\ell}(z)\equiv\sum_{m=0}^{\infty}z^{m}P_{\text{full}}(\ell,m)\;.$ (2)

The corresponding wait-conditional PGF $g_{\ell}(z)$ is related to the full
PGF by

$G_{\ell}(z)=P_{\text{NW}}\delta_{\ell 0}+(1-P_{\text{NW}})g_{\ell}(z)\;.$ (3)

To obtain an expression for $g_{\ell}(z)$, we introduce
$\zeta\equiv\zeta_{\pm}(z)$ as the two branches that solve the $z$-dependent
quadratic equation

$\zeta^{2}-(1+r-r_{\text{lo}}z)\zeta+r_{\text{hi}}=0\;.$ (4)

Thus, we can write

$\zeta_{\pm}(z)=\tfrac{1}{2}\left[b(z)\pm\sqrt{\Delta(z)}\right]\;,$ (5)

where

$\Delta(z)\equiv b^{2}(z)-4r_{\text{hi}}\;,\quad b(z)\equiv
1+r-r_{\text{lo}}z\;,$ (6)

and it is useful to note the quadratic identities

$\zeta_{+}(z)+\zeta_{-}(z)=b(z)\;,\quad\zeta_{+}(z){\cdot}\zeta_{-}(z)=r_{\text{hi}}\;,$
(7)

which also imply that

$[1-\zeta_{+}(z)]{\cdot}[1-\zeta_{-}(z)]=r_{\text{lo}}(z-1)\;.$ (8)

Notable special values are $\zeta_{+}(1)=1$, $\zeta_{-}(1)=r_{\text{hi}}$.
Following Cohen [5], the wait-conditional PGF is given by

$g_{\ell}(z)=g_{\text{lo}}(z){\cdot}\left[1-\zeta_{-}(z)\right]\zeta_{-}^{\ell}(z)\;,\quad
g_{\text{lo}}(z)=\frac{1-r}{\zeta_{+}(z)-r}\;.$ (9)

The quantity $g_{\text{lo}}(z)$ represents the PGF for the low-priority
marginal distribution, as can be seen from

$\sum_{\ell=0}^{\infty}g_{\ell}(z)=g_{\text{lo}}(z)\;.$ (10)

These results also follow as a special case of the multi-variate joint PGF for
the general multi-level non-preemptive queue derived in [30].

The discriminant $\Delta(z)$ can be factorized according to

$\Delta(z)=r_{\text{lo}}^{2}(z-x_{+})(z-x_{-})\;,$ (11)

with real strictly positive roots

$x_{\pm}=1+(1\pm\sqrt{r_{\text{hi}}})^{2}/r_{\text{lo}}\;.$ (12)

It is useful to introduce their difference

$x_{\text{dif}}\equiv x_{+}-x_{-}=4\sqrt{r_{\text{hi}}}/r_{\text{lo}}\;.$ (13)

One should observe that $\Delta(x)<0$ for $x_{-}<x<x_{+}$, which implies the
presence of a cut on this finite real interval, and that
$\Delta(1/r)=(r^{2}-r_{\text{hi}})^{2}/r^{2}>0$. Thus, the pole at $x=1/r$
does not lie in the open cut interval $(x_{-},x_{+})$.

##  3 Marginal Distributions

The high and low priority marginal PMFs are given, respectively, by

$P_{\text{hi}}(\ell)=\sum_{m=0}^{\infty}P(\ell,m)\;,\quad
P_{\text{lo}}(m)=\sum_{\ell=0}^{\infty}P(\ell,m)\;.$ (14)

An immediate consequence of (9) is the simple geometric form for the high-
priority marginal

$P_{\text{hi}}(\ell)=g_{\ell}(1)=(1-r_{\text{hi}})r_{\text{hi}}^{\ell}\;.$
(15)

Thus, only the low-priority marginal $p_{n}\equiv P_{\text{lo}}(n)$ is of
interest for further analysis. By Cauchy’s integral theorem, we have

$p_{n}=\oint_{\mathcal{C}}\frac{dz}{2\pi
i}\,\frac{1}{z^{n+1}}g_{\text{lo}}(z)\;,$ (16)

where the integration contour $\mathcal{C}$ is taken to be an anti-clockwise
circle centred on the origin, of radius $\eta$ within the radius of
convergence of the function $g_{\text{lo}}(z)$, namely $0<\eta<1/r$. Our
approach is to evaluate this integral directly.

We begin by expanding the contour to infinity, while avoiding the pole and cut
by wrapping around them in the usual way. However, in order to deform the
circular contour out to infinity, we must be able to continue the integrand to
the entire $z$-plane in a way that it remains analytic. For any given
$z\in\mathbb{C}$, the solution of the corresponding quadratic equation has two
distinct branches, but it is a priori arbitrary how to assign the plus and
minus signs, that multiply the discriminant term, to the branches. The most
general form of $\zeta_{\pm}(z)$ is

$\zeta_{\pm}(z)=\tfrac{1}{2}\left[b(z)\pm\xi(z)\left(\Delta(z)\right)^{1/2}\right]\;,$
(17)

where $\xi(z)$ is an arbitrary function such that $\xi^{2}(z)\equiv 1$, and
the principal branch of the square-root function is adopted. Analyticity of
$\zeta_{\pm}(z)$ everywhere in the complex $z$-plane away from the positive
real axis is achieved by the choice

$\xi(z)=\left\\{\begin{array}[]{ccc}+1&\quad\text{for}&\mathop{\rm
Re}\nolimits z\leq(x_{+}+x_{-})/2\\\ -1&\quad\text{for}&\mathop{\rm
Re}\nolimits z>(x_{+}+x_{-})/2\end{array}\right.\;.$ (18)

With this definition, we observe that

$\zeta_{-}(1/r)=\left\\{\begin{array}[]{ccc}r&\quad\text{for}&r_{\text{hi}}>r^{2}\\\
r_{\text{hi}}/r&\quad\text{for}&r_{\text{hi}}\leq r^{2}\end{array}\right.\;,$
(19)

which implies that there is no pole when $r_{\text{hi}}\geq r^{2}$. The
behaviour of $\zeta_{+}(1/r)$ is the opposite and follows directly from (7).
With the choice in (18), we also obtain the uniform asymptotic behaviour

$\zeta_{+}(z)\mathrel{\mathop{\sim}_{|z|\to\infty}}-r_{\text{lo}}z\;,\quad\zeta_{-}(z)\mathrel{\mathop{\sim}_{|z|\to\infty}}-\frac{r_{\text{hi}}}{r_{\text{lo}}z}\;,$
(20)

which is sufficient to ensure the vanishing of the contour at infinity.

After deformation of the contour, $p_{n}$ may be expressed as the sum of a
pole contribution and a cut integral
$p_{n}=P_{\text{pol}}(n)+P_{\text{cut}}(n)$, where

$P_{\text{pol}}(n)=\Theta(r^{2}-r_{\text{hi}})\left[1-\frac{r(1-r)}{r_{\text{lo}}}\right]{\cdot}(1-r)r^{n-1}\;,$
(21)

and

$P_{\text{cut}}(n)=\frac{1-r}{2\pi
r}\int_{x_{-}}^{x_{+}}\frac{dx}{x^{n+1}}\,\frac{\sqrt{(x_{+}-x)(x-x_{-})}}{x-1/r}\;.$
(22)

The change of integration variable $x\mapsto u$, $0\leq u\leq 1$, according to

$x=x_{-}+(x_{+}-x_{-})u\;,$ (23)

yields

$P_{\text{cut}}(n)=\frac{1-r}{2\pi
rx_{\text{dif}}^{n}}\int_{0}^{1}du\,\frac{\sqrt{u(1-u)}}{(u+a)^{n+1}(u+b)}\;,$
(24)

with

$a\equiv x_{-}/x_{\text{dif}}\;,\quad
b\equiv(x_{-}-1/r)/x_{\text{dif}}\;,\quad x_{\text{dif}}\equiv x_{+}-x_{-}\;,$
(25)

so that $1/x_{\text{dif}}=r(a-b)$. These parameters may also be conveniently
expressed in terms of the model inputs as

$a=\frac{1+r}{4\sqrt{r_{\text{hi}}}}-\frac{1}{2}\;,\quad
b=\frac{1}{4}\left(\frac{r}{\sqrt{r_{\text{hi}}}}+\frac{\sqrt{r_{\text{hi}}}}{r}\right)-\frac{1}{2}\;.$
(26)

In terms of Appell’s $F_{1}$ bivariate hypergeometric function, we may write

$\displaystyle\begin{aligned}
P_{\text{cut}}(n)&=\frac{1-r}{16rx_{\text{dif}}^{n}}{\cdot}\frac{1}{a^{n+1}b}{\cdot}\frac{1}{B(3/2,3/2)}\int_{0}^{1}du\,\frac{\sqrt{u(1-u)}}{\left(1+u/a\right)^{n+1}\left(1+u/b\right)}\\\
&=\frac{(1-r)r^{n-1}}{16ab}{\cdot}\left(1-b/a\right)^{n}F_{1}(3/2,n+1,1,3;-1/a,-1/b)\;,\end{aligned}$
(27)

where it is useful to note that $1-b/a=1/(rx_{-})$, and with $B(\mu,\nu)$
denoting the beta function. This follows from Picard’s integral representation
[24]

$F_{1}(\lambda,\rho,\sigma,\lambda+\mu;u,v)=\frac{1}{B(\mu,\lambda)}\int_{0}^{1}dx\,\frac{x^{\lambda-1}(1-x)^{\mu-1}}{(1-ux)^{\rho}(1-vx)^{\sigma}}\;,$
(28)

as reproduced in [13, sections 3.211, 9.18]. When, as in the present
application, the arguments $u,v$ are real and negative, the Appell
$F_{1}$-hypergeometric function can be efficiently evaluated by means of a
Gauss-Jacobi quadrature rule [10, 11]. Furthermore, in the present
application, where $(\mu,\lambda)=(3/2,3/2)$, the Gauss-Jacobi quadrature
specializes to a Gauss-Chebyshev quadrature of the second kind. To see this,
we make the further change of integration variable $v=2u-1$, which casts the
integral into the form

$P_{\text{cut}}(n)=C_{n}\int_{-1}^{+1}\frac{dv}{\pi/2}\,\frac{\sqrt{1-v^{2}}}{(v+\alpha)^{n+1}(v+\beta)}\;,\quad
C_{n}\equiv\frac{1-r}{4r}{\cdot}\left(\frac{2}{x_{\text{dif}}}\right)^{n}\;,$
(29)

where we have introduced the parameters

$\alpha\equiv 2a+1\;,\quad\beta\equiv 2b+1\;.$ (30)

Thus,

$\alpha=\frac{1+r}{2\sqrt{r_{\text{hi}}}}\;,\quad\beta=\frac{1}{2}\left(\frac{r}{\sqrt{r_{\text{hi}}}}+\frac{\sqrt{r_{\text{hi}}}}{r}\right)\;.$
(31)

It is known that Gauss-Chebyshev quadrature of the second kind is equivalent
to the trapezoidal rule on the unit circle [4]. Accordingly, if we set
$v=\cos(\pi\tau)$, so that $0\leq\tau\leq 1$, and subsequently uniformly
discretize the unit interval into $L$ sub-intervals each terminated at the
point $\tau_{k}=k/L$, $k=1,2,\ldots,L$, then we obtain the $L$-point
quadrature rule for the cut integral given by

$P_{\text{cut}}(n)\simeq
C_{n}\sum_{k=1}^{L}\frac{w_{k}}{[\cos(\pi\tau_{k})+\alpha]^{n+1}{\cdot}[\cos(\pi\tau_{k})+\beta]}\;,\quad
w_{k}\equiv(2/L){\cdot}\sin^{2}(\pi\tau_{k})\;.$ (32)

This structure mirrors the discussion of the distribution of the partial busy
period for a multi-server Markovian queue found in [32]. An alternative way of
expressing this, along the lines of [32], is

$P_{\text{cut}}(n)\simeq\frac{x_{\text{dif}}}{4}{\cdot}\sum_{k=1}^{L}\frac{W_{k}}{\left[x_{\text{dif}}(U_{k}+a)\right]^{n+1}}\;,$
(33)

with

$U_{k}\equiv\cos^{2}(\pi\tau_{k}/2)\;,\quad
W_{k}\equiv\frac{2(1-r)}{rL}{\cdot}\frac{1-U_{k}}{1+b/U_{k}}\;.$ (34)

The case $\beta=1$, equivalently $b=0$, is problematic for the computation of
the integral in (29) via the Gauss-Chebyshev quadrature of the second kind as
the denominator vanishes at an integration end-point, where the weight
function also vanishes. Thus, a removable singularity appears. At this
critical point, which corresponds to parameter choice $r_{\text{hi}}=r^{2}$ —
the same point at which the pole contribution disappears — (29) may be cast
into a form suitable for Gauss-Chebyshev quadrature of the first kind:

$P_{\text{cut}}(n)=C_{n}\int_{-1}^{+1}\frac{dv}{\pi/2}\,\frac{1}{\sqrt{1-v^{2}}}{\cdot}\frac{1-v}{(v+\alpha)^{n+1}}\;.$
(35)

The loss of accuracy away from the critical point but close to it may be
addressed by noting that there is no impediment to evaluating (29) by Gauss
quadrature of the first kind in the general case, noting, in particular, that
it does not place a node at $v=-1$. Furthermore, it is known that Gauss-
Chebyshev quadrature of the first kind is equivalent to the mid-point rule on
the unit circle. In fact, the quadrature rule expressed by (33) and (34) holds
equally for the midpoint rule, provided that one now chooses
$\tau_{k}=(k-1/2)/L$ for $k=1,2,\ldots,L$. The advantage of the associations
with the trapezoidal and mid-point rules is that they can both be applied
iteratively, in order to attain a pre-assigned level of accuracy [26].
Generally, in a Gaussian quadrature scheme, the nodes and weights are
generated for a pre-selected quadrature order $L$, and it is difficult to
estimate a priori the accuracy that will result for any particular problem.
While one can successively increase the order to refine the integration grid,
thus achieving greater accuracy, the nodes from the previous iteration cannot
be re-used. The mapping of a Gaussian quadrature to the trapezoidal rule
allows node re-use if the grid is doubled on each refinement and, for the mid-
point rule, node re-use results from tripling the grid on each refinement
[25]. While this makes the trapezoidal rule more efficient than the mid-point
rule in reaching a desired level of convergence, in the present problem, a
price is paid for the singular behaviour around the critical point. It should
be appreciated that the trapezoidal and mid-point rules do not generally
implement Gaussian quadratures. The fact that it is so in the present case
leads to an exponential rate of convergence [29].

Finally, we may note that, in the special case $r_{\text{hi}}=r^{2}$, we have

$p_{n}=\tensor[_{2}]{F}{{}_{1}}\left(n+1,1/2;2;-4r/(1-r)\right){\cdot}r^{n}\;,$
(36)

recalling that the pole contribution vanishes here, and using the result that
the Gauss hypergeometric function has integral representation

$\tensor[_{2}]{F}{{}_{1}}(a,b;c;z)=\frac{1}{B(b,c-b)}\int_{0}^{1}dt\,\frac{t^{b-1}(1-t)^{c-b-1}}{(1-zt)^{a}}\;.$
(37)

###  3.1 Asymptotic Limits

Let us consider the integral

$\mathcal{I}_{n}(a,b)\equiv\int_{0}^{1}du\,\frac{\sqrt{u(1-u)}}{(u+a)^{n+1}(u+b)}\mathrel{\mathop{\sim}_{n\to\infty}}\frac{1}{b}\int_{0}^{1}du\,\frac{\sqrt{u}}{(u+a)^{n+1}}\;,$
(38)

and continue the limiting procedure with

$\displaystyle\begin{aligned}
\int_{0}^{1}du\,\frac{\sqrt{u}}{(u+a)^{n+1}}&\;\;=\;\;\frac{2}{a^{n-1/2}}\int_{0}^{a^{2}}dv\,\frac{v^{2}}{(1+v^{2})^{n+1}}\\\
&\mathrel{\mathop{\sim}_{n\to\infty}}\frac{2}{a^{n-1/2}}\int_{0}^{\infty}dv\,\frac{v^{2}}{(1+v^{2})^{n+1}}\\\
&\;\;=\;\;\frac{\sqrt{\pi}}{2a^{n-1/2}}{\cdot}\frac{\Gamma(n-1/2)}{\Gamma(n+1)}\;.\end{aligned}$
(39)

Hence, using the result that
$\Gamma(n+c)\mathrel{\mathop{\sim}_{n\to\infty}}\Gamma(n)n^{c}$, we obtain

$\mathcal{I}_{n}(a,b)\mathrel{\mathop{\sim}_{n\to\infty}}\frac{\sqrt{\pi}}{2ba^{n-1/2}}{\cdot}\frac{1}{n^{3/2}}\;,$
(40)

and so we conclude that

$P_{\text{cut}}(n)\mathrel{\mathop{\sim}_{n\to\infty}}\frac{1-r}{4rb}\sqrt{\frac{a}{\pi
n^{3}}}{\cdot}\left[(1-b/a)r\right]^{n}\;,$ (41)

provided $b\neq 0$. The situation where $b=0$ occurs when
$r_{\text{hi}}=r^{2}$. In this case, we must consider

$\mathcal{I}_{n}(a)\equiv\int_{0}^{1}du\,\frac{1}{(u+a)^{n+1}}\sqrt{\frac{1-u}{u}}\mathrel{\mathop{\sim}_{n\to\infty}}\int_{0}^{1}\frac{du}{\sqrt{u}}\,\frac{1}{(u+a)^{n+1}}\;.$
(42)

Now,

$\displaystyle\begin{aligned}
\int_{0}^{1}\frac{du}{\sqrt{u}}\,\frac{1}{(u+a)^{n+1}}&\;\;=\;\;\frac{2}{a^{n+1/2}}\int_{0}^{a^{2}}dv\,\frac{1}{(1+v^{2})^{n+1}}\\\
&\mathrel{\mathop{\sim}_{n\to\infty}}\frac{2}{a^{n+1/2}}\int_{0}^{\infty}dv\,\frac{1}{(1+v^{2})^{n+1}}\\\
&\;\;=\;\;\frac{\sqrt{\pi}}{a^{n+1/2}}{\cdot}\frac{\Gamma(n+1/2)}{\Gamma(n+1)}\;,\end{aligned}$
(43)

and so it follows that

$\mathcal{I}_{n}(a)\mathrel{\mathop{\sim}_{n\to\infty}}\frac{1}{a^{n+1/2}}\sqrt{\frac{\pi}{n}}\;.$
(44)

Therefore, when $r_{\text{hi}}=r^{2}$, in which case $a=(1-r)/(4r)$, we obtain

$P_{\text{cut}}(n)\mathrel{\mathop{\sim}_{n\to\infty}}\sqrt{\frac{1-r}{\pi
r}}{\cdot}\frac{r^{n}}{\sqrt{n}}\;.$ (45)

##  4 Waiting-Time Distribution

It is known that the queue-length marginals can be recovered from the waiting-
time marginals by appealing to the distributional form of Little’s law [3,
18], bearing in mind the caveat that it is applicable to priority queues only
within a given priority level. Let us denote the Poisson PMF for arrival rate
$\lambda$ and time $t$ by

$\Lambda_{n}(\lambda t)\equiv\frac{(\lambda t)^{n}}{n!}e^{-\lambda t}\;,$ (46)

for arrival numbers $n=0,1,2,\ldots$. Then, one has

$P_{\mathscr{L}}(n)=\int_{0}^{\infty}dt\,P_{\mathscr{W}}(t)\Lambda_{n}(r_{\text{lo}}t)\;,$
(47)

where $\mathscr{L},\mathscr{W}$ denote the queue-length and waiting-time
random variables (RVs), respectively. This can be inverted, hence enabling us
to derive the waiting-time probability density function (PDF)
$P_{\mathscr{W}}(t)$ as a by-product of our derivation of the queue-length
marginal PMF $P_{\mathscr{L}}(n)$. It is implicit in (47) that we are dealing
with the low-priority level and that we have, without loss of generality, set
$c\mu=1$, which means that $\lambda_{\text{lo}}=r_{\text{lo}}$. This simply
implies that time is being measured in units of $c\mu$.

Let us consider

$g_{\text{lo}}(z)=\sum_{n=0}^{\infty}z^{n}P_{\mathscr{L}}(n)=\int_{0}^{\infty}dt\,P_{\mathscr{W}}(t)K(z,r_{\text{lo}}t)\;,$
(48)

where $K(z,\lambda t)$ is the PGF of the Poisson distribution

$K(z,\lambda t)\equiv\sum_{n=0}^{\infty}z^{n}\Lambda_{n}(\lambda t)=e^{\lambda
t(z-1)}\;.$ (49)

It follows immediately that the moment generating function (MGF) of the
waiting-time RV,

$M_{\mathscr{W}}(z)\equiv\langle
e^{z\mathscr{W}}\rangle_{\mathscr{W}}=\int_{0}^{\infty}dt\,P_{\mathscr{W}}(t)e^{zt}\;,$
(50)

is given by $M_{\mathscr{W}}(z)=g_{\text{lo}}(1+z/r_{\text{lo}})$. The
waiting-time PDF can now be recovered via the inverse Laplace transform as

$P_{\mathscr{W}}(t)=\int_{\Gamma}\frac{dz}{2\pi
i}\,e^{-tz}g_{\text{lo}}(1+z/r_{\text{lo}})\;,$ (51)

with Bromwich contour $\Gamma$ such that $z=\gamma+iy$, for
$-\infty<y<+\infty$ with any given real constant $\gamma$ such that
$0<\gamma<(1/r-1)r_{\text{lo}}$. After a scaling and shift of the integration
variable $z$, and a subsequent deformation of the resulting Bromwich contour
so that it wraps around the pole and cut of the function $g_{\text{lo}}(z)$ on
the real axis, we arrive at the representation

$P_{\mathscr{W}}(t)dt=r_{\text{lo}}dt{\cdot}\int_{\mathcal{C}_{\text{pol}}\cup\mathcal{C}_{\text{cut}}}\frac{dz}{2\pi
i}\,e^{-r_{\text{lo}}t(z-1)}g_{\text{lo}}(z)\;,$ (52)

noting that $\mathop{\rm Re}\nolimits z>1$ everywhere on
$\mathcal{C}=\mathcal{C}_{\text{pol}}\cup\mathcal{C}_{\text{cut}}$. An
alternative derivation is presented in Appendix 10. While more complex, it has
the merit of establishing a novel representation for the waiting-time PDF in
terms of the queue-length PMF as a convergent series of Laguerre polynomials
whose coefficients are finite linear combinations of queue-length
probabilities, corresponding to a binomial transform.

Let us consider the closed contour integral, for an arbitrary analytic
function $f(z)$, given by

$\mathcal{I}[f]\equiv\oint_{\mathcal{C}}\frac{dz}{2\pi
i}\,f(z)g_{\text{lo}}(z)\;,$ (53)

with

$g_{\text{lo}}(z)=\frac{1-r}{\zeta_{+}(z)-r}=\frac{1}{r_{\text{lo}}r}{\cdot}\frac{1-r}{z-1/r}{\cdot}\left[\zeta_{-}(z)-r\right]\;,$
(54)

where we recall that

$\lim_{z\to
1/r}\left(z-1/r\right)g_{\text{lo}}(z)=\frac{1-r}{r^{2}}\left[\frac{r(1-r)}{r_{\text{lo}}}-1\right]\Theta(r^{2}-r_{\text{hi}})\;.$
(55)

The chosen contour deformation leads to the pole/cut decomposition

$\mathcal{I}[f]=\int_{\mathcal{C}_{\text{pol}}\cup\mathcal{C}_{\text{cut}}}\frac{dz}{2\pi
i}\,f(z)g_{\text{lo}}(z)=\mathcal{I}_{\text{pol}}[f]+\mathcal{I}_{\text{cut}}[f]\;,$
(56)

with

$\mathcal{I}_{\text{pol}}[f]=\frac{1-r}{r^{2}}\left[1-\frac{r(1-r)}{r_{\text{lo}}}\right]\Theta(r^{2}-r_{\text{hi}}){\cdot}f(1/r)\;,$
(57)

and

$\displaystyle\begin{aligned} \mathcal{I}_{\text{cut}}[f]&=\frac{1-r}{2\pi
r}\int_{x_{-}}^{x_{+}}dx\,\frac{\sqrt{(x_{+}-x)(x-x_{-})}}{x-1/r}f(x)\\\
&=\frac{(1-r)x_{\text{dif}}}{2\pi
r}\int_{0}^{1}du\,\frac{\sqrt{u(1-u)}}{u+b}f(x_{\text{dif}}(u+a))\;.\end{aligned}$
(58)

It should be noted that $\mathcal{C}_{\text{pol}}$ and
$\mathcal{C}_{\text{cut}}$ are both clockwise contours about the pole at
$z=1/r$ and cut along $z\in[x_{-},x_{+}]$, respectively. Application to the
waiting-time PDF entails considering

$f(z;t)=r_{\text{lo}}e^{-r_{\text{lo}}t(z-1)}\;,$ (59)

in which case the wait-conditional waiting-time PDF, given by
$P_{W}(t)=\mathcal{I}[f](t)$, becomes

$\displaystyle\begin{aligned}
P_{\mathscr{W}}(t)&=-(1-r)\left(\frac{r_{\text{hi}}}{r^{2}}-1\right)\Theta(r^{2}-r_{\text{hi}})e^{-tr_{\text{lo}}(1/r-1)}\\\
&\quad{}+\frac{2(1-r)\sqrt{r_{\text{hi}}}}{\pi
r}e^{-t(4\sqrt{r_{\text{hi}}}a-r_{\text{lo}})}\int_{0}^{1}du\,e^{-4t\sqrt{r_{\text{hi}}}{\cdot}u}{\cdot}\frac{\sqrt{u(1-u)}}{u+b}\;,\end{aligned}$
(60)

where we have observed that
$r_{\text{lo}}x_{\text{dif}}=4\sqrt{r_{\text{hi}}}$. The integral appearing
here is an Euler-type integral representation for a special case of the
Humbert series $\Phi_{1}$, given by [7] as

$\Phi_{1}(a,b,c;x,y)=\frac{1}{B(a,c-a)}\int_{0}^{1}dt\,e^{yt}{\cdot}\frac{t^{a-1}(1-t)^{c-a-1}}{(1-xt)^{b}}\;.$
(61)

The Humbert series, also known as confluent Appell functions, are bivariate
generalizations of Kummer’s confluent hypergeometric function
$\tensor[_{1}]{F}{{}_{1}}(a;c;y)$. In the present case, we have
$\tensor[_{1}]{F}{{}_{1}}(a;c;y)=\Phi_{1}(a,b,c;0,y)$ for any $b$. The
function $\Phi_{1}$ in the general case is efficiently computed via Gauss-
Jacobi quadrature [9]. In the present specific application, the iterative
Gauss-Chebyshev scheme described earlier can be utilized effectively.

It is straightforward to show explicitly that $P_{\mathscr{W}}(0)=1-r$. A
useful intermediate result is

$\left(\sqrt{b+1}-\sqrt{b}\right)^{2}=\left\\{\begin{array}[]{ccc}\sqrt{r_{\text{hi}}}/r&\quad\text{for}&r_{\text{hi}}\leq
r^{2}\;,\\\
r/\sqrt{r_{\text{hi}}}&\quad\text{for}&r_{\text{hi}}>r^{2}\;.\end{array}\right.$
(62)

By contrast, it requires considerable algebra to confirm explicitly that
$P_{\mathscr{W}}(t)$ integrates to unity; whereas it is easy to see this on
general grounds. Taking the $t$-integral inside the $z$-integral yields

$\int_{0}^{\infty}dt\,f(z;t)=\frac{1}{z-1}\;,$ (63)

provided $\mathop{\rm Re}\nolimits z>1$. Then, on deforming the integration
contour $\mathcal{C}\mapsto\mathcal{C}_{\eta}$, with the contour
$\mathcal{C}_{\eta}$ denoting the anti-clockwise circle around the origin of
radius $1<\eta<1/r$, we obtain

$\int_{0}^{\infty}dt\,P_{\mathscr{W}}(t)=\oint_{\mathcal{C}_{\eta}}\frac{dz}{2\pi
i}\,\frac{g_{\text{lo}}(z)}{z-1}=g_{\text{lo}}(1)=1\;.$ (64)

The wait-conditional low-priority waiting-time marginal PDF given by (60)
agrees with the result in [6] derived from a different starting point, based
on a first-principles derivation of the MGF. An alternative derivation of the
same MGF has been presented in [19]. This result can be used to represent
waiting-time marginals for all priority levels in the general multi-level
problem. Let $r_{\kappa}$, $\kappa=1,2,\ldots,K$ denote the level traffic
intensities for a $K$-level problem, ordered in descending priority, so that
$r_{\text{agg}}=r_{1}+r_{2}+\cdots+r_{K}<1$ is the total (aggregated) traffic
intensity. A common mean service rate $\mu$ is assumed. We also introduce the
cumulative quantities [14]

$\sigma_{\kappa}\equiv\sum_{j=1}^{\kappa}r_{j}\;,$ (65)

so that $\sigma_{K}=r_{\text{agg}}$, and make the associations with an
effective two-level problem according to

$r_{\text{lo}}=r_{\kappa}\;,\quad r_{\text{hi}}=\sigma_{\kappa-1}\;,\quad
r=r_{\text{lo}}+r_{\text{hi}}\;,$ (66)

Then $P_{\mathscr{W}_{\kappa}}(t)$, as given by (60), yields the wait-
conditional waiting-time PDF for the $\kappa$-th highest priority level. In
particular, for the highest priority level ($\kappa=1$), we would set
$r_{\text{lo}}=r_{1}$, $r_{\text{hi}}=0$, so that $r=r_{1}$, which leads to
the result, arising solely from the pole contribution, that

$P_{\mathscr{W}_{1}}(t)=(1-r_{1})e^{-(1-r_{1})t}\;,$ (67)

being an exponential distribution, as expected. Since the waiting-time
marginal PDFs have already appeared in the literature in a similar form [6],
we shall not discuss them further in the present work. As far as we are aware,
the only explicit representation of the waiting-time distributions that exist
in the literature is the integral form first derived in [6] and elaborated
upon here. In Appendix 10, we establish a different representation as a
Laguerre series.

##  5 Joint Distribution

For the wait-conditional joint queue-length distribution, the pole/cut
decomposition reads

$P(\ell,m)=P_{\text{pol}}(\ell,m)+P_{\text{cut}}(\ell,m)\;,$ (68)

with

$\displaystyle\begin{aligned} P_{\text{pol}}(\ell,m)&=-\lim_{z\to
1/r}\left(z-1/r\right)g_{\ell}(z)/z^{m+1}\;,\\\
P_{\text{cut}}(\ell,m)&=\frac{1}{2\pi
i}\int_{x_{-}}^{x_{+}}\frac{dx}{x^{m+1}}\,\left[g_{\ell}(x+i\epsilon)-g_{\ell}(x-i\epsilon)\right]\;,\end{aligned}$
(69)

where $\epsilon$ is a positive infinitesimal. On inserting the explicit form
of the PGF, the cut contribution reads

$P_{\text{cut}}(\ell,m)=\frac{1-r}{2\pi
i}\int_{x_{-}}^{x_{+}}dx\,\frac{i\sqrt{|\Delta(x)|}}{x^{m+1}(1-rx)}\chi_{\ell}(x)-\frac{1-r}{2\pi
i}\int_{x_{-}}^{x_{+}}dx\,\frac{i\sqrt{|\Delta(x)|}}{x^{m}(1-rx)}\chi_{\ell+1}(x)\;,$
(70)

where

$\chi_{\ell}(x)\equiv-\frac{\zeta^{\prime\ell}_{-}(x)-\zeta^{\prime\ell}_{+}(x)}{i\sqrt{|\Delta(x)|}}\;,\quad\zeta^{\prime}_{\pm}(x)\equiv\lim_{\epsilon\to
0^{+}}\zeta_{\pm}(x-i\epsilon)\;,$ (71)

is an analytic function. On the cut, we can write

$\displaystyle\begin{aligned}
\zeta_{\pm}(x)&=\frac{r_{\text{lo}}}{4}\left(\sqrt{x_{+}-x}\pm
i\sqrt{x-x_{-}}\right)^{2}\\\ &=\sqrt{r_{\text{hi}}}\left(\sqrt{1-u}\pm
i\sqrt{u}\right)^{2}\;.\end{aligned}$ (72)

We also have

$|\Delta(x)|=r_{\text{lo}}^{2}x_{\text{dif}}^{2}u(1-u)=16r_{\text{hi}}u(1-u)\;.$
(73)

This leads to the result

$\displaystyle\begin{aligned}
\chi_{\ell}(x(u))&=(-1)^{\ell-1}r_{\text{hi}}^{(\ell-1)/2}{\cdot}\frac{\left(y+\sqrt{y^{2}-1}\right)^{\ell}-\left(y-\sqrt{y^{2}-1}\right)^{\ell}}{2\sqrt{y^{2}-1}}\\\
&=(-1)^{\ell-1}r_{\text{hi}}^{(\ell-1)/2}{\cdot}U_{\ell-1}(y)\\\
&=(-1)^{\ell-1}r_{\text{hi}}^{(\ell-1)/2}{\cdot}U_{\ell-1}(2u-1)\\\
&=r_{\text{hi}}^{(\ell-1)/2}{\cdot}\frac{\sin(2\ell\asin\sqrt{u})}{\sin(2\asin\sqrt{u})}\;,\end{aligned}$
(74)

where we had temporarily set $y=2u-1$ and used some standard representations
for the Chebyshev polynomial $U_{\ell}(y)$ of the second kind of order $\ell$.
One may observe that $\chi_{0}(x)=0$, $\chi_{1}(x)=1$. From the Chebyshev
generating function

$\sum_{\ell=0}^{\infty}t^{\ell}U_{\ell}(y)=\frac{1}{1-2ty+t^{2}}\;,$ (75)

one obtains

$\sum_{\ell=0}^{\infty}\chi_{\ell}(x(u))=\frac{1}{1+2\sqrt{r_{\text{hi}}}(2u-1)+r_{\text{hi}}}\;,$
(76)

which is useful for verifying that the low-priority marginal is recovered when
one sums over the high-priority indices $\ell$.

After some algebraic manipulation, we arrive at the expression

$\begin{split}P_{\text{cut}}(\ell,m)=&-C_{1}\int_{0}^{1}\frac{du}{\pi}{\cdot}\frac{\sqrt{u(1-u)}}{(u+a)^{m+1}(u+b)}{\cdot}\frac{\sin(2\ell\asin\sqrt{u})}{\sin(2\asin\sqrt{u})}\\\
&{}+C_{2}\int_{0}^{1}\frac{du}{\pi}{\cdot}\frac{\sqrt{u(1-u)}}{(u+a)^{m}(u+b)}{\cdot}\frac{\sin(2(\ell+1)\asin\sqrt{u})}{\sin(2\asin\sqrt{u})}\;,\end{split}$
(77)

with coefficients

$C_{1}\equiv\frac{2(1-r)}{r}\left(\frac{r_{\text{lo}}}{4}\right)^{m+1}{\cdot}r_{\text{hi}}^{(\ell-m-1)/2}\;,\quad
C_{2}\equiv\frac{2(1-r)}{r}\left(\frac{r_{\text{lo}}}{4}\right)^{m}{\cdot}r_{\text{hi}}^{(\ell-m+1)/2}\;.$
(78)

As for the marginal PMF, we set $u=(1+v)/2$, $v=\cos(\pi\tau)$, $0\leq\tau\leq
1$, so that

$u(\tau)=\left[1+\cos(\pi\tau)\right]/2=\cos^{2}(\pi\tau/2)\;,$ (79)

and

$\frac{du}{\pi}{\cdot}\frac{\sqrt{u(1-u)}}{u}=-d\tau{\cdot}[1-u(\tau)]\;.$
(80)

We also have

$C_{\ell}(u)\equiv\frac{\sin(2\ell\asin\sqrt{u})}{\sin(2\asin\sqrt{u})}=(-1)^{\ell-1}\frac{\sin(\ell\pi\tau)}{\sin(\pi\tau)}=(-1)^{\ell-1}\ell{\cdot}\frac{\mathop{\rm
sinc}\nolimits(\ell\pi\tau)}{\mathop{\rm sinc}\nolimits(\pi\tau)}\;,$ (81)

and we can write $P_{\text{cut}}(\ell,m)=\Lambda(\ell+1,m)-\Lambda(\ell,m+1)$,
where

$\Lambda(\ell,m)\equiv\frac{2(1-r)}{r}{\cdot}\frac{r_{\text{hi}}^{\ell/2}}{x_{-}^{m}}\int_{0}^{1}d\tau\,\frac{[1-u(\tau)]C_{\ell}(u(\tau))}{[1+u(\tau)/a]^{m}[1+b/u(\tau)]}\;.$
(82)

The $L$-point quadrature with nodes $\tau_{k}$, $k=1,2,\ldots,L$,
corresponding to the trapezoidal rule is obtained by introducing the
integration grid $\tau_{k}=k/L$, while the mid-point rule corresponds to
$\tau_{k}=(k-1/2)/L$. In both cases, the resulting quadrature rule may be
expressed as

$\Lambda(\ell,m)\simeq
r_{\text{hi}}^{\ell/2}\sum_{k=1}^{L}\frac{W_{k}C_{\ell}(U_{k})}{[x_{-}(1+U_{k}/a)]^{m}}\;,$
(83)

where the quadrature nodes and weights $U_{k}$, $W_{k}$, respectively, are
given by (34). As a sanity check, we may note that the identity for the
exclusively-low probabilities
$P_{\text{cut}}(0,m+1)=r_{\text{lo}}P_{\text{cut}}(m)$ implies that
$\Lambda(1,m+1)=r_{\text{lo}}P_{\text{cut}}(m)$. Recalling that
$x_{-}=ax_{\text{dif}}$, it is straightforward to verify that the quadrature
rules (83) and (33) applied to the left-hand and right-hand sides,
respectively, reproduce this equality. The explicit result for the pole
contribution is easily obtained as

$P_{\text{pol}}(\ell,m)=\Theta(r^{2}-r_{\text{hi}}){\cdot}(1-r)r^{m}\left(1-\frac{r_{\text{hi}}}{r^{2}}\right)\left(\frac{r_{\text{hi}}}{r}\right)^{\ell}\;.$
(84)

##  6 Exact Expressions

###  6.1 Marginal Distribution

In order to facilitate an exact representation of the low-priority marginal
PMF, we introduce a family of integrals to which we shall refer as the
marginal integrals. These are given by

$M_{n}(\alpha,\beta)\equiv\int_{0}^{\pi}d\theta\,\frac{1-\cos^{2}\theta}{(\alpha+\cos\theta)^{n+1}(\beta+\cos\theta)}\;,$
(85)

for $n=-1,0,\ldots$, so that
$M_{-1}(\alpha,\beta)=\pi(\beta-\sqrt{\beta^{2}-1})$. It is assumed that
$|\alpha|,|\beta|>1$. In terms of the marginal integral, we have

$P_{\text{cut}}(n)=\frac{(1-r)r^{n-1}}{2\pi}{\cdot}(\alpha-\beta)^{n}M_{n}(\alpha,\beta)\;,$
(86)

with the parameters $\alpha,\beta$ defined in (30). The marginal integral may
be evaluated in terms of Legendre polynomials $P_{n}(\chi)$ due to the
identity [13, section 3.661.4]

$\frac{1}{\pi}\int_{0}^{\pi}\frac{d\theta}{(\alpha+\cos\theta)^{n+1}}=\frac{1}{(\alpha^{2}-1)^{(n+1)/2}}P_{n}\left(\frac{\alpha}{\sqrt{\alpha^{2}-1}}\right)\;.$
(87)

for $\alpha>1$ and $n=-1,0,1,\ldots$, where we formally set
$P_{-1}(\chi)\equiv 1$. After a partial fraction decomposition of the
integrand, one obtains for $n=0,1,\ldots$,

$M_{n}(\alpha,\beta)=\frac{\pi}{(\alpha-\beta)^{n+1}}\sum_{k=-1}^{n}\psi_{k}(n)P_{k}\left(\frac{\alpha}{\sqrt{\alpha^{2}-1}}\right)\;,$
(88)

with

$\psi_{k}(n)=\frac{(\alpha-\beta)^{k}}{(\alpha^{2}-1)^{(k+1)/2}}\left[\beta^{2}-1+(\alpha^{2}-\beta^{2})\delta_{k,n}-(\alpha-\beta)^{2}\delta_{k,n-1}\right]\;,$
(89)

for $k=1,2,\ldots,n$, and

$\displaystyle\begin{aligned} \psi_{-1}(n)&=(\beta-\alpha)\delta_{n,0}\;,\\\
\psi_{0}(n)&=(\beta^{2}-1)\left(\frac{1}{\sqrt{\alpha^{2}-1}}-\frac{1}{\sqrt{\beta^{2}-1}}\right)+\frac{\alpha^{2}-\beta^{2}}{\sqrt{\alpha^{2}-1}}\delta_{n,0}-\frac{(\alpha-\beta)^{2}}{\sqrt{\alpha^{2}-1}}\delta_{n,1}\;.\end{aligned}$
(90)

An immediate consequence is that

$M_{0}(\alpha,\beta)/\pi=\frac{\psi_{-1}(0)+\psi_{0}(0)}{\alpha-\beta}=\frac{\alpha+\beta}{\sqrt{\alpha^{2}-1}+\sqrt{\beta^{2}-1}}-1\;,$
(91)

which, as expected, is symmetric in $\alpha$ and $\beta$.

The representation of the marginal PMF arising from (88) is clearly singular
in the limit $\alpha-\beta\to 0^{+}$, noting that, in the present queueing
application, we always have $\alpha>\beta>1$. This limit is encountered as
$r_{\text{lo}}\to 0$ or, equivalently, $\nu\to 1$. It is also singular in the
limit $\alpha\to 1^{+}$, which is encountered when $r_{\text{lo}}\to 0$,
$r_{\text{hi}}\to 1$ simultaneously. The representation may be re-arranged and
scaled to make it numerically robust throughout the entire ergodic region. The
final algorithm is detailed in the implementation summary at the end of this
section.

###  6.2 Joint Distribution

In order to facilitate an exact representation of the joint marginal PMF, we
introduce a family of integrals to which we shall refer as the joint
integrals. These are defined by

$\displaystyle\begin{aligned}
J_{m}^{\ell}(\alpha,\beta)&\equiv\int_{0}^{\pi}d\theta\,\frac{1-\cos^{2}\theta}{(\alpha+\cos\theta)^{m+1}(\beta+\cos\theta)}U_{\ell}(\cos\theta)\\\
&=\int_{0}^{\pi}d\theta\,\frac{\sin\theta{\cdot}\sin((\ell+1)\theta)}{(\alpha+\cos\theta)^{m+1}(\beta+\cos\theta)}\;.\end{aligned}$
(92)

A relationship with the marginal integrals is given by

$\displaystyle\begin{aligned}
J_{m}^{0}(\alpha,\beta)&=M_{m}(\alpha,\beta)\;,\\\
J_{m}^{1}(\alpha,\beta)&=2M_{m-1}(\alpha,\beta)-2\alpha
M_{m}(\alpha,\beta)\;.\end{aligned}$ (93)

In the general case, the recurrence relation for the Chebyshev polynomial
$U_{\ell}(\cos\theta)$ implies the recurrence relation for the joint integral

$J_{m}^{\ell+1}(\alpha,\beta)=2J_{m-1}^{\ell}(\alpha,\beta)-2\alpha
J_{m}^{\ell}(\alpha,\beta)-J_{m}^{\ell-1}(\alpha,\beta)\;,$ (94)

for $\ell=0,1,\ldots$, where we take the seed values

$J_{m}^{-1}(\alpha,\beta)\equiv 0\;,\quad
J_{m}^{0}(\alpha,\beta)=M_{m}(\alpha,\beta)\;,$ (95)

for $m=-1,0,1,\ldots$, noting that
$M_{-1}(\alpha,\beta)=\pi(\beta-\sqrt{\beta^{2}-1})$.

Next, we introduce the intermediate integrals

$\displaystyle\begin{aligned}
I_{m}^{\ell}(\alpha,\beta)&\equiv\int_{0}^{\pi}d\theta\,\frac{\cos(\ell\theta)}{(\alpha+\cos\theta)^{m+1}(\beta+\cos\theta)}\;,\\\
I_{m}^{\ell}(\alpha)&\equiv\int_{0}^{\pi}d\theta\,\frac{\cos(\ell\theta)}{(\alpha+\cos\theta)^{m+1}}\;,\end{aligned}$
(96)

for $m=-1,0,\ldots$, $\ell\in\mathbb{Z}$. One may observe that
$M_{n}(\alpha,\beta)=\left[I_{n}^{0}(\alpha,\beta)-I_{n}^{2}(\alpha,\beta)\right]/2$
and that $I_{-1}^{\ell}(\alpha,\beta)=I_{0}^{\ell}(\beta)$,
$I_{-1}^{\ell}(\alpha)=\pi\delta_{\ell 0}$. A partial fraction decomposition
yields

$I_{m}^{\ell}(\alpha,\beta)=c_{m}I_{0}^{\ell}(\beta)-\sum_{k=0}^{m}c_{k}I_{m-k}^{\ell}(\alpha)\;,$
(97)

where we have set $c_{m}\equiv 1/(\alpha-\beta)^{m+1}$. The joint integral may
be expressed in terms of the intermediate integrals as

$J_{m}^{\ell}(\alpha,\beta)=I_{m}^{\ell}(\alpha,\beta)+\beta
I_{m}^{\ell+1}(\alpha,\beta)-I_{m}^{\ell+1}(\alpha)\;.$ (98)

In terms of the joint integrals, we can write

$P_{\text{cut}}(\ell,m)=\frac{(1-r)r_{\text{lo}}}{2\pi
r}\left[\frac{(-\sqrt{r_{\text{hi}}})^{\ell}}{(x_{\text{dif}}/2)^{m-1}}J_{m-1}^{\ell}(\alpha,\beta)-\frac{(-\sqrt{r_{\text{hi}}})^{\ell-1}}{(x_{\text{dif}}/2)^{m}}J_{m}^{\ell-1}(\alpha,\beta)\right]\;.$
(99)

###  6.3 $Q$-Functions

To proceed further, we define a family of functions $Q_{m}^{\ell}(\chi)$, that
generalize the Legendre polynomials and are related closely to the associated
Legendre functions, by means of the integral representation

$Q_{m}^{\ell}(\chi)\equiv\frac{(-1)^{\ell}}{\pi}\left(\frac{\chi+1}{\chi-1}\right)^{\ell/2}\int_{0}^{\pi}d\theta\,\frac{\cos(\ell\theta)}{(\chi+\sqrt{\chi^{2}-1}\cos\theta)^{m+1}}\;.$
(100)

This may be compared with the integral representation of the associated
Legendre function given by [13, section 8.711.2]

$P_{m}^{\ell}(\chi)=(-1)^{\ell}\frac{\Gamma(m+1)}{\Gamma(m-\ell+1)}\int_{0}^{\pi}\frac{d\theta}{\pi}\,\frac{\cos(\ell\theta)}{(\chi+\sqrt{\chi^{2}-1}\cos\theta)^{m+1}}\;,$
(101)

noting also that the Legendre polynomials are recovered according to
$P_{m}(\chi)=P_{m}^{0}(\chi)$. Clearly $Q_{-1}^{\ell}(\chi)=\delta_{\ell 0}$,
$Q_{0}^{\ell}(\chi)=1$ and $Q_{m}^{0}(\chi)=P_{m}(\chi)$. It also follows that

$Q_{m}^{1}(\chi)=(\alpha+\sqrt{\alpha^{2}-1})\left[\alpha
P_{m}(\chi)-\sqrt{\alpha^{2}-1}P_{m-1}(\chi)\right]\;,$ (102)

for $m=0,1,\ldots$. From their definition and the trigonometric identity

$\cos[(\ell+1)\theta]+\cos[(\ell-1)\theta]=2\cos(\ell\theta)\cos\theta\;,$
(103)

it is easily shown that the $Q$-functions satisfy the recurrence relation

$Q_{m}^{\ell+1}(\chi)=\frac{2\chi}{\chi-1}Q_{m}^{\ell}(\chi)-\frac{2}{\chi-1}Q_{m-1}^{\ell}(\chi)-\frac{\chi+1}{\chi-1}Q_{m}^{\ell-1}(\chi)\;,$
(104)

for $m=0,1,\ldots$, $\ell\in\mathbb{Z}$. An explicit expansion of the
$Q$-function is given by

$Q_{m}^{\ell}(\chi)=\sum_{k=0}^{m}\binom{m+\ell}{k+\ell}\binom{k+m}{m}\left(\frac{\chi-1}{2}\right)^{k}\;.$
(105)

The relationship with associated Legendre functions is established by noting
that, for $\ell\leq m$,

$Q_{m}^{\ell}(\chi)=\frac{(m-\ell)!}{m!}\left(\frac{\chi+1}{\chi-1}\right)^{\ell/2}P_{m}^{\ell}(\chi)\;.$
(106)

Another recurrence relation, which follows directly from the recurrence
relations for the associated Legendre functions, is given by

$Q_{m+1}^{\ell}(\chi)=\frac{2m+1}{m+1}\chi
Q_{m}^{\ell}(\chi)-\frac{m^{2}-\ell^{2}}{m(m+1)}Q_{m-1}^{\ell}(\chi)\;,$ (107)

for $m=1,2,\ldots$, with the first few functions being

$Q_{0}^{\ell}(\chi)=1\;,\quad Q_{1}^{\ell}(\chi)=\chi+\ell\;,\quad
Q_{2}^{\ell}(\chi)=\tfrac{3}{2}\chi^{2}+\tfrac{3}{2}\ell\chi+\tfrac{1}{2}(\ell^{2}-1)\;.$
(108)

This recurrence is particularly useful for practical purposes since, while
superficially appearing to be a two-dimensional recurrence, it is actually a
family of independent one-dimensional recurrences in $m$, each indexed by the
value of $\ell$.

The generating function for the $Q$-functions is found to be

$\sum_{m=0}^{\infty}z^{m}Q_{m}^{\ell}(\chi)=\frac{1}{\sqrt{z^{2}-2\chi
z+1}}\left(\frac{\chi-z-\sqrt{z^{2}-2\chi z+1}}{\chi-1}\right)^{\ell}\;.$
(109)

Using

$\chi=\frac{\alpha}{\sqrt{\alpha^{2}-1}}\;,\quad
z=\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\;,$ (110)

and noting that $z^{2}-2\chi z+1=(\beta^{2}-1)/(\alpha^{2}-1)$, this may be
cast in the alternative form

$\sum_{m=0}^{\infty}z^{m}Q_{m}^{\ell}(\chi)=\sqrt{\frac{\alpha^{2}-1}{\beta^{2}-1}}{\cdot}\frac{d_{\ell}(\beta)}{d_{\ell}(\alpha)}\;,$
(111)

where we have set

$d_{\ell}(\alpha)\equiv(\alpha-\sqrt{\alpha^{2}-1})^{\ell}=\frac{1}{(\alpha+\sqrt{\alpha^{2}-1})^{\ell}}\;,$
(112)

and analogously for $d_{\ell}(\beta)$.

The single-argument intermediate integral can be expressed in terms of the
$Q$-function as

$\displaystyle\begin{aligned}
I_{m}^{\ell}(\alpha)&=\frac{\pi(-1)^{\ell}}{(\alpha+\sqrt{\alpha^{2}-1})^{\ell}(\alpha^{2}-1)^{(m+1)/2}}Q_{m}^{\ell}\left(\frac{\alpha}{\sqrt{\alpha^{2}-1}}\right)\\\
&=\pi(-1)^{\ell}\left(\frac{\chi-1}{\chi+1}\right)^{\ell/2}{\cdot}(\chi^{2}-1)^{(m+1)/2}Q_{m}^{\ell}(\chi)\;.\end{aligned}$
(113)

Consequently, we may express the joint integral in terms of the $Q$-functions,
according to

$\begin{split}\frac{(-1)^{\ell}J^{\ell}_{m}(\alpha,\beta)}{\pi
c_{m}}&=d_{\ell+1}(\beta)+d_{\ell+1}(\alpha)\left(\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\right)^{m+1}Q_{m}^{\ell+1}(\chi)\\\
&\quad{}-\frac{1}{\sqrt{\alpha^{2}-1}}\sum_{k=0}^{m}\left(\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\right)^{k}\left[d_{\ell}(\alpha)Q_{k}^{\ell}(\chi)-\beta
d_{\ell+1}(\alpha)Q_{k}^{\ell+1}(\chi)\right]\;.\end{split}$ (114)

###  6.4 Implementation

We conclude this section by giving complete, self-contained specifications of
the algorithms for obtaining exact closed-form algebraic expressions for the
marginal and joint queue-length PMFs that are also numerically efficient and
robust. Let us first recall that

$\alpha=\frac{r+1}{2\sqrt{r_{\text{hi}}}}\;,\quad\beta=\frac{1}{2}\left(\frac{r}{\sqrt{r_{\text{hi}}}}+\frac{\sqrt{r_{\text{hi}}}}{r}\right)\;,$
(115)

and that

$\chi=\frac{\alpha}{\sqrt{\alpha^{2}-1}}\;,\quad
z=\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\;,$ (116)

so that

$\chi=\frac{1+r}{\sqrt{(1+r)^{2}-4r_{\text{hi}}}}\;,\quad
z=\frac{1-r_{\text{hi}}/r}{\sqrt{(1+r)^{2}-4r_{\text{hi}}}}\;,\quad
t\equiv\chi z\;.$ (117)

It is useful to note that

$d_{1}(\alpha)=\frac{1}{2\sqrt{r_{\text{hi}}}}\left[1+r-\sqrt{(1+r)^{2}-4r_{\text{hi}}}\right]=\zeta_{-}(0)/\sqrt{r_{\text{hi}}}\;,$
(118)

and that

$d_{1}(\beta)=\min\left(r/\sqrt{r_{\text{hi}}},\sqrt{r_{\text{hi}}}/r\right)\;.$
(119)

It is also useful to note the bounds on $z$ given by

$0\leq\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\leq\frac{1}{1+r}<1\;,$ (120)

which shows that, while $z$ never becomes large, it can be very small. The
lower and upper endpoints for $z$ correspond to $r_{\text{lo}}=0$ and
$r_{\text{hi}}=0$, respectively. We also have

$\displaystyle\begin{aligned} 1&\leq\chi\leq(1+r)/(1-r)\;,\\\ 0&\leq\chi z\leq
1/(1+r)\;.\end{aligned}$ (121)

Table 6.4 enumerates the singular limits of the marginal and joint integrals.
The algorithms described in the following two subsections implement scalings
and algebraic manipulations to mitigate their effects, and render the final
expression for the marginal and joint PMF numerically robust across the entire
ergodic region.

Singular Regions  $\alpha,\beta$ $r_{\text{lo}},r_{\text{hi}}$ $z$ $\beta\to
1^{+}$ $r_{\text{hi}}\to r^{2}$ $\alpha-\beta\to 0^{+}$ $r_{\text{lo}}\to 0$
$z\to 0$ $\alpha-\beta\to\infty$ $r_{\text{hi}}\to 0$ $z\to 1/(1+r)$
$\alpha\to 1^{+}$ $r_{\text{lo}}\to 0,r_{\text{hi}}\to 1$ $z\to 0$

####  6.4.1 Marginal Distribution

We introduce scaled Legendre polynomials $\hat{P}_{n}(t;z)\equiv
z^{n}P_{n}(t/z)$ that are computed according to the recurrence relation

$\hat{P}_{n}(t;z)=(2-1/n)t\hat{P}_{n-1}(t;z)-(1-1/n)z^{2}\hat{P}_{n-2}(t;z)\;,$
(122)

for $n=1,2,\ldots$, subject to the boundary values

$\hat{P}_{-1}(t;z)=1/z\;,\quad\hat{P}_{0}(t;z)=1\;.$ (123)

The scaled marginal integral

$M_{\text{scl}}(n)\equiv\pi^{-1}(\alpha-\beta)^{n}M_{n}(\alpha,\beta)\;,$
(124)

in terms of which the cut contribution to the low-priority marginal PMF is
expressed as

$P_{\text{cut}}(n)=\frac{(1-r)r^{n-1}}{2}{\cdot}M_{\text{scl}}(n)\;,$ (125)

is computed from the scaled Legendre polynomials as follows. With

$D\equiv z^{2}-2\chi z+1=(\beta^{2}-1)/(\alpha^{2}-1)\;,$ (126)

we can write

$M_{\text{scl}}(n)=-z\hat{P}_{n-1}(t;z)+(2\chi-z)\hat{P}_{n}(t;z)-\frac{\sqrt{D}}{z}\biggl{[}1-\sqrt{D}\sum_{k=0}^{n}\hat{P}_{k}(t;z)\biggr{]}\;.$
(127)

It is instructive to note the asymptotic behaviour

$1-\sqrt{D}\sum_{k=0}^{n}\hat{P}_{k}(t;z)\mathrel{\mathop{\sim}_{n\to\infty}}\sqrt{\frac{\beta+1}{\beta-1}}z\hat{P}_{n}(t;z)\;,$
(128)

provided $\beta\neq 1$. When $\beta=1$, the foregoing expression is
identically unity for all $n<\infty$ since $\sqrt{D}=0$. This enables us to
establish the bounding property

$0<\frac{\sqrt{D}}{z}\biggl{[}1-\sqrt{D}\sum_{k=0}^{n}\hat{P}_{k}(t;z)\biggr{]}<\sqrt{\frac{\beta+1}{\beta-1}}{\cdot}\sqrt{D}\hat{P}_{n}(t;z)\;.$
(129)

The generating function for the Legendre polynomials implies that

$\sqrt{D}\sum_{k=0}^{\infty}\hat{P}_{k}(t;z)=1\;.$ (130)

Use of this identity leads directly to an alternative expression for the
scaled marginal integral, given by

$M_{\text{scl}}(n)=-z\hat{P}_{n-1}(t;z)+(2\chi-z)\hat{P}_{n}(t;z)-\frac{D}{z}\sum_{k=n+1}^{\infty}\hat{P}_{k}(t;z)\;.$
(131)

In the numerical implementation, we construct the last term of (127) and test
whether it satisfies the bound given by (129). If the test is passed, then the
scaled marginal integral is evaluated via (127). If it fails (due to
arithmetic underflow), then the scaled marginal integral is evaluated instead
via (131), where the sum is computed by sequentially adding terms until the
desired level of convergence is attained.

####  6.4.2 Joint Distribution

To deal with the joint distribution, we introduce scaled $Q$-functions
according to $\hat{Q}^{\ell}_{m}(t;z)\equiv z^{m}Q^{\ell}_{m}(t/z)$. Let

$A_{m}\equiv\frac{2m+1}{m+1}{\cdot}t\;,\quad
B_{m}^{\ell}\equiv\frac{m^{2}-\ell^{2}}{m(m+1)}{\cdot}z^{2}\;.$ (132)

Then,

$\hat{Q}_{m+1}^{\ell}(t,z)=A_{m}\hat{Q}_{m}^{\ell}(t,z)-B_{m}^{\ell}\hat{Q}_{m-1}^{\ell}(t,z)\;,$
(133)

for $\ell=0,1,\ldots$, $m=1,2,\ldots$, subject to the boundary values

$\hat{Q}_{0}^{\ell}(t,z)=1\;,\quad\hat{Q}_{1}^{\ell}(t,z)=t+\ell z\;.$ (134)

We use the scaled $Q$-functions to construct the scaled joint integrals

$J_{\text{scl}}(\ell,m)\equiv\pi^{-1}(-1)^{\ell}(\alpha-\beta)^{m+1}J_{m}^{\ell}(\alpha,\beta)\;,$
(135)

as

$\begin{split}J_{\text{scl}}(\ell,m)&=d_{\ell+1}(\beta)+d_{\ell+1}(\alpha)\biggl{[}z\hat{Q}_{m}^{\ell+1}(\chi
z,z)\\\ &\quad{}-(1+\chi)\biggl{(}\sum_{k=0}^{m}\hat{Q}_{k}^{\ell}(\chi
z,z)-\frac{1-z/\chi}{1+1/\chi}\sum_{k=0}^{m}\hat{Q}_{k}^{\ell+1}(\chi
z,z)\biggr{)}\biggr{]}\;,\end{split}$ (136)

and we may note the generalization of (130), given by

$\sqrt{D}\sum_{k=0}^{\infty}\hat{Q}_{k}^{\ell}(t;z)=\frac{d_{\ell}(\beta)}{d_{\ell}(\alpha)}\;.$
(137)

One can also confirm that
$J_{\text{scl}}(\ell,0)=d_{\ell+1}(\beta)-d_{\ell+1}(\alpha)$.

Given the definition of the scaled joint integrals in (135), the cut
contribution to the wait-conditional joint PMF becomes

$P_{\text{cut}}(\ell,m)=(1-r)\sqrt{r_{\text{hi}}}\left[r_{\text{hi}}^{\ell/2}r^{m-1}J_{\text{scl}}(\ell,m-1)-r_{\text{hi}}^{(\ell-1)/2}r^{m}J_{\text{scl}}(\ell-1,m)\right]\;.$
(138)

As a sanity check, one can confirm, with the aid of (118) and (119), that this
expression leads to the results

$\displaystyle\begin{aligned}
P_{\text{cut}}(0,0)&=(1-r){\cdot}\min\left(1,r_{\text{hi}}/r^{2}\right)\,,\\\
P_{\text{cut}}(0,m+1)&=r_{\text{lo}}{\cdot}P_{\text{cut}}(m)\;,\\\
P_{\text{cut}}(\ell,0)&=P_{\text{xhi}}(\ell)-P_{\text{pol}}(\ell,0)\;,\end{aligned}$
(139)

where the cut contribution to the low-priority marginal PMF
$P_{\text{cut}}(m)$ is given by (125), $P_{\text{pol}}(\ell,0)$ follows from
(84) and $P_{\text{xhi}}(\ell)\equiv P(\ell,0)=(1-r)\zeta_{-}^{\ell}(0)$ [32].
An intermediate step in deriving the last relationship is

$P_{\text{cut}}(\ell,0)=(1-r)r_{\text{hi}}^{\ell/2}\left[d_{\ell}(\alpha)-d_{\ell}(\beta)+(\sqrt{r_{\text{hi}}}/r){\cdot}d_{\ell+1}(\beta)\right]\;,$
(140)

which is established with the aid of the identity

$I_{0}^{\ell}(\gamma)+\gamma
I_{0}^{\ell+1}(\gamma)=\pi(-1)^{\ell}d_{\ell+1}(\gamma)\;,$ (141)

with $\gamma=\alpha,\beta$. In fact, it may be noted that

$P_{\text{pol}}(\ell,m)=(1-r)r^{m-1}r_{\text{hi}}^{\ell/2}\left[rd_{\ell}(\beta)-\sqrt{r_{\text{hi}}}d_{\ell+1}(\beta)\right]\;.$
(142)

Finally, we mention that the large-$\ell$ asymptotics of the scaled joint
integrals are given by

$J_{\text{scl}}(\ell,m)\mathrel{\mathop{\sim}_{\ell\to\infty}}d_{\ell+1}(\beta)-\frac{(\ell+1)^{m}}{m!}\left(\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\right)^{m}d_{\ell+1}(\alpha)\;.$
(143)

A derivation is presented in Appendix 11. The large-$\ell$ behaviour of the
joint PMF $P(\ell,m)$ for fixed $m$ follows directly from this.

##  7 Numerical Results

In the present work, we have introduced two approaches to the practical
calculation of both the marginal and joint queue-length distributions. We
shall test the numerical performance of these two methods against each other
and against two previously developed methods, namely, the quadratic-recurrence
method from [31] and the Fast Fourier Transform (FFT) method from [32]. We
subject the algorithms to a number of tests, each with an associated measure
of performance (MOP) as follows:

x Pairwise comparisons among four algorithms for the low-priority marginal
distribution of the queue length as a function of the fraction of high-
priority arrivals $\nu$, with total traffic intensity $r=0.95$. The values on
the vertical axis indicate the number of decimal places of agreement.

x Aggregation tests for the joint probability distribution of the queue
lengths as a function of the fraction of high-priority arrivals $\nu$, with
total traffic intensity $r=0.9$, comparing four algorithms as indicated. The
values on the vertical axis indicate the number of decimal places of
agreement.

x Xhi-tests for the joint probability distribution of the queue lengths as a
function of the fraction of high-priority arrivals $\nu$, with total traffic
intensity $r=0.9$, comparing four algorithms as indicated. The values on the
vertical axis indicate the number of decimal places of agreement.

x Xlo-tests for the joint probability distribution of the queue lengths as a
function of the fraction of high-priority arrivals $\nu$, with total traffic
intensity $r=0.9$, comparing four algorithms as indicated. The values on the
vertical axis indicate the number of decimal places of agreement.

##### Pairwise Test:

We compare the low-priority marginal PMFs computed via the various methods
with each other in a pairwise manner. To compare method A with method B, we
adopt the MOP

$\Xi_{\text{pair}}\equiv-\max_{m\geq
0}\left\\{\log_{10}\left(|\ln(P_{\text{lo}}^{(\text{A})}(m))-\ln{(}P_{\text{lo}}^{(\text{B})}(m))|\right)\right\\}\;,$
(144)

where the maximum is taken over all values $0\leq m\leq n_{\text{lim}}$ such
that

$\min\\{P_{\text{lo}}^{(\text{A})}(m),P_{\text{lo}}^{(\text{B})}(m)\\}>p_{\text{lim}}>0\;.$
(145)

Since we are working in double-precision arithmetic222All computation is
performed in Matlab R2020a, which implements IEEE Standard 754 for double
precision., all MOPs of this kind are capped at a maximum allowed value of
$16$. The interpretation of $\Xi_{\text{pair}}$ (and similarly for all of the
subsequent MOPs) is that it indicates the number of decimal places of
numerical agreement in the worst case. Finite values for the limiting
probability level $p_{\text{lim}}$, typically around $10^{-12}$, are adopted
as one cannot expect any digital algorithm to work infinitely far into the
tail. Finite values of the limiting queue occupancy $n_{\text{lim}}$ are adopt
to keep computational load under control.

##### Aggregation Test:

The aggregated queue-length distribution describes the total number of
entities in the queue, regardless of priority level. The exact result for the
aggregated PMF is [31, 32]

$P_{\text{agg}}^{\text{(ex)}}(k)=(1-r)r^{k}\;,$ (146)

for $k=0,1,\ldots$. We construct the aggregated probabilities from the
numerically computed joint PMF according to

$P_{\text{agg}}(k)=\sum_{m=0}^{k}P(k-m,m)$ (147)

and test against the exact result using the MOP

$\Xi_{\text{agg}}\equiv-\max_{k\geq
0}\left\\{\log_{10}\left(|\ln(P_{\text{agg}}(k))-\ln{(}P_{\text{agg}}^{(\text{ex})}(k))|\right)\right\\}\;,$
(148)

where the maximum is taken over all values $0\leq k\leq n_{\text{lim}}$ such
that $P_{\text{agg}}^{(\text{ex})}(k)>p_{\text{lim}}>0$.

##### Xhi-Test:

The exact result for the exclusively-high queue-length probabilities is [32]

$P_{\text{xhi}}^{\text{(ex)}}(\ell)\equiv
P(\ell,0)=(1-r){\cdot}\left[\frac{1+r-\sqrt{(1+r)^{2}-4r_{\text{hi}}}}{2}\right]^{\ell}\;,$
(149)

for $\ell=0,1,\ldots$. We construct the exclusively-high probabilities from
the numerically computed joint PMF according to
$P_{\text{xhi}}(\ell)=P(\ell,0)$ and test against the exact result using the
MOP

$\Xi_{\text{xhi}}\equiv-\max_{\ell\geq
0}\left\\{\log_{10}\left(|\ln(P_{\text{xhi}}(\ell))-\ln{(}P_{\text{xhi}}^{(\text{ex})}(\ell))|\right)\right\\}\;,$
(150)

where the maximum is taken over all values $0\leq\ell\leq n_{\text{lim}}$ such
that $P_{\text{xhi}}^{(\text{ex})}(\ell)>p_{\text{lim}}>0$.

##### Xlo-Test:

An exact relationship between that exclusively-low queue probabilities and the
low-priority marginal PMF is given by [31, 32]

$P_{\text{xlo}}^{(\text{ex})}(m)\equiv
P(0,m)=(1-r)\delta_{m0}+(1-\delta_{m0})r_{\text{lo}}{\cdot}P_{\text{lo}}(m-1)\;,$
(151)

for $m=0,1,2,\ldots$, where we can formally set $P_{\text{lo}}(-1)\equiv 0$.
We construct the exclusively-low probabilities from the numerically computed
joint PMF according to $P_{\text{xlo}}(m)=P(0,m)$ and test against the low-
priority marginal using the MOP

$\Xi_{\text{xlo}}\equiv-\max_{m>0}\left\\{\log_{10}\left(|\ln(P_{\text{xlo}}(m))-\ln{(}r_{\text{lo}}P_{\text{lo}}(m-1))|\right)\right\\}\;,$
(152)

where the maximum is taken over all values $0<m\leq n_{\text{lim}}$ such that
$P_{\text{xlo}}^{(\text{ex})}(m)>p_{\text{lim}}>0$.

In Figures 7–7, the MOP values relevant to each test are displayed on the
vertical axis against the full range of high-priority arrival fraction
(hifrac) $0\leq\nu\leq 1$ on the horizontal axis. For the marginal
distribution, Figure 7 presents results of the pairwise tests for $r=0.95$
with $p_{\text{lim}}=10^{-12}$. The blue curve compares the Legendre
polynomial method of Section 6.4.1 and the quadratic recurrence method from
[31]. The green curve compares the Gaussian quadrature method obtained from
(21) and (33) of Section 3 and the aforementioned recurrence method. The red
curve compares the quadrature and Legendre methods. The cyan curve compares
the recurrence method and the FFT method from [32].

Here, and subsequently, $n_{\text{max}}$ denotes the maximum queue-length
actually sampled in achieving the limiting probability level $p_{\text{lim}}$.
The FFT method is configured to a relative error of $10^{-12}$, at least for
small values of $\nu$, and this is reflected in all of the graphs. The fact
that the blue and red curves dip in accuracy at the critical point ($\nu=r$),
but the green curve does not, indicates that only the Legendre method
struggles in this region, while the quadratic-recurrence and quadrature
methods are able to maintain accuracy. The fact that the blue and green curves
(where the quadratic recurrence is the common method) align closely but lie
below the red curve suggests that the quadrature and Legendre methods are
slightly more accurate, at least away from the critical point (indicated by
the dashed black line).

In Figures 7–7, that test the joint distribution, the red curves correspond to
the Gaussian quadrature method of Section 5. The blue curves correspond to the
Legendre polynomial method of Section 6.4.2. The green curves result from the
quadratic recurrence method from [31], and the cyan curves from the FFT method
described in [32]. Figure 7 presents results of the aggregation test for
$r=0.9$ with $p_{\text{lim}}=10^{-9}$. Figure 7 presents results of the xhi-
test for $r=0.9$ with $p_{\text{lim}}=10^{-12}$. Figure 7 presents results of
the xlo-test for $r=0.9$ with $p_{\text{lim}}=10^{-12}$. All methods exhibit
good levels of numerical performance.

x Two-dimensional level-set mapping for the joint probability distribution of
the queue lengths for total traffic intensity $r=0.9$ and fraction of high-
priority arrivals at the critical point $\nu=r$.

x Error histogram, comparing the quadrature and quadratic-recurrence methods,
for the joint probability distribution of the queue lengths for total traffic
intensity $r=0.9$ and fraction of high-priority arrivals at the critical point
$\nu=r$.

x CDF sections for the joint probability distribution of the queue lengths for
total traffic intensity $r=0.9$ and fraction of high-priority arrivals at the
critical point $\nu=r$.

In order to illustrate the qualitative features of the queue-length
distribution, in Figure 7, we present a two-dimensional contour plot with
level curves for the joint queue-length PMF computed via the quadrature method
at the critical point ($\nu=r$) for $r=0.9$. This is computationally the
worst-case scenario for a given total traffic intensity $r$. The level curves
are displayed on a (base-$10$) logarithmic scale down to a limiting
probability level of $p_{\text{lim}}=10^{-12}$, with the low-priority
occupancy on the horizontal axis and the high-priority occupancy on the
vertical axis. We also benchmark this case against the quadratic-recurrence
method, with the results for the number of decimal places of agreement
presented as a histogram in Figure 7. Finally, in Figure 7, we plot one-
dimensional sections of the cumulative distribution function (CDF) for various
constant values of the high-priority queue occupancy, i.e.

$F_{\ell}(n)\equiv\sum_{m=0}^{n}P(\ell,m)\;.$ (153)

The curves correspond to high-priority queue-lengths $\ell=0,1,\ldots\,14$,
with the uppermost black curve representing $\ell=100$.

##  8 Conclusions

The analytic approach to the calculation of the joint and marginal queue-
length distributions for the two-level non-preemptive priority queue studied
here leads to integral representations that are amenable to both efficient
numerical quadrature rules and exact expressions as finite sums of easily
computable special functions that generalize the Legendre polynomials. One
benefit of integral forms obtained here is that the large queue-length
asymptotic behaviour becomes apparent. Another is that there is scope for
extending the methodology to cater for additional complications, such as
multiple arrival classes.

##  9 No-Wait Probability

Let $c$ denote the number of servers, and let $\mathscr{N}_{\text{sys}}$ be
the RV for system occupancy. It is well-known that
$\Pr(\mathscr{N}_{\text{sys}}=n)=Ar^{n}$ for all $n\geq c$, for some constant
$A$ and, by definition, we have
$P_{\text{NW}}=\Pr(\mathscr{N}_{\text{sys}}\leq c-1)$ for the no-wait
probability. Therefore,

$\displaystyle\begin{aligned}
P_{\text{NW}}&=\sum_{n=0}^{c-1}\Pr(\mathscr{N}_{\text{sys}}=n)\\\
&=1-\sum_{n=c}^{\infty}\Pr(\mathscr{N}_{\text{sys}}=n)\\\ &=1-Ar^{c}/(1-r)\\\
&=1-\Pr(\mathscr{N}_{\text{sys}}=c)/(1-r)\;.\end{aligned}$ (154)

Furthermore,
$\Pr(\mathscr{N}_{\text{sys}}=c)=[(rc)^{c}/c!]\Pr(\mathscr{N}_{\text{sys}}=0)$,
and [12]

$\frac{1}{\Pr(\mathscr{N}_{\text{sys}}=0)}=\sum_{n=0}^{c}\frac{(rc)^{n}}{n!}+\frac{(rc)^{c}}{c!(1/r-1)}\;.$
(155)

The result for $P_{\text{NW}}$ follows, and may be compactly stated as

$\frac{1}{1-P_{\text{NW}}}=1+\frac{(1-r)c!}{(rc)^{c}}\sum_{n=0}^{c-1}\frac{(rc)^{n}}{n!}\;.$
(156)

##  10 Laguerre Series for the Waiting Time

Let us write $p_{n}\equiv P_{\mathscr{L}}(n)$ so that

$p_{n}=\frac{1}{n!}{\cdot}\left.\frac{d^{n}}{dp^{n}}g_{\text{lo}}(z)\right|_{z=0}=\oint_{\mathcal{C}_{\eta}}\frac{dz}{2\pi
i}\,\frac{g_{\text{lo}}(z)}{z^{n+1}}\;,$ (157)

where the closed anti-clockwise traversed integration contour
$\mathcal{C}_{\eta}$, initially taken to be a circle centred on the origin of
radius $0<\eta<1/r$, is deformed into the Bromwich contour $\Gamma$ on which
$z=\delta+iy$ with $-\infty<y<+\infty$, for some $1<\delta<1/r$ so that
$\mathop{\rm Re}\nolimits z>1$. Then, we consider the Binomial transform of
the $p_{n}$ sequence

$\displaystyle\begin{aligned}
\gamma_{n}&\equiv\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}p_{k}\\\
&=\int_{\Gamma}\frac{dz}{2\pi
i}\,g_{\text{lo}}(z)\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{z^{k+1}}\\\
&=\int_{\Gamma}\frac{dz}{2\pi
i}\,\frac{g_{\text{lo}}(z)}{z}{\cdot}(1-1/z)^{n}\;.\end{aligned}$ (158)

We now recall the generating function for the Laguerre polynomials $L_{n}(x)$,
given by

$\sum_{n=0}^{\infty}\tau^{n}L_{n}(x)=\frac{1}{1-\tau}e^{-\tau x/(1-\tau)}\;.$
(159)

On setting $\tau=1-1/z$ and $x=r_{\text{lo}}t$, this reads

$\sum_{n=0}^{\infty}\left(1-1/z\right)^{n}L_{n}(r_{\text{lo}}t)=ze^{-r_{\text{lo}}t(z-1)}\;,$
(160)

from which it follows that

$\sum_{n=0}^{\infty}\gamma_{n}L_{n}(r_{\text{lo}}t)=\int_{\Gamma}\frac{dz}{2\pi
i}\,g_{\text{lo}}(z)e^{-r_{\text{lo}}t(z-1)}\;.$ (161)

Now, (47) can be written as

$P_{\mathscr{L}}(k)=\int_{0}^{\infty}dt^{\prime}\,P^{\prime}_{\mathscr{W}}(t^{\prime})\frac{t^{\prime
k}}{k!}e^{-t^{\prime}}\;,$ (162)

where $t^{\prime}=r_{\text{lo}}t$ and
$P^{\prime}_{\mathscr{W}}(t^{\prime})dt^{\prime}=P_{\mathscr{W}}(t)dt$. We
also have

$L_{n}(t)=\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{k}}{k!}t^{k}\;.$ (163)

These may be combined to yield

$\int_{0}^{\infty}dt^{\prime}\,e^{-t^{\prime}}L_{n}(t^{\prime})P^{\prime}_{\mathscr{W}}(t^{\prime})=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\int_{0}^{\infty}dt^{\prime}\,P^{\prime}_{\mathscr{W}}(t^{\prime})\frac{t^{\prime
k}}{k!}e^{-t^{\prime}}=\gamma_{n}\;.$ (164)

The completeness relation [20]

$\sum_{n=0}^{\infty}L_{n}(t)L_{n}(t^{\prime})=e^{(t+t^{\prime})/2}{\cdot}\delta(t-t^{\prime})\;,$
(165)

$\delta(t)$ being the Dirac delta function, implies the convergent Laguerre-
series representation

$P^{\prime}_{\mathscr{W}}(t^{\prime})=\sum_{n=0}^{\infty}\gamma_{n}L_{n}(t^{\prime})\;,$
(166)

and we conclude, via (161), that

$P_{\mathscr{W}}(t)dt=r_{\text{lo}}dt{\cdot}\int_{\Gamma}\frac{dz}{2\pi
i}\,g_{\text{lo}}(z)e^{-r_{\text{lo}}t(z-1)}\;.$ (167)

The Bromwich contour $\Gamma$, with $\real z>1$, may now be deformed as
$\Gamma\mapsto\mathcal{C}_{\text{pol}}\cup\mathcal{C}_{\text{cut}}$ so that it
wraps around the pole and cut of $g_{\text{lo}}(z)$. This serves to reproduce
(52).

The convergence claim for (166) follows from the theorem that the formal
expansion

$\gamma(t)\equiv\sum_{k=0}^{\infty}\gamma_{k}L_{k}(t)$ (168)

converges in $L_{2}$ if and only if $\|\gamma\|^{2}_{L_{2}}<\infty$, where

$\|\gamma\|^{2}_{L_{2}}\equiv\int_{0}^{\infty}dt\,e^{-t}|\gamma(t)|^{2}=\sum_{k=0}^{\infty}|\gamma_{k}|^{2}\;.$
(169)

In the present application, we have $\gamma(t)=P^{\prime}_{\mathscr{W}}(t)$
for which, trivially,

$\|P^{\prime}_{\mathscr{W}}\|_{L_{2}}^{2}=\int_{0}^{\infty}dt\,e^{-t}\left[P^{\prime}_{\mathscr{W}}(t)\right]^{2}<\int_{0}^{\infty}dt\,e^{-t}=1\;.$
(170)

One should note that, while the coefficients $\gamma_{n}$ of the Laguerre
series (166) are constructed as (potentially) numerically problematic
alternating combinatorial sums, the last equality in (158) shows that they can
be efficiently computed by means of the very same iterative quadrature method
derived for the queue-length probabilities $p_{n}$. Explicitly,

$\begin{split}\gamma_{n}&=\Theta(r^{2}-r_{\text{hi}})\left[1-\frac{r(1-r)}{r_{\text{lo}}}\right]{\cdot}\frac{(1-r)^{n+1}}{r}\\\
&\quad{}+\frac{1-r}{2\pi
r}\int_{0}^{1}du\,\sqrt{u(1-u)}{\cdot}\frac{(u+c)^{n}}{(u+a)^{n+1}(u+b)}\;,\end{split}$
(171)

where $c\equiv(x_{-}-1)/x_{\text{dif}}=(1-r)a+rb$.

##  11 Joint Asymptotics

To determine the large-$\ell$ behaviour of the joint integral (92), we use the
symmetry properties and periodicity of the integrand to write it as an
integral over the interval $\theta\in[0,2\pi)$. We then promote the
integration variable $\theta$ to the complex plane, and deform the contour
from $\theta=0$ to $\theta=2\pi$ into three straight legs: (i) from $0$ to
$0+i\infty$, (ii) from $0+i\infty$ to $2\pi+i\infty$, and (iii) from
$2\pi+i\infty$ to $2\pi$, making small excursions on the vertical legs around
the points $\theta=\acos(\alpha),\acos(\beta)$. The contribution from the
horizontal contour vanishes, while the vertical contours cancel apart from the
excursions around the singular points. As a result, we obtain the
decomposition
$J_{\text{scl}}(\ell,m)=J_{\text{scl}}^{\alpha}(\ell,m)+J_{\text{scl}}^{\beta}(\ell,m)$
where, on writing $\theta=it$, we have

$J_{\text{scl}}^{\alpha,\beta}(\ell,m)=(\beta-\alpha)^{m+1}\oint_{C_{\alpha,\beta}}\frac{dt}{2\pi
i}\,\frac{\sinh(t)e^{-\ell t}}{(\cosh(t)-\alpha)^{m+1}(\cosh(t)-\beta)}\;,$
(172)

where $C_{\alpha}$ is an infinitesimal anti-clockwise circle about the point
$t=\operatorname{acosh}(\alpha)$, and likewise for $C_{\beta}$ about the point
$t=\operatorname{acosh}(\beta)$.

It is clear that $J_{\text{scl}}^{\beta}(\ell,m)=d_{\ell+1}(\beta)$. To
evaluate $J_{\text{scl}}^{\alpha}(\ell,m)$, we set $\alpha\equiv\cosh(s)$ so
that

$\displaystyle\begin{aligned}
J_{\text{scl}}^{\alpha}(\ell{-}1,m)&=\left(\frac{\beta-\alpha}{2}\right)^{m+1}\\!\\!e^{-\ell
s}\oint_{C_{\alpha}}\frac{dt}{2\pi
i}\,\frac{e^{-\ell(t-s)}}{\sinh^{m+1}((t-s)/2)}\\\
&\quad{}\times\frac{\sinh(t)}{\sinh^{m+1}((t+s)/2)[\cosh(t)-\beta]}\\\
&=\left(\frac{\beta-\alpha}{2}\right)^{m+1}\frac{e^{-\ell
s}}{\ell}\oint_{C_{0}}\frac{dt}{2\pi
i}\,\frac{e^{-t}}{\sinh^{m+1}(t/(2\ell))}\\\
&\quad{}\times\frac{\sinh(s+t/\ell)}{\sinh^{m+1}(s+t/(2\ell))[\cosh(s+t/\ell)-\beta]}\;.\end{aligned}$
(173)

Taking the limit yields

$J_{\text{scl}}^{\alpha}(\ell{-}1,m)\mathrel{\mathop{\sim}_{\ell\to\infty}}\frac{(\beta-\alpha)^{m+1}\ell^{m}e^{-\ell
s}}{\sinh^{m}(s)[\cosh(s)-\beta]}\oint_{C_{0}}\frac{dt}{2\pi
i}\,\frac{e^{-t}}{t^{m+1}}=-\frac{\ell^{m}}{m!}\left(\frac{\alpha-\beta}{\sqrt{\alpha^{2}-1}}\right)^{m}d_{\ell}(\alpha)\,,$
(174)

where $\mathcal{C}_{0}$ is an infinitesimal anti-clockwise circle about the
origin. The result in the main text follows. A by-product of this discussion
is an alternative representation of the joint integral, given by

$J_{\text{scl}}(\ell,m)=d_{\ell+1}(\beta)+\oint_{C_{\alpha}}\frac{dw}{2\pi
i}\,\frac{d_{\ell+1}(w)}{w-\beta}\left(\frac{\beta-\alpha}{w-\alpha}\right)^{m+1}\;.$
(175)

## Acknowledgments

The authors gratefully acknowledge useful discussions with Dr. Stephen
Bocquet.

## References

  * Almehdawe et al. [2013] Almehdawe E, Jewkes B, Q-M He (2013) A Markovian queueing model for ambulance offload delays. _European Journal of Operational Research_ 226(3):602–614.
  * Bertsimas and Mourtzinou [1997] Bertsimas D, Mourtzinou G (1997) Multiclass queueing systems in heavy traffic: An asymptotic approach based on distributional and conservation laws. _Operations Research_ 45(3):470–487.
  * Bertsimas and Nakazato [1995] Bertsimas D, Nakazato D (1995) The distributional Little’s law and its applications. _Operations Research_ 43(2):298–310.
  * Chawla [1970] Chawla M (1970) Estimation of errors of Gauss-Chebyshev quadratures. _The Computer Journal_ 13(1):107–109.
  * Cohen [1956] Cohen J (1956) Certain delay problems for a full availability trunk group loaded by two traffic sources. _Philips Telecommunications Review_ 16(3):105–113.
  * Davis [1966] Davis R (1966) Waiting time distribution of a multi-server, priority queuing sytem. _Operations Research_ 14(1):133–136.
  * Erdélyi [1940] Erdélyi A (1940) Some confluent hypergeometric functions of two variables. _Proceedings of the Royal Society of Edinburgh_ 60(3):344–361.
  * Gail et al. [1988] Gail H, Hantler S, Taylor B (1988) Analysis of non-preemptive priority multi-server queue. _Advances in Applied Probability_ 20(4):852–879.
  * Gautschi [2002] Gautschi W (2002) Gauss quadrature approximations to hypergeometric functionsand confluent hypergeometric functions. _Journal of Computational and Applied Mathematics_ 139(1):173–187.
  * Gautschi [2004] Gautschi W (2004) _Orthogonal Polynomials: Computation and Approximation_ (Oxford, UK: Oxford University Press), first edition.
  * Gautschi [2005] Gautschi W (2005) Orthogonal polynomials (in Matlab). _Journal of Computational and Applied Mathematics_ 178(1):425–434.
  * Gnedenko and Kovalenko [1989] Gnedenko B, Kovalenko I (1989) _Introduction to Queueing Theory_ (Boston, MA, USA: Birkhäuser), second edition.
  * Gradshteyn and Ryzhik [2007] Gradshteyn I, Ryzhik I (2007) _Table of Integrals, Series, and Products_ (Burlington, MA, USA: Academic Press), seventh edition.
  * Holley [1954] Holley J (1954) Waiting line subject to priorities. _Operations Research_ 2(3):341–343.
  * Hou and Zhao [2020] Hou J, Zhao X (2020) Using a priority queueing approach to improve emergency department performance. _Journal of Management Analytics_ 7(1):28–43.
  * Kao and Narayanan [1990] Kao E, Narayanan K (1990) Computing steady state probabilities of a non-preemptive priority queue. _INFORMS Journal on Computing_ 2(3):211–218.
  * Karlin and McGregor [1958] Karlin S, McGregor J (1958) Many server queueing processes with Poisson input and exponential service times. _Pacific Journal of Mathematics_ 8(1):87–118.
  * Keilson and Servi [1988] Keilson J, Servi L (1988) A distributional form of Little’s law. _Operations Research Letters_ 7(5):223–227.
  * Kella and Yechiali [1985] Kella O, Yechiali U (1985) Waiting times in the non-preemptive priority M/M/$c$ queue. _Communications in Statistics. Stochastic Models_ 1(2):257–262.
  * Li and Wong [2008] Li Y, Wong R (2008) Integral and series representations of the Dirac delta function. _Communications in Pure and Applied Analysis_ 7(2):229–247.
  * Orman [1995] Orman A (1995) _Models for Scheduling a Multifunction Phased Array Radar System_. Ph.D. thesis, University of Southampton, Southampton, UK.
  * Orman et al. [1996] Orman A, Potts C, Shahani A, Moore A (1996) Scheduling for a multifunction phased array radar system. _European Journal of Operational Research_ 90(1):13–25.
  * Pestalozzi [1964] Pestalozzi G (1964) Priority rules for runway use. _Operations Research_ 12(6):941–950.
  * Picard [1881] Picard E (1881) Sur une extension aux fonctions de deux variables du problème de Riemann relatif aux fonctions hypergéométriques. _Annales Scientifiques de l’École Normale Supérieure, Serie 2_ 10:305–322.
  * Press et al. [2007] Press W, Teukolsky S, Vetterling W, Flannery B (2007) _Numerical Recipes: The Art of Scientific Computing_ (Cambridge, UK: Cambridge University Press), third edition.
  * Rice [1980] Rice S (1980) Distribution of quadratic forms in normal random variables – evaluation by numerical integration. _SIAM Journal on Scientific and Statistical Computing_ 1(4):438–448.
  * Sleptchenko [2003] Sleptchenko A (2003) Multi-class, multi-server, queues with non-preemptive priorities. Report Eurandom Vol. 2003016, EURANDOM, Eindhoven University of Technology, The Netherlands, URL https://pure.tue.nl/ws/portalfiles/portal/3290904/567759.pdf.
  * Taylor and Templeton [1980] Taylor I, Templeton J (1980) Waiting time in a multi-server cutoff-priority queue, and its application to an urban ambulance service. _Operations Research_ 28(5):1168–1188.
  * Trefethen and Weideman [2014] Trefethen L, Weideman J (2014) The exponentially convergent trapezoidal rule. _SIAM Review_ 56(3):385–458.
  * Zuk and Kirszenblat [2023a] Zuk J, Kirszenblat D (2023a) Exact results for the distribution of the partial busy period for a multi-server queue, URL http://dx.doi.org/10.2139/ssrn.4569933, submitted for publication to _Queueing Sytems_.
  * Zuk and Kirszenblat [2023b] Zuk J, Kirszenblat D (2023b) Explicit results for the distributions of queue lengths for a non-preemptive two-level priority queue, URL http://dx.doi.org/10.2139/ssrn.4574550, submitted for publication to _Annals of Operations Research_.
  * Zuk and Kirszenblat [2023c] Zuk J, Kirszenblat D (2023c) Joint queue-length distribution for the non-preemptive multi-server multi-level Markovian priority queue, URL http://dx.doi.org/10.2139/ssrn.4621771, submitted for publication to _Annals of Operations Research_.