Learning_Bayesian_Networks.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Learning Bayesian Networks.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/kristiyandd/ai-combinator/blob/master/Learning_Bayesian_Networks.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pr7Z3IA9npXC",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 202
        },
        "outputId": "64d7f0a3-f8af-4a49-d4e9-4e3afde48866"
      },
      "source": [
        "!pip install pgmpy\n",
        "!pip install pandas\n",
        "!pip install numpy"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: pgmpy in /usr/local/lib/python3.6/dist-packages (0.1.7)\n",
            "Requirement already satisfied: scipy>=1.0.0 in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.3.0)\n",
            "Requirement already satisfied: networkx<1.12,>=1.11 in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.11)\n",
            "Requirement already satisfied: numpy>=1.14.0 in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.16.4)\n",
            "Requirement already satisfied: decorator>=3.4.0 in /usr/local/lib/python3.6/dist-packages (from networkx<1.12,>=1.11->pgmpy) (4.4.0)\n",
            "Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (0.24.2)\n",
            "Requirement already satisfied: python-dateutil>=2.5.0 in /usr/local/lib/python3.6/dist-packages (from pandas) (2.5.3)\n",
            "Requirement already satisfied: pytz>=2011k in /usr/local/lib/python3.6/dist-packages (from pandas) (2018.9)\n",
            "Requirement already satisfied: numpy>=1.12.0 in /usr/local/lib/python3.6/dist-packages (from pandas) (1.16.4)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.6/dist-packages (from python-dateutil>=2.5.0->pandas) (1.12.0)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (1.16.4)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GzNfOBdPiINs",
        "colab_type": "text"
      },
      "source": [
        "# **Learning Bayesian Networks**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "haooIWMtiSQP",
        "colab_type": "text"
      },
      "source": [
        "Previous notebooks showed how Bayesian networks economically encode a probability distribution over a set of variables, and how they can be used e.g. to predict variable states, or to generate new samples from the joint distribution. This section will be about obtaining a Bayesian network, given a set of sample data. Learning a Bayesian network can be split into two problems:\n",
        "\n",
        "  **Parameter learning**: Given a set of data samples and a DAG that captures the dependencies between the variables, estimate the (conditional) probability distributions of the individual variables.\n",
        "\n",
        "  **Structure learning**: Given a set of data samples, estimate a DAG that captures the dependencies between the variables.\n",
        "\n",
        "  This notebook aims to illustrate how parameter learning and structure learning can be done with pgmpy. Currently, the library supports:\n",
        "  \n",
        "\n",
        "Parameter learning for discrete nodes:\n",
        "\n",
        "*   Maximum Likelihood Estimation\n",
        "*   Bayesian Estimation\n",
        "    \n",
        "Structure learning for discrete, fully observed networks:\n",
        "    \n",
        "*    Score-based structure estimation (BIC/BDeu/K2 score; exhaustive search, hill climb/tabu search)\n",
        "*   Constraint-based structure estimation (PC)\n",
        "*   Hybrid structure estimation (MMHC)\n",
        "\n",
        "\n",
        "# **Parameter Learning **\n",
        "\n",
        "  Suppose we have the following data:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rmqcjyXdiH32",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 269
        },
        "outputId": "98da23ef-673c-41a5-dcf0-1c0a305b1c1b"
      },
      "source": [
        "\n",
        "import pandas as pd\n",
        "data = pd.DataFrame(data={'fruit': [\"banana\", \"apple\", \"banana\", \"apple\", \"banana\",\"apple\", \"banana\", \n",
        "                                    \"apple\", \"apple\", \"apple\", \"banana\", \"banana\", \"apple\", \"banana\",], \n",
        "                          'tasty': [\"yes\", \"no\", \"yes\", \"yes\", \"yes\", \"yes\", \"yes\", \n",
        "                                    \"yes\", \"yes\", \"yes\", \"yes\", \"no\", \"no\", \"no\"], \n",
        "                          'size': [\"large\", \"large\", \"large\", \"small\", \"large\", \"large\", \"large\",\n",
        "                                    \"small\", \"large\", \"large\", \"large\", \"large\", \"small\", \"small\"]})\n",
        "print(data)"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "     fruit tasty   size\n",
            "0   banana   yes  large\n",
            "1    apple    no  large\n",
            "2   banana   yes  large\n",
            "3    apple   yes  small\n",
            "4   banana   yes  large\n",
            "5    apple   yes  large\n",
            "6   banana   yes  large\n",
            "7    apple   yes  small\n",
            "8    apple   yes  large\n",
            "9    apple   yes  large\n",
            "10  banana   yes  large\n",
            "11  banana    no  large\n",
            "12   apple    no  small\n",
            "13  banana    no  small\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PyZuFK0ljuxU",
        "colab_type": "text"
      },
      "source": [
        "We know that the variables relate as follows:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-qb8XKJIh1fO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "\n",
        "from pgmpy.models import BayesianModel\n",
        "\n",
        "model = BayesianModel([('fruit', 'tasty'), ('size', 'tasty')])  # fruit -> tasty <- size"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gJuxmhzGj0Ht",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "Parameter learning is the task to estimate the values of the conditional probability distributions (CPDs), for the variables fruit, size, and tasty."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2WglsabviF7O",
        "colab_type": "text"
      },
      "source": [
        "State counts\n",
        "To make sense of the given data, we can start by counting how often each state of the variable occurs. If the variable is dependent on parents, the counts are done conditionally on the parents states, i.e. for seperately for each parent configuration:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7Jbq8VxNj6Mq",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 185
        },
        "outputId": "5f2a5e53-edc2-4944-eeb5-9ea1fb0eb4ac"
      },
      "source": [
        "from pgmpy.estimators import ParameterEstimator\n",
        "pe = ParameterEstimator(model, data)\n",
        "print(\"\\n\", pe.state_counts('fruit'))  # unconditional\n",
        "print(\"\\n\", pe.state_counts('tasty'))  # conditional on fruit and size"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "\n",
            "         fruit\n",
            "apple       7\n",
            "banana      7\n",
            "\n",
            " fruit apple       banana      \n",
            "size  large small  large small\n",
            "tasty                         \n",
            "no      1.0   1.0    1.0   1.0\n",
            "yes     3.0   2.0    5.0   0.0\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v93TdPWdj4AY",
        "colab_type": "text"
      },
      "source": [
        "We can see, for example, that as many apples as bananas were observed and that 5 large bananas were tasty, while only 1 was not."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VaeRKm7qj-SS",
        "colab_type": "text"
      },
      "source": [
        "**Maximum Likelihood Estimation**\n",
        "\n",
        "\n",
        "A natural estimate for the CPDs is to simply use the relative frequencies, with which the variable states have occured. We observed 7 apples among a total of 14 fruits, so we might guess that about 50% of fruits are apples.\n",
        "\n",
        "This approach is Maximum Likelihood Estimation (MLE). According to MLE, we should fill the CPDs in such a way, that $P(\\text{data}|\\text{model})$ is maximal. This is achieved when using the relative frequencies. See [1], section 17.1 for an introduction to ML parameter estimation. pgmpy supports MLE as follows:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "njSaSDoFkcvK",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 252
        },
        "outputId": "fbcf3c73-1f8e-495e-9fdf-c27736f1a8f1"
      },
      "source": [
        "from pgmpy.estimators import MaximumLikelihoodEstimator\n",
        "mle = MaximumLikelihoodEstimator(model, data)\n",
        "print(mle.estimate_cpd('fruit'))  # unconditional\n",
        "print(mle.estimate_cpd('tasty'))  # conditional"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "+---------------+-----+\n",
            "| fruit(apple)  | 0.5 |\n",
            "+---------------+-----+\n",
            "| fruit(banana) | 0.5 |\n",
            "+---------------+-----+\n",
            "+------------+--------------+--------------------+---------------------+---------------+\n",
            "| fruit      | fruit(apple) | fruit(apple)       | fruit(banana)       | fruit(banana) |\n",
            "+------------+--------------+--------------------+---------------------+---------------+\n",
            "| size       | size(large)  | size(small)        | size(large)         | size(small)   |\n",
            "+------------+--------------+--------------------+---------------------+---------------+\n",
            "| tasty(no)  | 0.25         | 0.3333333333333333 | 0.16666666666666666 | 1.0           |\n",
            "+------------+--------------+--------------------+---------------------+---------------+\n",
            "| tasty(yes) | 0.75         | 0.6666666666666666 | 0.8333333333333334  | 0.0           |\n",
            "+------------+--------------+--------------------+---------------------+---------------+\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q0ey1LYUkeZs",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "mle.estimate_cpd(variable) computes the state counts and divides each cell by the (conditional) sample size. The mle.get_parameters()-method returns a list of CPDs for all variable of the model.\n",
        "\n",
        "The built-in fit()-method of BayesianModel provides more convenient access to parameter estimators:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aaWwb3bjkh5f",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 67
        },
        "outputId": "b97fc3f8-d7c3-4df9-9e6a-32e8d525f602"
      },
      "source": [
        "# Calibrate all CPDs of `model` using MLE:\n",
        "model.fit(data, estimator=MaximumLikelihoodEstimator)"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "WARNING:root:Replacing existing CPD for fruit\n",
            "WARNING:root:Replacing existing CPD for size\n",
            "WARNING:root:Replacing existing CPD for tasty\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SkmI3Kn_ki8K",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "While very straightforward, the ML estimator has the problem of overfitting to the data. In above CPD, the probability of a large banana being tasty is estimated at 0.833, because 5 out of 6 observed large bananas were tasty. Fine. But note that the probability of a small banana being tasty is estimated at 0.0, because we observed only one small banana and it happened to be not tasty. But that should hardly make us certain that small bananas aren't tasty! We simply do not have enough observations to rely on the observed frequencies. If the observed data is not representative for the underlying distribution, ML estimations will be extremly far off.\n",
        "\n",
        "When estimating parameters for Bayesian networks, lack of data is a frequent problem. Even if the total sample size is very large, the fact that state counts are done conditionally for each parents configuration causes immense fragmentation. If a variable has 3 parents that can each take 10 states, then state counts will be done seperately for 10^3 = 1000 parents configurations. This makes MLE very fragile and unstable for learning Bayesian Network parameters. A way to mitigate MLE's overfitting is Bayesian Parameter Estimation."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "24DHWylrkm4q",
        "colab_type": "text"
      },
      "source": [
        "**Bayesian Parameter Estimation**\n",
        "\n",
        "\n",
        "The Bayesian Parameter Estimator starts with already existing prior CPDs, that express our beliefs about the variables before the data was observed. Those \"priors\" are then updated, using the state counts from the observed data. See [1], Section 17.3 for a general introduction to Bayesian estimators.\n",
        "\n",
        "One can think of the priors as consisting in pseudo state counts, that are added to the actual counts before normalization. Unless one wants to encode specific beliefs about the distributions of the variables, one commonly chooses uniform priors, i.e. ones that deem all states equiprobable.\n",
        "\n",
        "A very simple prior is the so-called K2 prior, which simply adds 1 to the count of every single state. A somewhat more sensible choice of prior is BDeu (Bayesian Dirichlet equivalent uniform prior). For BDeu we need to specify an equivalent sample size N and then the pseudo-counts are the equivalent of having observed N uniform samples of each variable (and each parent configuration). In pgmpy:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5Ehug-m1ksU_",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 168
        },
        "outputId": "bc0ba2b8-dbd4-4fe5-d029-d5d246c3e6ec"
      },
      "source": [
        "from pgmpy.estimators import BayesianEstimator\n",
        "est = BayesianEstimator(model, data)\n",
        "\n",
        "print(est.estimate_cpd('tasty', prior_type='BDeu', equivalent_sample_size=10))"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "+------------+---------------------+--------------------+--------------------+---------------------+\n",
            "| fruit      | fruit(apple)        | fruit(apple)       | fruit(banana)      | fruit(banana)       |\n",
            "+------------+---------------------+--------------------+--------------------+---------------------+\n",
            "| size       | size(large)         | size(small)        | size(large)        | size(small)         |\n",
            "+------------+---------------------+--------------------+--------------------+---------------------+\n",
            "| tasty(no)  | 0.34615384615384615 | 0.4090909090909091 | 0.2647058823529412 | 0.6428571428571429  |\n",
            "+------------+---------------------+--------------------+--------------------+---------------------+\n",
            "| tasty(yes) | 0.6538461538461539  | 0.5909090909090909 | 0.7352941176470589 | 0.35714285714285715 |\n",
            "+------------+---------------------+--------------------+--------------------+---------------------+\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5rODmDHiku43",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "The estimated values in the CPDs are now more conservative. In particular, the estimate for a small banana being not tasty is now around 0.64 rather than 1.0. Setting equivalent_sample_size to 10 means that for each parent configuration, we add the equivalent of 10 uniform samples (here: +5 small bananas that are tasty and +5 that aren't).\n",
        "\n",
        "BayesianEstimator, too, can be used via the fit()-method. Full example:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lO8azF9pkv8H",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 487
        },
        "outputId": "02354b73-4b83-4dcf-ceba-feb0dfd8bdfb"
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "from pgmpy.models import BayesianModel\n",
        "from pgmpy.estimators import BayesianEstimator\n",
        "\n",
        "# generate data\n",
        "data = pd.DataFrame(np.random.randint(low=0, high=2, size=(5000, 4)), columns=['A', 'B', 'C', 'D'])\n",
        "model = BayesianModel([('A', 'B'), ('A', 'C'), ('D', 'C'), ('B', 'D')])\n",
        "\n",
        "model.fit(data, estimator=BayesianEstimator, prior_type=\"BDeu\") # default equivalent_sample_size=5\n",
        "for cpd in model.get_cpds():\n",
        "    print(cpd)"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "+------+----------+\n",
            "| A(0) | 0.497203 |\n",
            "+------+----------+\n",
            "| A(1) | 0.502797 |\n",
            "+------+----------+\n",
            "+------+---------------------+--------------------+\n",
            "| A    | A(0)                | A(1)               |\n",
            "+------+---------------------+--------------------+\n",
            "| B(0) | 0.5056258790436006  | 0.4852970395390423 |\n",
            "+------+---------------------+--------------------+\n",
            "| B(1) | 0.49437412095639943 | 0.5147029604609576 |\n",
            "+------+---------------------+--------------------+\n",
            "+------+--------------------+--------------------+---------------------+---------------------+\n",
            "| A    | A(0)               | A(0)               | A(1)                | A(1)                |\n",
            "+------+--------------------+--------------------+---------------------+---------------------+\n",
            "| D    | D(0)               | D(1)               | D(0)                | D(1)                |\n",
            "+------+--------------------+--------------------+---------------------+---------------------+\n",
            "| C(0) | 0.4804440822978203 | 0.4865213082259663 | 0.5069091648039017  | 0.49805636540330417 |\n",
            "+------+--------------------+--------------------+---------------------+---------------------+\n",
            "| C(1) | 0.5195559177021797 | 0.5134786917740337 | 0.49309083519609836 | 0.5019436345966958  |\n",
            "+------+--------------------+--------------------+---------------------+---------------------+\n",
            "+------+--------------------+--------------------+\n",
            "| B    | B(0)               | B(1)               |\n",
            "+------+--------------------+--------------------+\n",
            "| D(0) | 0.4788263762855414 | 0.5029697089685211 |\n",
            "+------+--------------------+--------------------+\n",
            "| D(1) | 0.5211736237144585 | 0.4970302910314789 |\n",
            "+------+--------------------+--------------------+\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_T9pPn62k0OQ",
        "colab_type": "text"
      },
      "source": [
        "# **Structure Learning**\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "To learn model structure (a DAG) from a data set, there are two broad techniques:\n",
        "\n",
        "*   score-based structure learning\n",
        "*   constraint-based structure learning\n",
        "\n",
        "The combination of both techniques allows further improvement:\n",
        "\n",
        "*   hybrid structure learning\n",
        "\n",
        "\n",
        "We briefly discuss all approaches and give examples.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-N3IL_H2lF8m",
        "colab_type": "text"
      },
      "source": [
        "# **Score-based Structure Learning**\n",
        "\n",
        "\n",
        "This approach construes model selection as an optimization task. It has two building blocks:\n",
        "\n",
        "A scoring function $s_D\\colon M \\to \\mathbb R$ that maps models to a numerical score, based on how well they fit to a given data set $D$.\n",
        "A search strategy to traverse the search space of possible models $M$ and select a model with optimal score.\n",
        "\n",
        "\n",
        "**Scoring functions**\n",
        "\n",
        "\n",
        "Commonly used scores to measure the fit between model and data are Bayesian Dirichlet scores such as BDeu or K2 and the Bayesian Information Criterion (BIC, also called MDL). See [1], Section 18.3 for a detailed introduction on scores. As before, BDeu is dependent on an equivalent sample size."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4jheNi_tlOVL",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 118
        },
        "outputId": "ce697c1b-a450-4eb2-aa14-958c2ca5a2dd"
      },
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "from pgmpy.estimators import BdeuScore, K2Score, BicScore\n",
        "from pgmpy.models import BayesianModel\n",
        "\n",
        "# create random data sample with 3 variables, where Z is dependent on X, Y:\n",
        "data = pd.DataFrame(np.random.randint(0, 4, size=(5000, 2)), columns=list('XY'))\n",
        "data['Z'] = data['X'] + data['Y']\n",
        "\n",
        "bdeu = BdeuScore(data, equivalent_sample_size=5)\n",
        "k2 = K2Score(data)\n",
        "bic = BicScore(data)\n",
        "\n",
        "model1 = BayesianModel([('X', 'Z'), ('Y', 'Z')])  # X -> Z <- Y\n",
        "model2 = BayesianModel([('X', 'Z'), ('X', 'Y')])  # Y <- X -> Z\n",
        "\n",
        "\n",
        "print(bdeu.score(model1))\n",
        "print(k2.score(model1))\n",
        "print(bic.score(model1))\n",
        "\n",
        "print(bdeu.score(model2))\n",
        "print(k2.score(model2))\n",
        "print(bic.score(model2))"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "-13937.222029340224\n",
            "-14328.010361283588\n",
            "-14293.26130380845\n",
            "-20901.894717858864\n",
            "-20928.729232000325\n",
            "-20945.941054985335\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hPcM5iOFlQWa",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "While the scores vary slightly, we can see that the correct model1 has a much higher score than model2. Importantly, these scores decompose, i.e. they can be computed locally for each of the variables given their potential parents, independent of other parts of the network:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CW1ZAIGDlSlF",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 67
        },
        "outputId": "bd080f2d-093a-4b8d-aaa4-3dd31fe33975"
      },
      "source": [
        "\n",
        "print(bdeu.local_score('Z', parents=[]))\n",
        "print(bdeu.local_score('Z', parents=['X']))\n",
        "print(bdeu.local_score('Z', parents=['X', 'Y']))"
      ],
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "-9188.60772547207\n",
            "-6991.873950785468\n",
            "-57.11920936734987\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2KEzCraWlUPg",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "**Search strategies**\n",
        "\n",
        "\n",
        "The search space of DAGs is super-exponential in the number of variables and the above scoring functions allow for local maxima. The first property makes exhaustive search intractable for all but very small networks, the second prohibits efficient local optimization algorithms to always find the optimal structure. Thus, identifiying the ideal structure is often not tractable. Despite these bad news, heuristic search strategies often yields good results.\n",
        "\n",
        "If only few nodes are involved (read: less than 5), ExhaustiveSearch can be used to compute the score for every DAG and returns the best-scoring one:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TyjbrUfGlX4u",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 487
        },
        "outputId": "287a0e9f-481d-40d9-9abf-1b321ebf423f"
      },
      "source": [
        "from pgmpy.estimators import ExhaustiveSearch\n",
        "\n",
        "es = ExhaustiveSearch(data, scoring_method=bic)\n",
        "best_model = es.estimate()\n",
        "print(best_model.edges())\n",
        "\n",
        "print(\"\\nAll DAGs by score:\")\n",
        "for score, dag in reversed(es.all_scores()):\n",
        "    print(score, dag.edges())"
      ],
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[('X', 'Z'), ('Y', 'Z')]\n",
            "\n",
            "All DAGs by score:\n",
            "-14293.26130380845 [('X', 'Z'), ('Y', 'Z')]\n",
            "-14327.280577499805 [('Y', 'X'), ('Z', 'X'), ('Z', 'Y')]\n",
            "-14327.280577499805 [('Y', 'Z'), ('Y', 'X'), ('Z', 'X')]\n",
            "-14327.280577499805 [('X', 'Z'), ('Y', 'Z'), ('Y', 'X')]\n",
            "-14327.280577499805 [('X', 'Y'), ('Z', 'X'), ('Z', 'Y')]\n",
            "-14327.280577499807 [('X', 'Y'), ('X', 'Z'), ('Z', 'Y')]\n",
            "-14327.280577499807 [('X', 'Y'), ('X', 'Z'), ('Y', 'Z')]\n",
            "-16492.280585496126 [('X', 'Y'), ('Z', 'Y')]\n",
            "-16492.57119054296 [('Y', 'X'), ('Z', 'X')]\n",
            "-18746.631168250824 [('Z', 'X'), ('Z', 'Y')]\n",
            "-18746.631168250824 [('Y', 'Z'), ('Z', 'X')]\n",
            "-18746.631168250824 [('X', 'Z'), ('Z', 'Y')]\n",
            "-20911.631176247138 [('Z', 'Y')]\n",
            "-20911.63117624714 [('Y', 'Z')]\n",
            "-20911.921781293975 [('Z', 'X')]\n",
            "-20911.921781293975 [('X', 'Z')]\n",
            "-20945.650449938497 [('Y', 'X'), ('Z', 'Y')]\n",
            "-20945.650449938497 [('Y', 'Z'), ('Y', 'X')]\n",
            "-20945.650449938497 [('X', 'Y'), ('Y', 'Z')]\n",
            "-20945.941054985335 [('X', 'Z'), ('Y', 'X')]\n",
            "-20945.941054985335 [('X', 'Y'), ('Z', 'X')]\n",
            "-20945.941054985335 [('X', 'Y'), ('X', 'Z')]\n",
            "-23076.921789290296 []\n",
            "-23110.94106298165 [('Y', 'X')]\n",
            "-23110.941062981652 [('X', 'Y')]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HCj58MtblaLe",
        "colab_type": "text"
      },
      "source": [
        "Once more nodes are involved, one needs to switch to heuristic search. HillClimbSearch implements a greedy local search that starts from the DAG start (default: disconnected DAG) and proceeds by iteratively performing single-edge manipulations that maximally increase the score. The search terminates once a local maximum is found."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M5smiEZ1lc1y",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "408bf3cb-1c6f-4580-ec8f-294696f4cb4b"
      },
      "source": [
        "from pgmpy.estimators import HillClimbSearch\n",
        "\n",
        "# create some data with dependencies\n",
        "data = pd.DataFrame(np.random.randint(0, 3, size=(2500, 8)), columns=list('ABCDEFGH'))\n",
        "data['A'] += data['B'] + data['C']\n",
        "data['H'] = data['G'] - data['A']\n",
        "\n",
        "hc = HillClimbSearch(data, scoring_method=BicScore(data))\n",
        "best_model = hc.estimate()\n",
        "print(best_model.edges())"
      ],
      "execution_count": 20,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[('A', 'H'), ('A', 'C'), ('A', 'B'), ('C', 'B'), ('G', 'H')]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zoOdyHyGledk",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "[('A', 'C'), ('A', 'B'), ('C', 'B'), ('G', 'A'), ('G', 'H'), ('H', 'A')]\n",
        "\n",
        "\n",
        "The search correctly identifies e.g. that B and C do not influnce H directly, only through A and of course that D, E, F are independent.\n",
        "\n",
        "\n",
        "\n",
        "To enforce a wider exploration of the search space, the search can be enhanced with a tabu list. The list keeps track of the last n modfications; those are then not allowed to be reversed, regardless of the score. Additionally a white_list or black_list can be supplied to restrict the search to a particular subset or to exclude certain edges. The parameter max_indegree allows to restrict the maximum number of parents for each node."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AwSeYYmNlhK9",
        "colab_type": "text"
      },
      "source": [
        "**Constraint-based Structure Learning**\n",
        "\n",
        "A different, but quite straightforward approach to build a DAG from data is this:\n",
        "\n",
        "  1.Identify independencies in the data set using hypothesis tests\n",
        "  \n",
        "  2.Construct DAG (pattern) according to identified independencies\n",
        "\n",
        "**(Conditional) Independence Tests**\n",
        "\n",
        "Independencies in the data can be identified using chi2 conditional independence tests. To this end, constraint-based estimators in pgmpy have a test_conditional_independence(X, Y, Zs)-method, that performs a hypothesis test on the data sample. It allows to check if X is independent from Y given a set of variables Zs:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "31N1SwT8lsSL",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 101
        },
        "outputId": "3c5b1cc6-9ff5-4149-f083-8d50ed259797"
      },
      "source": [
        "from pgmpy.estimators import ConstraintBasedEstimator\n",
        "\n",
        "data = pd.DataFrame(np.random.randint(0, 3, size=(2500, 8)), columns=list('ABCDEFGH'))\n",
        "data['A'] += data['B'] + data['C']\n",
        "data['H'] = data['G'] - data['A']\n",
        "data['E'] *= data['F']\n",
        "\n",
        "est = ConstraintBasedEstimator(data)\n",
        "\n",
        "print(est.test_conditional_independence('B', 'H'))          # dependent\n",
        "print(est.test_conditional_independence('B', 'E'))          # independent\n",
        "print(est.test_conditional_independence('B', 'H', ['A']))   # independent\n",
        "print(est.test_conditional_independence('A', 'G'))          # independent\n",
        "print(est.test_conditional_independence('A', 'G',  ['H']))  # dependent"
      ],
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "(618.1365846887918, 9.790437115287228e-114, True)\n",
            "(6.08642201164528, 0.8675360809549586, True)\n",
            "(8.716141373222557, 0.9999999983732946, True)\n",
            "(11.514738671308502, 0.9317668416134064, True)\n",
            "(4619.0, 0.0, True)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_YCLXxVolzJR",
        "colab_type": "text"
      },
      "source": [
        "test_conditional_independence() returns a tripel (chi2, p_value, sufficient_data), consisting in the computed chi2 test statistic, the p_value of the test, and a heuristig flag that indicates if the sample size was sufficient. The p_value is the probability of observing the computed chi2 statistic (or an even higher chi2 value), given the null hypothesis that X and Y are independent given Zs.\n",
        "\n",
        "\n",
        "This can be used to make independence judgements, at a given level of significance:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JqFYK5pkl2Wj",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 101
        },
        "outputId": "86052eaa-762c-4d55-ffb0-1ee59c59a064"
      },
      "source": [
        "\n",
        "def is_independent(X, Y, Zs=[], significance_level=0.05):\n",
        "    return est.test_conditional_independence(X, Y, Zs)[1] >= significance_level\n",
        "\n",
        "print(is_independent('B', 'H'))\n",
        "print(is_independent('B', 'E'))\n",
        "print(is_independent('B', 'H', ['A']))\n",
        "print(is_independent('A', 'G'))\n",
        "print(is_independent('A', 'G', ['H']))"
      ],
      "execution_count": 22,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "False\n",
            "True\n",
            "True\n",
            "True\n",
            "False\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mQGv89s5l5e2",
        "colab_type": "text"
      },
      "source": [
        "**DAG (pattern) construction**\n",
        "\n",
        "\n",
        "With a method for independence testing at hand, we can construct a DAG from the data set in three steps:\n",
        "\n",
        "1.Construct an undirected skeleton - estimate_skeleton()\n",
        "\n",
        "2.Orient compelled edges to obtain partially directed acyclid graph (PDAG; I-equivalence class of DAGs) - skeleton_to_pdag()\n",
        "\n",
        "3.Extend DAG pattern to a DAG by conservatively orienting the remaining edges in some way - pdag_to_dag()\n",
        "\n",
        "\n",
        "Step 1.&2. form the so-called PC algorithm, see [2], page 550. PDAGs are DirectedGraphs, that may contain both-way edges, to indicate that the orientation for the edge is not determined."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GudZn7FImFzK",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 67
        },
        "outputId": "3bef4566-edfe-4b3a-f6a7-e04118a868ce"
      },
      "source": [
        "skel, seperating_sets = est.estimate_skeleton(significance_level=0.01)\n",
        "print(\"Undirected edges: \", skel.edges())\n",
        "\n",
        "pdag = est.skeleton_to_pdag(skel, seperating_sets)\n",
        "print(\"PDAG edges:       \", pdag.edges())\n",
        "\n",
        "model = est.pdag_to_dag(pdag)\n",
        "print(\"DAG edges:        \", model.edges())"
      ],
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Undirected edges:  [('A', 'B'), ('A', 'C'), ('A', 'H'), ('E', 'F'), ('G', 'H')]\n",
            "PDAG edges:        [('A', 'H'), ('B', 'A'), ('C', 'A'), ('E', 'F'), ('F', 'E'), ('G', 'H')]\n",
            "DAG edges:         [('A', 'H'), ('B', 'A'), ('C', 'A'), ('F', 'E'), ('G', 'H')]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qBopQo5-mIHn",
        "colab_type": "text"
      },
      "source": [
        "The estimate()-method provides a shorthand for the three steps above and directly returns a BayesianModel:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WK5CBokcmJse",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "8cb99d2d-317e-4776-df98-dd829a8c54fa"
      },
      "source": [
        "print(est.estimate(significance_level=0.01).edges())"
      ],
      "execution_count": 24,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[('A', 'H'), ('B', 'A'), ('C', 'A'), ('F', 'E'), ('G', 'H')]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fEHCbBzHmKok",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "The estimate_from_independencies()-method can be used to construct a BayesianModel from a provided set of independencies (see class documentation for further features & methods):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vMWvNsJbmOb6",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "10d0c3e8-5846-47e2-e0f9-2fce84bc291a"
      },
      "source": [
        "from pgmpy.independencies import Independencies\n",
        "\n",
        "ind = Independencies(['B', 'C'],\n",
        "                     ['A', ['B', 'C'], 'D'])\n",
        "ind = ind.closure()  # required (!) for faithfulness\n",
        "\n",
        "model = ConstraintBasedEstimator.estimate_from_independencies(\"ABCD\", ind)\n",
        "\n",
        "print(model.edges())"
      ],
      "execution_count": 25,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[('A', 'D'), ('B', 'D'), ('C', 'D')]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wfQnqikzmNjd",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "PC PDAG construction is only guaranteed to work under the assumption that the identified set of independencies is faithful, i.e. there exists a DAG that exactly corresponds to it. Spurious dependencies in the data set can cause the reported independencies to violate faithfulness. It can happen that the estimated PDAG does not have any faithful completions (i.e. edge orientations that do not introduce new v-structures). In that case a warning is issued."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fmP0iyuHmSGN",
        "colab_type": "text"
      },
      "source": [
        "# **Hybrid Structure Learning**\n",
        "The MMHC algorithm [3] combines the constraint-based and score-based method. It has two parts:\n",
        "\n",
        "1.Learn undirected graph skeleton using the constraint-based construction procedure MMPC\n",
        "\n",
        "2.Orient edges using score-based optimization (BDeu score + modified hill-climbing)\n",
        "\n",
        "We can perform the two steps seperately, more or less as follows:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hZKTfLEJmZmL",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "outputId": "bf31f6d0-3fa1-4ee2-a31c-2bfc571f2ef7"
      },
      "source": [
        "from pgmpy.estimators import MmhcEstimator\n",
        "\n",
        "data = pd.DataFrame(np.random.randint(0, 3, size=(2500, 8)), columns=list('ABCDEFGH'))\n",
        "data['A'] += data['B'] + data['C']\n",
        "data['H'] = data['G'] - data['A']\n",
        "data['E'] *= data['F']\n",
        "\n",
        "mmhc = MmhcEstimator(data)\n",
        "skeleton = mmhc.mmpc()\n",
        "print(\"Part 1) Skeleton: \", skeleton.edges())\n",
        "\n",
        "# use hill climb search to orient the edges:\n",
        "hc = HillClimbSearch(data, scoring_method=BDeuScore(data))\n",
        "model = hc.estimate(tabu=10, white_list=skeleton.to_directed().edges())\n",
        "print(\"Part 2) Model:    \", model.edges())"
      ],
      "execution_count": 26,
      "outputs": [
        {
          "output_type": "error",
          "ename": "ImportError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-26-f0422d536974>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mpgmpy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mestimators\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mMmhcEstimator\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2500\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ABCDEFGH'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'A'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'B'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'C'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'H'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'G'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'A'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mImportError\u001b[0m: cannot import name 'MmhcEstimator'",
            "",
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0;32m\nNOTE: If your import is failing due to a missing package, you can\nmanually install dependencies using either !pip or !apt.\n\nTo view examples of installing some common dependencies, click the\n\"Open Examples\" button below.\n\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kVqZN1ZVmdi7",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "MmhcEstimator.estimate() is a shorthand for both steps and directly estimates a BayesianModel.\n",
        "\n",
        "# **Conclusion**\n",
        "\n",
        "This notebook aimed to give an overview of pgmpy's estimators for learning Bayesian network structure and parameters. For more information about the individual functions see their docstring documentation. If you used pgmpy's structure learning features to satisfactorily learn a non-trivial network from real data, feel free to drop us an eMail via the mailing list or just open a Github issue. We'd like to put your network in the examples-section!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P-EzBV39mi8z",
        "colab_type": "text"
      },
      "source": [
        "# **References**\n",
        "[1] Koller & Friedman, Probabilistic Graphical Models - Principles and Techniques, 2009\n",
        "\n",
        "[2] Neapolitan, Learning Bayesian Networks, 2003\n",
        "\n",
        "[3] Tsamardinos et al., The max-min hill-climbing BN structure learning algorithm, 2005"
      ]
    }
  ]
}