MCMC based modelling of queuing systems from empirical data

. Markov chain Monte Carlo (MCMC) modelling technique requires one to be able to construct a proposal density. There is no universal way to achieve this. This paper considers the universal proposal selection technique based on the kernel density estimate. Two channel queuing system with a priority was modelled using this technique. Empirical data (the observed service times) and the rates of arrival processes are all the information used for simulating the system.


Introduction
This paper considers the MCMC approach to modelling service times in the M/G/2/∞ system with a priority. The system is depicted in Fig. 1. The probability distribution of the service times is evaluated from the data observed. The kernel density estimate (KDE) was used for this purpose. It could be difficult to draw samples from a distribution of a complicated or unknown analytical form. Having the empirical data it is possible, to construct a nonparametrical probability density estimate for it. Custom approach of MCMC method helps to sample from such a distribution.
By performing the modelling of the system the mean number of customers L (i) and the mean waiting time W (i) , i ∈ {1; 2}, are obtained. Fig. 2 illustrates a sample run of a two channel queuing system. The results are compared to the theoretical system characteristics according to the formulas (1) and (2) [2]. , ,

Markov Chain Monte Carlo (MCMC)
Suppose a researcher must generate random service times x i , i = 1, n which are distributed according to the π(·). The idea of MCMC is to construct a Markov chain Every Markov chain can be determined through an initial state and a transition kernel (3). It is known that the stationary distribution is unique if Markov chain is ergodic: π(y) = x∈Ω π(x)P (y|x), ∀y ∈ Ω.
Having ergodic and discrete Markov chain with the discrete stationary probabilities π(x i ), the equation (4) holds. Total number of (n − 1) equations and n(n − 1) unknown transition kernel probabilities are apparent.
Thus there exist an infinite number of transition kernels representing the stationary distribution π(x). Any of those transition kernels can be constructed and used for generating x i . One of the most widely used methods for constructing such Markov chain is Metropolis-Hastings algorithm [1]. It is implemented as follows. At first an optional transition kernel Q(.y|x) is chosen. Then there exists a probability α for chosen kernel Q being equal to transition kernel P : By substituting (5) into (4), we have: π(x)Q(y|x)α(y|x) = π(y)Q(x|y)α(x|y), ∀x = y.
By solving (6) and taking in mind the higher acceptance ratio when sampling random numbers [3], it is shown that: Sampling of each x i is performed in 4 steps. Firstly a candidate point x i is drawn from the proposal distribution. Then the probability α i , indicating that this point is also distributed by the target density, is calculated. The next step is to draw From (7) it is evident that π(x) can be determined up to a multiplicative con-

Nonparametric Probability Density Function (PDF) estimation
Consider a sample consisting of random independent and identically distributed values x i . The kernel density estimate was chosen in order to evaluate the probability density of where k(·) is the kernel function, h is its width [4]. The triangular kernel function is useful if the data has sharp edged distribution. Gaussian kernel (9) makes the plot of the estimate's PDF very smooth.
Such probability density estimation could be interpreted as the assignment of kernel density function to each x i plus the weighted sum of all the other assignations. The contribution of any other x j to the probability value at x i decreases if x i − x j increases. The only drawback of such estimation is the necessity of using all the points from the sample while evaluating the probability at a particular point.
Further research requires the cumulative density function of KDE. By integrating the (8) we get (10) where K(·) is the cumulative of k(·).

Custom scheme for modelling service times
The first step of modelling a M/G/2/∞ queue is to observe the real system and pick the empirical service times x i (Fig. 3a)). Here both of the arrival processes are considered to have exponential distributions. The empirical data will be used to calculate both: KDE and cumulative KDE functions. So the second step is to evaluate (8) and (10) by finding optimal h according to the empirical data. (8) will serve as the target distribution in Metropolis-Hastings method.
The key feature of the scheme presented is the proposal density q(·) for KDE. There are a number of techniques to construct it, because the problem is to find probability distribution similar in shape to the target distribution. In this paper the proposal density is considered to be a piecewise linear distribution, which acts as a histogram of KDE.
The researcher can divide the possible service times into intervals and set probabilities for them according to the areas below the KDE in those intervals. This is a fast KDE approximation, but we propose a different approach instead. Calculating the area below KDE in any interval involves approximating it by an angled line. We bypass this by calculating the cumulative KDE values at the endpoints of the intervals (Fig. 3b)) and thus knowing the values for proposal density. By doing this, the approximation of cumulative KDE is obtained and is used for sampling from the proposal distribution.  Next follows the sampling process itself. It is performed by simply drawing a number u i and reflecting it to the domain using obtained approximation of the cumulative KDE (which is also an angled line). The resulting value x i is accepted to the sample according to the (7).
The intervals between arrival times for each queue are generated by inverse CDF of exponential distribution. Now the implementation of the queuing system model in a chosen programming language remains.

The results of M/G/2/∞ system simulation
The queuing system with parameters λ 1 = 0.7, λ 2 = 2.2 and log-normally distributed service times has been modelled with the technique discussed above. The parameters of the service time distribution were µ logn = −1.5 and σ logn = 0.4. 5 evenly spaced values of service time x were used for constructing the approximation of the cumulative service time distribution function. By using the inverse log-normal distribution function, 100 empirical service times were generated.
According to the table above, the relative error of the modelling is unacceptable. This is due to the constant kernel width and the fact that KDE was performed on the support (0; ∞). Using variable KDE width would decrease relative errors to a certain degree. On the other hand, randomly sampled service times do not follow log-normal distribution identically. In other words, our approach simulates the queuing system as is, i.e. simulation is completely based upon empirical data.    4 shows the histograms of random service times which were generated by the proposed approach. Here the process of generating x i is divided into 2 parts. Firstly the proposal density is sampled. Then each x i is accepted to the sample or not. It is noticeable that MCMC is similar to well known rejection sampling.