Sampling example using PyMC3

PyMC3 is a powerful Python Bayesian framework that relies on Theano to perform high-speed computations (see the information box at the end of this paragraph for the installation instructions). It implements all the most important continuous and discrete distributions, and performs the sampling process mainly using the No-U-Turn and Metropolis-Hastings algorithms. For all the details about the API (distributions, functions, and plotting utilities), I suggest visiting the documentation home page http://docs.pymc.io/index.html, where it's also possible to find some very intuitive tutorials.

The example we want to model and simulate is based on this scenario: a daily flight from London to Rome has a scheduled departure time at 12:00 am, and a standard flight time of two hours. We need to organize the operations at the destination airport, but we don't want to allocate resources when the plane hasn't landed yet. Therefore, we want to model the process using a Bayesian network and considering some common factors that can influence the arrival time. In particular, we know that the onboarding process can be longer than expected, as well as the refueling one, even if they are carried out in parallel. London air traffic control can also impose a delay, and the same can happen when the plane is approaching Rome. We also know that the presence of rough weather can cause another delay due to a change of route. We can summarize this analysis with the following plot:

Considering our experience, we decide to model the random variables using the following distributions:

• Passenger onboarding ∼ Wald(μ=0.5, λ=0.2)
• Refueling ∼ Wald(μ=0.25, λ=0.5)
• Departure traffic delay ∼ Wald(μ=0.1, λ=0.2)
• Arrival traffic delay ∼ Wald(μ=0.1, λ=0.2)
• Departure time = 12 + Departure traffic delay + max(Passenger onboarding, Refueling)
• Rough weather ∼ Bernoulli(p=0.35)
• Flight time ∼ Exponential(λ=0.5 - (0.1 · Rough weather)) (The output of a Bernoulli distribution is 0 or 1 corresponding to False and True)
• Arrival time = Departure time + Flight time + Arrival traffic delay

The probability density functions are:

Departure Time and Arrival Time are functions of random variables, and the parameter λ of Flight Time is also a function of Rough Weather. 

Even if the model is not very complex, the direct inference is rather inefficient, and therefore we want to simulate the process using PyMC3.

The first step is to create a model instance:

import pymc3 as pm
 
 model = pm.Model()

From now on, all operations must be performed using the context manager provided by the model variable. We can now set up all the random variables of our Bayesian network:

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, 
                                                 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 + 
                                    arrival_traffic_delay)

We have imported two namespaces, pymc3.distributions.continuous and pymc3.distributions.discrete, because we are using both kinds of variable. Wald and exponential are continuous distributions, while Bernoulli is discrete. In the first three rows, we declare the variables passenger_onboarding, refueling, and departure_traffic_delay. The structure is always the same: we need to specify the class corresponding to the desired distribution, passing the name of the variable and all the required parameters.

The departure_time variable is declared as pm.Deterministic. In PyMC3, this means that, once all the random elements have been set, its value becomes completely determined. Indeed, if we sample from departure_traffic_delay, passenger_onboarding, and refueling, we get a determined value for departure_time. In this declaration, we've also used the utility function pmm.switch, which operates a binary choice based on its first parameter (for example, if A > B, return A, else return B).

The other variables are very similar, except for flight_time, which is an exponential variable with a parameter λ, which is a function of another variable (rough_weather). As a Bernoulli variable outputs 1 with probability p and 0 with probability 1 - p, λ = 0.4 if there's rough weather, and 0.5 otherwise.

Once the model has been set up, it's possible to simulate it through a sampling process. PyMC3 picks the best sampler automatically, according to the type of variables. As the model is not very complex, we can limit the process to 500 samples:

nb_samples = 500

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

The output can be analyzed using the built-in pm.traceplot() function, which generates the plots for each of the sample's variables. The following graph shows the detail of one of them:

The right column shown the samples generated for the random variable (in this case, the arrival time), while the left column shows the relative frequencies. This plot can be useful to have a visual confirmation of our initial ideas; in fact, the arrival time has the majority of its mass concentrated in the interval 14:00 to 16:00 (the numbers are always decimal, so it's necessary to convert the times); however, we should integrate to get the probabilities. Instead, through the pm.summary() function, PyMC3 provides a statistical summary that can help us in making the right decisions. In the following snippet, the output containing the summary of a single variable is shown:

pm.summary(samples)

...

Arrival time:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  15.174           2.670            0.102            [12.174, 20.484]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  12.492         13.459         14.419         16.073         22.557

For each variable, it contains mean, standard deviation, Monte Carlo error, 95% highest posterior density interval, and the posterior quantiles. In our case, we know that the plane will land at about 15:10 (15.174).

This is only a very simple example to show the power of Bayesian networks. For deep insight, I suggest the book Introduction to Statistical Decision Theory, Pratt J., Raiffa H., Schlaifer R., The MIT Press, where it's possible to study different Bayesian applications that are out of the scope of this book.