Bayesian Queueing Networks?

I am interested in using PyMC to estimate parameters in queueing network models. To date I have found the following tools outside PyMC:

  • The queue standard library has queue.Queue, queue.LifoQueue, and queue.PriorityQueue. Combined with a tool to handle a network of queues, like Networkx, this
  • "CiW is a discrete event simulation library for open queueing networks. Its core features include the capability to simulate networks of queues, multiple customer classes, and implementation of Type I blocking for restricted networks. A number of other features are also implemented, including priorities, baulking, reneging, schedules, batch arrivals, service disciplines, dynamic customer classes, and deadlock detection.
  • Queueing-tool is a ready-to-use queueing network library.
  • SimPy is a discrete event simulation tool with an emphasis on asynchronous even handling.

None of the above tools are Bayesian per se. My guess is that using queue.PriorityQueue might be the easiest path for integration with PyMC b/c it does not have any 3rd party dependencies. SimPy is likewise built on CPython with few dependencies. However, CiW and Queuing-tool are already ready-made implementations of queue networks.

Is there a way to implement PyMC queueing networks or patch existing tools to work with PyMC?

This blog post provides a minimal example of a Markov chain in PyMC3. However, modern queueing networks are often more detailed than just the Markovian state transitions between states. See queueing theory for an overview.

Objective: Implementing Markov Chains model
Author: Venkatesh Nagilla
Date: 2020-08-09

import pymc3 as pm
model = pm.Model()

import pymc3.distributions.continuous as pmc
import pymc3.distributions.discrete as pmd
import pymc3.math as pmm

with model:
    passenger_onboarding = pmc.Wald('Passenger Onboarding', mu=0.5, lam=0.2)
    refueling = pmc.Wald('Refueling', mu=0.25, lam=0.5)
    departure_traffic_delay = pmc.Wald('Departure Traffic Delay', mu=0.1, lam=0.2)

    departure_time = pm.Deterministic('Departure Time',
                                      12.0 + departure_traffic_delay +
                                      pmm.switch(passenger_onboarding >= refueling,

    rough_weather = pmd.Bernoulli('Rough Weather', p=0.35)

    flight_time = pmc.Exponential('Flight Time', lam=0.5 - (0.1 * rough_weather))
    arrival_traffic_delay = pmc.Wald('Arrival Traffic Delay', mu=0.1, lam=0.2)

    arrival_time = pm.Deterministic('Arrival time',
                                    departure_time +
                                    flight_time +

nb_samples = 500
with model:
    samples = pm.sample(draws=nb_samples, random_seed=1000)