{"id":375,"date":"2020-08-13T09:38:52","date_gmt":"2020-08-13T09:38:52","guid":{"rendered":"https:\/\/data-science.gotoauthority.com\/2020\/08\/13\/markov-chain-monte-carlo-simulation-for-airport-queuing-network\/"},"modified":"2020-08-13T09:38:52","modified_gmt":"2020-08-13T09:38:52","slug":"markov-chain-monte-carlo-simulation-for-airport-queuing-network","status":"publish","type":"post","link":"https:\/\/wealthrevelation.com\/data-science\/2020\/08\/13\/markov-chain-monte-carlo-simulation-for-airport-queuing-network\/","title":{"rendered":"Markov Chain Monte Carlo Simulation For Airport Queuing Network"},"content":{"rendered":"
\n
\"Markov\"Markov<\/span><\/div>\n
\n

Today we\u2019ll introduce a new set of algorithms called Markov Chain Monte Carlo which won’t fall in supervised learning algorithms<\/a>. In this blog post we will walk through, what Markov Chains are and where we can use it.<\/p>\n

We will introduce the main family of algorithms, known collectively as Markov Chain Monte Carlo (MCMC), that allows us to approximate the posterior distribution as calculated by Bayes’ Theorem<\/strong>.\u00a0<\/p>\n

In particular, we consider the Metropolis Algorithm<\/strong>, which is easily stated and relatively straightforward to understand. It serves as a useful starting point when learning about MCMC before delving into more sophisticated algorithms such as Metropolis-Hastings, Gibbs Samplers, and Hamiltonian Monte Carlo.<\/p>\n

Once we have described how MCMC works, we will carry it out using the open-source PyMC3 library<\/strong>, which takes care of many of the underlying implementation details, allowing us to concentrate on Bayesian modeling.<\/p>\n<\/div>\n

\n
\n
\n

Markov Chain Monte Carlo Simulation For Airport Queuing Network<\/p>\n<\/div>\n

\n\t\t\t
\n\t\t\t\t<\/i>
\n\t\t\t\tClick to Tweet<\/span>
\n\t\t\t<\/span>\n\t\t<\/p>\n<\/div>\n<\/div>\n

\n

Before we drive further let\u2019s look at what you are going to learn by the end of the article.<\/p>\n<\/div>\n

\n

Markov chain Monte Carlo analogy<\/h2>\n

Before getting started we\u2019ll try to understand the analogy behind Markov Chains. When we are getting into a learning curve in the field of analytics we have various divisions like first we\u2019ll start with forecasting and then linear regression after we\u2019ll get into classification algorithms<\/a> which are non-parametric models.\u00a0<\/p>\n

After this curve,\u00a0 we\u2019ll get into Neural Networks like CNN, R-CNN, Auto-encoders so on, and so forth.<\/p>\n

Once these are done now we will get into the stage of Markov Chains(MC) and Hidden Markov Chains (HMC) which are purely stochastic models in order to make a statement on random predictions.\u00a0<\/p>\n

Let\u2019s understand it by example,<\/p>\n

Decision tree algorithms<\/a> can say whether they are going to buy it or not. Random forest<\/a> says what are the various conditions need to be satisfied in order to make a statement. The logistic algorithm<\/a> can say yes\/no statements using sigmoid equations. When getting into CNN they can recognize the images and by using RNN they can do sequential tasking.<\/p>\n

Where we use Markov chains<\/h2>\n

All the above set of algorithms categories are intended to predict the feature using the historical data, But if we want to predict like if you\u2019re sitting in a restaurant and you\u2019re waiting for the waiter in order to take up the order right.<\/p>\n

So which algorithm be used in order to make the statement?.<\/p>\n

Whether Virat Kohli<\/strong> going to hit a six or not,\u00a0<\/p><\/blockquote>\n

which algorithm we need to use? In these scenarios we can\u2019t go into deep learning algorithms, we can\u2019t do forecasting right. So we want to make a statement on the above statements based on the current event we want to make a prediction about the next ball.\u00a0<\/p>\n

So for all these instant based predictions, the only route we have is Markov Chains and Hidden Markov Chains<\/strong>. Markov Chains are one of the powerful algorithms where we\u2019re able to extract or able to make a statement on random events.<\/p>\n

For these kinds of events, we prefer Markov Chains. Before we learn about markov chains, we need to learn about bayes’s rule.<\/p>\n

Let’s spend some time now.<\/p>\n

Bayes\u2019s Rule<\/h2>\n

If we recall Bayes\u2019s Rule:<\/p>\n

<\/a><\/p>\n

Where:<\/strong><\/p>\n

\nP(D) is evidence. This is the probability of the data as determined by summing (or integrating) across all possible values of
\n<\/a> , weighted by how strongly we believe in those particular values of
\n<\/a>.\n<\/div>\n

We can see that we need to calculate the evidence P(D). In order to achieve this we need to evaluate the following integral, which integrates over all possible values of, the parameters:<\/p>\n

<\/a><\/p>\n

The fundamental problem is that we are often unable to evaluate this integral analytically and so we must turn to a numerical approximation method instead.\u00a0<\/p>\n

An additional problem is that our models might require a large number of parameters. This means that our prior distributions could potentially have a large number of dimensions.<\/p>\n

This in turn means that our posterior distributions will also be high dimensional. Hence, we are in a situation where we have to numerically evaluate an integral in a potentially very large dimensional space.<\/p>\n

This means we are in a situation often described as the Curse of Dimensionality. Informally, this means that the volume of a high-dimensional space is so vast that any available data becomes extremely sparse within that space and hence leads to problems of statistical significance.<\/p>\n

Practically, in order to gain any statistical significance, the volume of data needed must grow exponentially with the number of dimensions.<\/p>\n

Such problems are often extremely difficult to tackle unless they are approached in an intelligent manner. The motivation behind Markov Chain Monte Carlo\u2019s methods is that they perform an intelligent search within a high dimensional space and thus Bayesian Models in high dimensions become easy to control.<\/p>\n

The basic idea is to sample from the posterior distribution by combining a \u201crandom search<\/strong>\u201d with a mechanism for intelligently \u201cjumping\u201d around, but in a manner that ultimately doesn\u2019t depend on where we started from.<\/p>\n

Hence Markov Chain Monte Carlo methods are memoryless searches performed with intelligent jumps.<\/p>\n

The Metropolis Algorithm\u00a0<\/h2>\n

There is a large family of Algorithms that perform MCMC. Most of these algorithms can be expressed at a high level as follows:\u00a0<\/p>\n

    \n
  1. Begin the algorithm at the current position in parameter space.\u00a0<\/li>\n
  2. Propose a \u201cjump\u201d to a new position in parameter space.\u00a0<\/li>\n
  3. Accept or reject the jump probabilistically using the prior information and available data.\u00a0<\/li>\n
  4. If the jump is accepted, move to the new position and return to step 1.<\/li>\n
  5. If the jumps are rejected, stay where you are and return to step 1.<\/li>\n
  6. After a set of a number of jumps has occurred, return all accepted positions.\u00a0<\/li>\n<\/ol>\n

    The main difference between MCMC algorithms<\/strong> occurs in how you jump as well as how you decide whether to jump.\u00a0<\/p>\n

    The Metropolis algorithm uses a normal distribution to propose a jump. This normal distribution has a mean value \u03bc<\/strong> which is equal to the current position and takes a “proposal width” for its standard deviation \u03c3.<\/strong><\/p>\n

    A normal distribution is a good choice for such a proposal distribution (for continuous parameters), it is more likely to select points nearer to the current position than further away. However, it will occasionally choose points further away, allowing the space to be explored.<\/p>\n

    Once the jump has been proposed, we need to decide (in a probabilistic<\/strong> manner) whether it is a good move to jump to the new position. How do we do this? We calculate the ratio of the proposal distribution of the new position and the proposal distribution at the current position to determine the probability of moving, p:<\/p>\n

    <\/a><\/p>\n

    About PyMC3<\/h3>\n

    PyMC3<\/a> is a Python package for Bayesian statistical modeling and probabilistic machine learning which focuses on advanced Markov chain Monte Carlo and variational fitting algorithms. It is a rewrite from scratch of the previous version of the PyMC software.<\/p>\n

    Simulation using PyMC3<\/h3>\n

    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.<\/p>\n

    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.<\/p>\n

    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.\u00a0<\/p>\n

    We can summarise<\/strong> this analysis with the following plot<\/p>\n<\/div>\n

    \"markov\"markov<\/span><\/div>\n
    \n

    Bayesian network representing the air traffic<\/strong> control problem<\/p>\n

    Considering our experience, we decide to model the random variables using the following distributions:<\/p>\n

      \n
    1. Passenger onboarding ~ Wald(\u00b5 = 0.5, \u03bb = 0.2)<\/li>\n
    2. Refueling ~ Wald( \u00b5 = 0.25, \u03bb = 0.5)<\/li>\n
    3. Departure traffic delay ~ Wald(\u00b5 = 0.1, \u03bb = 0.2)<\/li>\n
    4. Arrival traffic delay ~ Wald(\u00b5 = 0.1, \u03bb = 0.2)<\/li>\n
    5. Departure time = 12 + Departure traffic delay + max(Passenger onboarding, Refueling)<\/li>\n
    6. Rough weather ~ Bernoulli(p =0.35)<\/li>\n
    7. Flight time ~ Exponential(\u03bb = 0.5 – (0.1 . Rough weather))(The output of a Bernoulli distribution is 0 or 1 corresponding to False and True)<\/li>\n
    8. Arrival time = Departure time + Flight time + Arrival traffic delay<\/li>\n<\/ol>\n

      Departure time and Arrival time are functions of random variables, and the parameter \u03bb of Flight time is also a function of Rough Weather.\u00a0<\/p>\n

      Markov Chain Monte Carlo Simulation with PyMC3<\/h2>\n

      Even if the model is not very complex, the direct inference is rather inefficient, and therefore we want to simulate the process using PyMC3.<\/p>\n

      The first step is to create a model instance:<\/p>\n<\/div>\n

      \"pymc3\"pymc3<\/span><\/div>\n
      \n

      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<\/p>\n<\/div>\n

      \"Markov\"Markov<\/span><\/div>\n
      \n

      We have imported two namespaces, pymc3.distributions.continuous and pymc3.distributions.discrete because we are using both kinds of variables. <\/p>\n

      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. <\/p>\n

      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.<\/p>\n

      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. <\/p>\n

      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).<\/p>\n

      The other variables are very similar, except for flight_time, which is an exponential variable with a parameter \u03bb, which is a function of another variable (rough_weather). As a Bernoulli variable outputs 1 with probability p and 0 with probability 1 – p, \u03bb = 0.4 if there’s rough weather, and 0.5 otherwise.<\/p>\n

      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:<\/p>\n<\/div>\n

      \"mark\"mark<\/span><\/div>\n
      \n

      The output can be analyzed using the built-in pm.traceplot() function, which generates the plots for each of the sample\u2019s variables. The following graph shows the detail of one of them:<\/p>\n<\/div>\n

      \"arrive\"arrive<\/span><\/div>\n
      \n

      Distribution and samples for the arrival time random variable<\/strong>\u00a0<\/p>\n

      PyMC3 provides a statistical summary that can help us in making the right decisions using pm.summary()<\/strong>. In the following snippet, the output containing the summary of a single variable is shown:<\/p>\n<\/div>\n

      \"PyMC3\"PyMC3<\/span><\/div>\n
      \n

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

      This is only a very simple example to show the power of Bayesian networks.<\/p>\n

      List of parametric and non-parametric Algorithms<\/h2>\n

      Machine learning algorithms can be classified as two distinct groups: parametric and non-parametric<\/p>\n

      We can classify algorithms as non-parametric when models become more complex if the number of samples in the training set increases. Vice versa, a model would be parametric if the model becomes stable when the number of examples in the training set increases.\u00a0<\/p>\n

      In simple terms, we can say parametric has a functional form while non-parametric has no functional form.\u00a0<\/p>\n

      Functional form comprises a simple formula like y = f(x). So if you input a value, you are to get a fixed output value. It means if the data set is changed or being changed there is not much variation in the results. But in non-parametric algorithms, a small change in data sets can result in a large change in results.<\/p>\n