Day 25 - Sparse Candidate Algorithm/Starter/notebook.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-21932d3f3255749a",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "# Day 25  notebook\n",
                "\n",
                "The objectives of this notebook are to practice\n",
                "\n",
                "* computing the marginal distribution for a subset of variables\n",
                "* computing the Kullback\u2013Leibler divergence\n",
                "* selecting candidate parents for variables, as is done in the Sparse Candidate Algorithm"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-67cefe199f36d7c8",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Modules used for this assignment"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-badf267606df8231",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# standard library modules\n",
                "import math             # for log\n",
                "import random           # for seed\n",
                "import collections      # for Counter\n",
                "\n",
                "# course modules\n",
                "import bayesian_network # for the BayesianNetwork class"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-c6b1c490fb7cc4ce",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## Updates to the `bayesian_network` module\n",
                "\n",
                "The `BayesianNetwork` class has been updated a bit since last time.  Please note the following new methods:\n",
                "1. `joint_dist`: Returns the joint probability distribution represented by the network.\n",
                "2. `estimate_parameters`: Estimates (and sets) the parameters of the model given a data set using maximum likelihood.\n",
                "\n",
                "In addition, the `bayesian_network` module has a few functions for creating some of the example datasets, which we will use in this activity:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-3b1f130850d64952",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "scrolled": true,
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "lac_operon_network = bayesian_network.make_lac_operon_network()\n",
                "flight_weather_network = bayesian_network.make_flight_weather_network()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-b1e3cb08024c220b",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "We will use some simulated datasets from each of these example datasets:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-725cfd6321f50468",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "random.seed(5)\n",
                "lac_operon_dataset = [lac_operon_network.sample() for _ in range(1000)]\n",
                "random.seed(42)\n",
                "flight_weather_dataset = [flight_weather_network.sample() for _ in range(100)]"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-0c6a3ec2d3525e7a",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "In this activity we will need to work with the joint (and marginal) distributions estimated from these data sets, which we will refer to as the \"empirical\" distributions.  We will borrow a function from the Day 10 notebook that we used for estimating a joint distribution:"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-bd5ce9aee09f8b0c",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def estimate_joint_dist(observations):\n",
                "    \"\"\"Returns a joint distribution for the random variables measured in observations.\n",
                "    \n",
                "    Args:\n",
                "        observations: a list of tuples, where the ith element of each tuple represents the\n",
                "            observed value of the ith random variable.  \n",
                "    Returns:\n",
                "        A joint distribution in the form of a dictionary with random variable configurations\n",
                "        as keys and probabilities as values.\n",
                "    \"\"\"\n",
                "    counter = collections.Counter(observations)\n",
                "    return {key: count \/ len(observations) for key, count in counter.items()}"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-669d6127e5bbf340",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "lac_operon_dataset_joint = estimate_joint_dist(lac_operon_dataset)\n",
                "flight_weather_dataset_joint = estimate_joint_dist(flight_weather_dataset)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-af0d57bd4ad633c5",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 1: Computing the marginal distribution for a *subset* of variables (1 POINT)\n",
                "In the Day 10 notebook you implemented a functon for computing the marginal distribution of a single variable given a joint distribution.  In general, one can compute a marginal distribution for any (strict) subset of the variables in a joint distribution by summing out all of the variables that are not in that subset.  This will be important for determining the marginal distribution for a pair of random variablees in the Sparse Candidate Algorithm.  Implement the `compute_marginal_dist` function below."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def compute_marginal_dist(joint_distribution, indices):\n",
                "    \"\"\"Returns the marginal distribution a subset of random variables in a joint distribution.\n",
                "    \n",
                "    Args:\n",
                "        joint_distribution: a distribution in the form of a dictionary with with random variable\n",
                "            configurations (tuples) as keys and probabilities as values.\n",
                "        indices: a tuple or list of indices of the random variables to keep in the marginal distribution.\n",
                "    Returns:\n",
                "        The marginal distribution (using the same dictionary based representation as the input joint distribution)\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "compute_marginal_dist",
                    "locked": true,
                    "points": 1,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for compute_marginal_dist\n",
                "def round_dist(dist, digits=5):\n",
                "    \"\"\"Returns a new distribution with probabilities rounded to the specified number of digits.\"\"\"\n",
                "    return {key: round(value, digits) for key, value in dist.items()}\n",
                "\n",
                "# marginal distribution of weather (index 1) and status (index 2) variables in flight_weather system\n",
                "weather_status_dist = {\n",
                "    ('rain', 'delayed'): 0.159,\n",
                "    ('rain', 'on-time'): 0.141,\n",
                "    ('snow', 'delayed'): 0.174,\n",
                "    ('snow', 'on-time'): 0.026,\n",
                "    ('sun', 'delayed'): 0.085,\n",
                "    ('sun', 'on-time'): 0.415}\n",
                "assert round_dist(compute_marginal_dist(flight_weather_network.joint_dist(), (1, 2))) == weather_status_dist\n",
                "\n",
                "# marginal distribution of airline (index 0) and status (index 2) variables in flight_weather system\n",
                "airline_status_dist = {\n",
                "    ('Delta', 'delayed'): 0.117,\n",
                "    ('Delta', 'on-time'): 0.183,\n",
                "    ('United', 'delayed'): 0.301,\n",
                "    ('United', 'on-time'): 0.399}\n",
                "assert round_dist(compute_marginal_dist(flight_weather_network.joint_dist(), (0, 2))) == airline_status_dist\n",
                "\n",
                "# marginal distribution of status (index 2) variable in flight_weather system\n",
                "status_dist = {('delayed',): 0.418, ('on-time',): 0.582}\n",
                "assert round_dist(compute_marginal_dist(flight_weather_network.joint_dist(), (2,))) == status_dist\n",
                "\n",
                "# marginal distribution of L (index 0), G (index 2), and Z (index 3) variables in lac_operon_network system\n",
                "L_G_Z_dist = {\n",
                "    ('absent', 'absent', 'absent'): 0.3033,\n",
                "    ('absent', 'absent', 'high'): 0.09149,\n",
                "    ('absent', 'absent', 'low'): 0.05521,\n",
                "    ('absent', 'present', 'absent'): 0.3033,\n",
                "    ('absent', 'present', 'high'): 0.05577,\n",
                "    ('absent', 'present', 'low'): 0.09093,\n",
                "    ('present', 'absent', 'absent'): 0.0085,\n",
                "    ('present', 'absent', 'high'): 0.03083,\n",
                "    ('present', 'absent', 'low'): 0.01067,\n",
                "    ('present', 'present', 'absent'): 0.0085,\n",
                "    ('present', 'present', 'high'): 0.01099,\n",
                "    ('present', 'present', 'low'): 0.03052}\n",
                "\n",
                "assert round_dist(compute_marginal_dist(lac_operon_network.joint_dist(), (0, 2, 6))) == L_G_Z_dist\n",
                "\n",
                "print(\"SUCCESS: compute_marginal_dist passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-c12725b4785da3d3",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 2: Computing the Kullback\u2013Leibler divergence (1 POINT)\n",
                "The Sparse Candidate algorithm uses Kullback\u2013Leibler divergence to identify candidate parents for each random variable.  Implement the computation of the Kullback\u2013Leibler divergence in the function `kl_divergence` below."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def kl_divergence(p, q):\n",
                "    \"\"\"Computes the Kullback\u2013Leibler divergence from Q to P, i.e., D_KL(P || Q)\n",
                "    \n",
                "    Args:\n",
                "        p and q: distributions over the same set of random variables.  Each distribution is \n",
                "                 represented as a dictionary with configurations (tuples) as keys and probabilities\n",
                "                 as values.\n",
                "    Returns:\n",
                "        The Kullback\u2013Leibler divergence as a floating point value.\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "kl_divergence",
                    "locked": true,
                    "points": 1,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# tests for kl_divergence\n",
                "true_weather_status_dist = {\n",
                "    ('rain', 'delayed'): 0.159,\n",
                "    ('rain', 'on-time'): 0.141,\n",
                "    ('snow', 'delayed'): 0.174,\n",
                "    ('snow', 'on-time'): 0.026,\n",
                "    ('sun', 'delayed'): 0.085,\n",
                "    ('sun', 'on-time'): 0.415}\n",
                "empirical_weather_status_dist = {\n",
                "    ('rain', 'delayed'): 0.07,\n",
                "    ('rain', 'on-time'): 0.2,\n",
                "    ('snow', 'delayed'): 0.19,\n",
                "    ('snow', 'on-time'): 0.08,\n",
                "    ('sun', 'delayed'): 0.1,\n",
                "    ('sun', 'on-time'): 0.36}\n",
                "assert round(kl_divergence(true_weather_status_dist, empirical_weather_status_dist), 5) == 0.08182\n",
                "\n",
                "true_status_dist = {('delayed',): 0.418, ('on-time',): 0.582}\n",
                "empirical_status_dist = {('delayed',): 0.36, ('on-time',): 0.64}\n",
                "assert round(kl_divergence(true_status_dist, empirical_status_dist), 5) == 0.00715\n",
                "\n",
                "true_L_G_Z_dist = {\n",
                "    ('absent', 'absent', 'absent'): 0.3033,\n",
                "    ('absent', 'absent', 'high'): 0.09149,\n",
                "    ('absent', 'absent', 'low'): 0.05521,\n",
                "    ('absent', 'present', 'absent'): 0.3033,\n",
                "    ('absent', 'present', 'high'): 0.05577,\n",
                "    ('absent', 'present', 'low'): 0.09093,\n",
                "    ('present', 'absent', 'absent'): 0.0085,\n",
                "    ('present', 'absent', 'high'): 0.03083,\n",
                "    ('present', 'absent', 'low'): 0.01067,\n",
                "    ('present', 'present', 'absent'): 0.0085,\n",
                "    ('present', 'present', 'high'): 0.01099,\n",
                "    ('present', 'present', 'low'): 0.03052}\n",
                "empirical_L_G_Z_dist = {\n",
                "    ('absent', 'absent', 'absent'): 0.321,\n",
                "    ('absent', 'absent', 'high'): 0.091,\n",
                "    ('absent', 'absent', 'low'): 0.058,\n",
                "    ('absent', 'present', 'absent'): 0.287,\n",
                "    ('absent', 'present', 'high'): 0.069,\n",
                "    ('absent', 'present', 'low'): 0.072,\n",
                "    ('present', 'absent', 'absent'): 0.009,\n",
                "    ('present', 'absent', 'high'): 0.038,\n",
                "    ('present', 'absent', 'low'): 0.005,\n",
                "    ('present', 'present', 'absent'): 0.01,\n",
                "    ('present', 'present', 'high'): 0.012,\n",
                "    ('present', 'present', 'low'): 0.028}\n",
                "assert round(kl_divergence(true_L_G_Z_dist, empirical_L_G_Z_dist), 5) == 0.00811\n",
                "\n",
                "print(\"SUCCESS: kl_divergence passed all tests!\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-5494c7b262697e99",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "source": [
                "## PROBLEM 3: Selecting candidate parents for variables (1 POINT)\n",
                "\n",
                "Using your `compute_marginal_dist` and `kl_divergence` functions, we can now implement the \"restrict\" step of the Sparse Candidate algorithm.  In this step, the parameters of the current network are estimated using the dataset, and then the empirical and network pairwise marginal distributions are compared with each other using KL-divergence.  Then a set of $k$ candidate parents are selected for each variable, which includes the current parents and the other variables for which the KL-divergence is highest.  Most of the implementation is provided for you, except for `candidate_parents`, which selects the candidate parents for a single variable.  Fill in the implementation for this function."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": false,
                    "grade_id": "cell-bc9b1afddf079179",
                    "locked": true,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "def matrix(num_rows, num_cols, value=None):\n",
                "    \"\"\"Constructs a matrix (a list of lists)\"\"\"\n",
                "    return [[value] * num_cols for i in range(num_rows)]\n",
                "\n",
                "def pair_marginal_divergences(empirical_joint_dist, network_joint_dist):\n",
                "    \"\"\"For each pair of variables, computes the KL-divergence between the empirical and network\n",
                "       derived marginal distributions for that pair of variables.\n",
                "\n",
                "    Args: \n",
                "        empirical_joint_dist: the empirical joint distribution (dictionary representation)\n",
                "        network_joint_dist:  the network joint distribution (dictionary representation)\n",
                "    Returns:\n",
                "        A matrix (list of lists) M with M[i][j] giving the KL-divergence between empirical and network\n",
                "        derived marginal distributions of variables i and j\n",
                "    \"\"\"\n",
                "    # the length of each configuration is the # of variables \n",
                "    # (we'll grab one configuration by taking first key when iterating over the empirical_joint_dist)\n",
                "    num_variables = len(next(iter(empirical_joint_dist))) \n",
                "    M = matrix(num_variables, num_variables, 0)\n",
                "    for i in range(num_variables):\n",
                "        for j in range(i + 1, num_variables):\n",
                "            empirical_marginal = compute_marginal_dist(empirical_joint_dist, (i, j))\n",
                "            network_marginal = compute_marginal_dist(network_joint_dist, (i, j))\n",
                "            M[i][j] = M[j][i] = kl_divergence(empirical_marginal, network_marginal)\n",
                "    return M\n",
                "\n",
                "def restrict(network, dataset, k):\n",
                "    \"\"\"Computes the candidate parents for all variables in the network, given a dataset.\n",
                "    \n",
                "    Args:\n",
                "        network: The Bayesian network\n",
                "        dataset: a list of observations (tuples)\n",
                "        k: the number of candidate parents to select for each variable\n",
                "    Returns:\n",
                "        A list of sets of the candidate parents for each variable.\n",
                "    \"\"\"\n",
                "    # compute maximum likelihood (ML) estimates for the network for the given data\n",
                "    network.estimate_parameters(dataset)\n",
                "    \n",
                "    # compute the joint distribution represented by the network, given the ML parameters\n",
                "    network_joint_dist = network.joint_dist()\n",
                "    \n",
                "    # compute the empirical joint distribution from the data\n",
                "    empirical_joint_dist = estimate_joint_dist(dataset)\n",
                "    \n",
                "    # For each pair of variables, compute the KL-divergence between the empirical and network\n",
                "    # marginal distributions for that pair of variables\n",
                "    M = pair_marginal_divergences(empirical_joint_dist, network_joint_dist)\n",
                "    \n",
                "    # Return a list of the candidate parent sets for each variable\n",
                "    return [candidate_parents(network, i, M[i], k) for i in range(network.num_vertices())]"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true
            },
            "outputs": [],
            "source": [
                "def candidate_parents(network, i, divergences, k):\n",
                "    \"\"\"Computes the k candidate parents for variable i in the network, given the computed divergences.\n",
                "    \n",
                "    The candidate parent set will always include the current parents of variable i.  Other variables will\n",
                "    be added to the candidate set based on the divergences.\n",
                "    Args:\n",
                "        network: The Bayesian network\n",
                "        i: the index of the variable for which to select parents\n",
                "        divergences: a list of the pairwise marginal divergences between variable i and all other variables\n",
                "        k: the number of candidate parents to select\n",
                "    Returns:\n",
                "        A set of the candidate parent variable indices.\n",
                "    \"\"\"\n",
                "    ###\n",
                "    ### YOUR CODE HERE\n",
                "    ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {
                "collapsed": true,
                "nbgrader": {
                    "grade": true,
                    "grade_id": "candidate_parents",
                    "locked": true,
                    "points": 1,
                    "schema_version": 1,
                    "solution": false
                },
                "editable": false,
                "deletable": false
            },
            "outputs": [],
            "source": [
                "# test candidate_parents\n",
                "empty_lac_operon_network = bayesian_network.make_empty_lac_operon_network()\n",
                "assert (restrict(empty_lac_operon_network, lac_operon_dataset, 1) == \n",
                "        [{4}, {4}, {5}, {5}, {6}, {2}, {4}])\n",
                "assert (restrict(empty_lac_operon_network, lac_operon_dataset, 2) ==\n",
                "        [{4, 6}, {4, 6}, {5, 6}, {5, 6}, {0, 6}, {2, 3}, {1, 4}])\n",
                "\n",
                "assert (restrict(lac_operon_network, lac_operon_dataset, 2) ==\n",
                "        [{5, 6}, {2, 6}, {1, 6}, {0, 6}, {0, 1}, {2, 3}, {4, 5}])\n",
                "assert (restrict(lac_operon_network, lac_operon_dataset, 3) ==\n",
                "        [{3, 5, 6}, {2, 5, 6}, {1, 4, 6}, {0, 4, 6}, {0, 1, 5}, {2, 3, 4}, {1, 4, 5}])\n",
                "\n",
                "print(\"SUCCESS: candidate_parents passed all tests!\")"
            ]
        }
    ],
    "metadata": {
        "celltoolbar": "Create Assignment",
        "kernelspec": {
            "display_name": "Python 3 [3.6]",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text\/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.6.4"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}