LChrisman/Queuing theory application of MetaLog

This page contains brainstorming-level investigation into the application of the Keelin (MetaLog) distribution to queueing theory. Right now it is just a place for notes, may contain mistakes, etc.

Motivation

The paper

• Baucells, Chrisman, Keelin and Xu (2025) , "On the properties of the Metalog distribution"

derives a closed-form formula for partial expectation.

Here I wonder whether this could be applied to a model of an M/G/1 queue, where service times have a MetaLog distribution, enabling a closed-form expression for partial expectations on the wait time distribution.

M/G/1 queue

Arrivals follow a Poisson process with an average arrival rate of [math]\displaystyle{ \lambda }[/math] (i.e., arrival rate is Exponential(λ)).

Service times are MetaLog distributed. And there is a single server.

Notation:

• [math]\displaystyle{ \lambda }[/math]: Arrival rate. Reciprocal of expected time between successive arrivals.
• [math]\displaystyle{ \mu }[/math]: Expected number of service completions per unit time. It is the reciprocal of mean service time.
• [math]\displaystyle{ \rho = \lambda/\mu }[/math]: The utilization.
• [math]\displaystyle{ W }[/math] = Mean wait time
• [math]\displaystyle{ L }[/math] = Mean # of customers in the system ([math]\displaystyle{ L = \lambda W }[/math]).
• [math]\displaystyle{ x }[/math] = Variable used for service-time.
• [math]\displaystyle{ p_S(x) }[/math] = PDF of service-time distribution (note: positive only).
• [math]\displaystyle{ F_S(x) }[/math] = CDF of the service-time distribution.
• [math]\displaystyle{ M_S(y) }[/math] = Quantile function of service-time distribution. (A MetaLog)

For a functional system, [math]\displaystyle{ \rho \lt 1 }[/math]. If [math]\displaystyle{ \rho\ge 1 }[/math] then jobs arrive faster than they are serviced.

The Laplace-Stieltjes transform of [math]\displaystyle{ M_S(y) }[/math] is

[math]\displaystyle{ \mathcal{L}\{M_S\}(s) = \int_0^1 e^{- s M_S(y)} dy = \int_{-\infty}^\infty e^{-s x} dF_S(x) = E[e^{-s x}] }[/math]

From the Pollaczek-Khinchine transform equation:

[math]\displaystyle{ \mathcal{L}\{W\} = {{(1-p) s \mathcal{L}\{p\}(s)} \over {s - \lambda (1-\mathcal{L}\{p\}(s))}} = {{(1-p)s} \over {s - \lambda + \lambda \mathcal{L}\{M_S \}(s)}} }[/math]

Waiting time distribution as a meta-log

The idea now is to fit a MetaLog to the inverse transform of [math]\displaystyle{ \mathcal{L}\{W\} }[/math], giving us a MetaLog representation of wait time that could then be used to compute partial expectations on waiting time.

Standard methods for inverting: Talbot, Stehfest, Fourier-series methods. Straightforward because [math]\displaystyle{ \mathcal{L}\{W\} }[/math] is known explicitly.

Alternative: Simple MC simulation of service times.

Possible advantages:

• Analytic moments
• Analytic expression for partial expectation.
• Easy random variate generation
• Faster sensitivity studies
• When service-time distribution is a rational function, [math]\displaystyle{ \mathcal{L}\{W\}(s) }[/math] is rational and its inverse is a finite mixture of exponentials. A MetaLog can match this exactly only with an infinite number of terms. Shows that an exact "solution" for the inverse as a MetaLog does not exist.

Laplace transforms of MetaLog

Two-term

[math]\displaystyle{ M(p) = a_0 + a_1 \ell(p) }[/math]

[math]\displaystyle{ \mathcal{L}\{M\}(s) = a_0 + {a_1 \over 2}\left( \phi\left({s\over 2}\right) - \phi\left({{s+1}\over 2}\right)\right) }[/math]

where [math]\displaystyle{ \phi(z) = {{\Gamma'(z)}\over{\Gamma(z)}} }[/math] is the digamma function.

k-terms

[math]\displaystyle{ M(y) = \sum_{i=1}^k a_i B_i(y) }[/math]

Note: Paper uses [math]\displaystyle{ Y_i(y) }[/math] where I use [math]\displaystyle{ B_i(y) }[/math] for the basis functions.

[math]\displaystyle{ \mathcal{L}\{M\}(s) = \sum_{i=1}^k a_i \mathcal{L}\{B_i\}(s) }[/math]

The poly terms

When [math]\displaystyle{ i \in \{1,4,5,8,9,...\} }[/math],

[math]\displaystyle{ B_i(y) = a_i (y - 0.5)^d }[/math] where [math]\displaystyle{ d = \left\lfloor{{i-1}\over 2}\right\rfloor }[/math].

Then

[math]\displaystyle{ \mathcal{L}\{B_i(y)\}(s) = \sum_{j=0}^d \left( \begin{array}{c}d \\ j\end{array} \right) \left( -{1\over 2}\right)^{m-j} {{\gamma( j+1,s)} \over {s^{j+1}}} }[/math]

where [math]\displaystyle{ \gamma(s,x) }[/math] is the incomplete gamma function.

The logit terms

When [math]\displaystyle{ i \in \{2, 3, 6, 7, 10, 11, ...\} }[/math],

[math]\displaystyle{ B_i(y) = a_i (y - 0.5)^d logit(y) }[/math]

Then

[math]\displaystyle{ \mathcal{L}\{B_i(y)\}(s) = \sum_{j=0}^m \left( \begin{array}{c}d \\ j\end{array} \right) \left( -{1\over 2}\right)^{m-j} (-1)^j {{d^j}\over{ds^j}} \left[ {{e^{-s}-1}\over s} ( \ln s + \gamma_E ) - {{E_1(s) + e^{-s}E_i(s)}\over s} \right] }[/math]

where

• [math]\displaystyle{ \gamma_E }[/math] is Euler's constant (0.5772156649...)
• [math]\displaystyle{ E_1(s) = \int_s^\infty {{e^{-t}}\over t} dt }[/math] is the exponential integral.
• [math]\displaystyle{ Ei(s) = -E_i(-s) }[/math] is the standard exponential integral.