class16_bayesnet_python3_template.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class 16:  Bayesian Networks\n",
    "## Infer a Bayesian network from a matrix of discretized phospho-flow cytometry data.\n",
    "### Based on supplementary data from the 2005 article by Karen Sachs et al. (Science v308, 2005). \n",
    "\n",
    "In this class exercise, we will use the fundamental theorem for the likelihood of a Bayesian network structure for categorical variables, in order to score the posterior probability of the network shown in the Sachs et al. article (Figure 3A) vs. the phospho-flow cytometry data that the same authors provided in their supplementary data. The phospho-flow cytometry data have been already discretized for you (see \"class16_bayesnet_dataprep_R.ipynb\").  We will need to implement a single-vertex log-likelihood function using Theorem 1 from the article by Cooper & Herskovits in *Machine Learning* (volume 9, pages 309-347, 1992)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this notebook, we will need a bunch of python modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import pandas\n",
    "import math\n",
    "import scipy.special\n",
    "import igraph\n",
    "import itertools\n",
    "import matplotlib.pyplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the tab-delimited data file of discretized phosphoprotein expression data (12 columns; first 11 columns are the expression levels -- \"low\", \"medium\", \"high\"; last column is the experiment identifier for the row; there are nine experiments). Print out the first six lines of the data frame, so you can see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>praf</th>\n",
       "      <th>pmek</th>\n",
       "      <th>plcg</th>\n",
       "      <th>PIP2</th>\n",
       "      <th>PIP3</th>\n",
       "      <th>p44.42</th>\n",
       "      <th>pakts473</th>\n",
       "      <th>PKA</th>\n",
       "      <th>PKC</th>\n",
       "      <th>P38</th>\n",
       "      <th>pjnk</th>\n",
       "      <th>expt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>high</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>low</td>\n",
       "      <td>low</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>medium</td>\n",
       "      <td>high</td>\n",
       "      <td>high</td>\n",
       "      <td>1_cd3cd28</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     praf    pmek    plcg PIP2    PIP3  p44.42 pakts473     PKA     PKC  \\\n",
       "0     low     low     low  low    high     low      low  medium  medium   \n",
       "1     low     low  medium  low     low  medium   medium  medium     low   \n",
       "2  medium    high  medium  low  medium  medium   medium  medium  medium   \n",
       "3  medium    high  medium  low     low     low      low  medium  medium   \n",
       "4     low  medium     low  low  medium  medium   medium     low     low   \n",
       "5     low     low  medium  low     low  medium   medium  medium  medium   \n",
       "\n",
       "      P38    pjnk       expt  \n",
       "0    high    high  1_cd3cd28  \n",
       "1     low    high  1_cd3cd28  \n",
       "2  medium  medium  1_cd3cd28  \n",
       "3  medium  medium  1_cd3cd28  \n",
       "4  medium    high  1_cd3cd28  \n",
       "5    high    high  1_cd3cd28  "
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g_discret_data = \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a log likelihood function for the vertex degree distribution. Compute a vector `degarray_float` of floating point values for the degree of each vertex. The log likelihood is the logarithm of the product (over all vertices) of the single-vertex power-law degree distribution, which we will define as 1/(1+k)^2. Remember that the log of the product is the sum of the logs, so you can just compute this as `numpy.sum` of `numpy.log` of `numpy.power(1.0 + degarray_float, -2)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_prob_network_prior(network):\n",
    "    degarray = \n",
    "    degarray_float = \n",
    "    degarray_float[:] = \n",
    "    return FILL IN HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a vectorized log-factorial function, since it doesn't seem to be a builtin in python. Use the fact that factorial(n) = gamma(n+1), where `gamma` means the gamma function. Note that we can compute log gamma as `scipy.special.gammaln`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lfactorial(n):\n",
    "    return FILL IN HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a log likelihood function for a single vertex, based on Theorem 1 in the 1992 article in *Machine Learning* by Cooper & Herskovits. Note:  we are using igraph's adjacency matrix format which is the transpose of Newman's adjacency matrix definition!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood_network_vertex(network, vertex, discret_data):\n",
    "    # network is a NxN numpy matrix (N = 11 vertices)\n",
    "    # vertex is an integer\n",
    "    # discret_data is \"g_discret_data\" (N = 11 vertices, M = 7466 samples)\n",
    "    parents_vertex = numpy.where(network[:,vertex]==1)[0].tolist()\n",
    "    all_vertices = parents_vertex.copy()\n",
    "    all_vertices.append(vertex)\n",
    "    df1 = discret_data.iloc[:, all_vertices] \n",
    "\n",
    "    # change the name of the vertex column to \"vertex\"\n",
    "    df1_column_names = df1.columns.tolist()\n",
    "    df1_column_names[len(parents_vertex)] = \"vertex\"\n",
    "    df1.columns = df1_column_names\n",
    "\n",
    "    # count the data, grouped by all columns (parents & vertex)\n",
    "    df1 = df1.groupby(df1.columns.values.tolist()).size().reset_index(name='count')\n",
    "    \n",
    "    # make a new column called \"count factorial\" that is the log of the factorial of the count\n",
    "    df1[\"countfactorial\"] = lfactorial(df1[\"count\"].values)\n",
    "    \n",
    "    # drop the \"count\" column, as we no longer need it\n",
    "    df1 = df1.drop(\"count\", 1)\n",
    "    \n",
    "    if len(parents_vertex) > 0:\n",
    "\n",
    "        # sum up log-factorial-counts values for all possible states of \"vertex\" and its parent vertices, \n",
    "        # for each possible combination of parent vertices\n",
    "        nijkdf = df1.groupby(by=df1.columns[list(range(0,len(parents_vertex)))].tolist(), \n",
    "                             as_index=False).sum()\n",
    "        \n",
    "        # count number of cells with each possible combination of states for its parent vertices\n",
    "        df3 = discret_data.iloc[:, parents_vertex]\n",
    "        nijdf = df3.groupby(df3.columns.values.tolist(), as_index=False).size().reset_index()\n",
    "        nijdf_col_names = nijdf.columns.values\n",
    "        nijdf_col_names[len(nijdf_col_names)-1] = \"count\"\n",
    "        nijdf.columns = nijdf_col_names\n",
    "        \n",
    "        # compute the log factorial of the counts\n",
    "        nijdf[\"countfactorial\"] = math.log(2) - lfactorial(2 + nijdf[\"count\"])\n",
    "        \n",
    "        # drop the \"count\" column as we no longer need it\n",
    "        nijdf = nijdf.drop(\"count\", 1)\n",
    "        \n",
    "        # merge the two log-factorial-count values from nijdf and nijkdf, into two columns in a single dataframe\n",
    "        nmerge = nijdf.merge(nijkdf, how=\"outer\", \n",
    "                     on=nijkdf.columns[0:(len(nijkdf.columns)-1)].values.tolist(),\n",
    "                     copy=False)\n",
    "        \n",
    "        # sum the log-factorial-count values from nijdf and nijkdf\n",
    "        llh_res = numpy.sum(nmerge[\"countfactorial_x\"] + nmerge[\"countfactorial_y\"])\n",
    "    else:\n",
    "        # we handle the case of no parent vertices specially, to simplify the code\n",
    "        M = discret_data.shape[0]\n",
    "        llh_res = math.log(2) - lfactorial(M + 2) + numpy.sum(df1[\"countfactorial\"].values)\n",
    "    return llh_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_posterior_prob_network(network, discret_data):\n",
    "    Nvert = network.shape[1]\n",
    "    lpvert_values = []\n",
    "    for i in range(0, Nvert):\n",
    "        lpvert_value = log_likelihood_network_vertex(network, i, discret_data)\n",
    "        lpvert_values.append(lpvert_value)\n",
    "    return COMPUTE THE SUM OF THE LOG-PRIOR-PROBABILITY AND THE LOG-LIKELIHOOD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a log-posterior-probability function for the whole graph, using the per-vertex likelihood and the network structure prior:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define an adjacency matrix for the \"real\" network shown in Fig. 3A of the Sachs et al. article (not including the \"missed\" edges which are the dotted arcs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "real_network_adj_np = numpy.zeros(shape=[11,11])\n",
    "molec_names = list(g_discret_data.columns)[0:11]\n",
    "real_network_adj_df = pandas.DataFrame(real_network_adj_np, index=molec_names, columns=molec_names)\n",
    "real_network_adj_df[\"PKA\"][\"PKC\"] = 1\n",
    "real_network_adj_df[\"praf\"][\"PKC\"] = 1\n",
    "real_network_adj_df[\"pjnk\"][\"PKC\"] = 1\n",
    "real_network_adj_df[\"P38\"][\"PKC\"] = 1\n",
    "real_network_adj_df[\"pjnk\"][\"PKA\"] = 1\n",
    "real_network_adj_df[\"P38\"][\"PKA\"] = 1\n",
    "real_network_adj_df[\"praf\"][\"PKA\"] = 1\n",
    "real_network_adj_df[\"pmek\"][\"PKA\"] = 1\n",
    "real_network_adj_df[\"p44.42\"][\"PKA\"] = 1  # p44.42 = ERK\n",
    "real_network_adj_df[\"pakts473\"][\"PKA\"] = 1\n",
    "real_network_adj_df[\"pakts473\"][\"p44.42\"] = 1\n",
    "real_network_adj_df[\"pmek\"][\"PKC\"] = 1\n",
    "real_network_adj_df[\"pmek\"][\"praf\"] = 1\n",
    "real_network_adj_df[\"p44.42\"][\"pmek\"] = 1\n",
    "real_network_adj_df[\"PIP2\"][\"plcg\"] = 1\n",
    "real_network_adj_df[\"PIP3\"][\"plcg\"] = 1\n",
    "real_network_adj_df[\"PIP2\"][\"PIP3\"] = 1\n",
    "real_network_adj = real_network_adj_df.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make an igraph network out of the adjacency matrix that you just created, and print the network summary and plot the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IGRAPH D--- 11 17 -- '"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the log posterior probability of the real network `real_network_adj` from Sachs *et al.* Figure 3A, using your function `log_posterior_prob_network`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-84131.03111133902\n"
     ]
    }
   ],
   "source": [
    "lp_real = \n",
    "print(lp_real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate 1000 random rewirings of the network -- eliminating any rewired digraphs that contain cycles -- and for each randomly rewired network compute the log posterior probability. Use `igraph.Graph.rewire_edges(prob=1, loops=False, multiple=False)` to randomly rewire the graph. Make sure to make a copy of the graph before you rewire it. Check that the rewired graph is a DAG using `igraph.Graph.is_dag()`. If it is a DAG, compute the log posterior probability of the network using `log_posterior_prob_network`. You can get the adjacency matrix for the random graph as a numpy array using `numpy.array(randgraph.get_adjacency().data)`  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "lprobs_rand = []\n",
    "for _ in itertools.repeat(None, 1000):\n",
    "    graphcopy = \n",
    "    REWIRE EDGES HERE\n",
    "    if CHECK FOR IS_DAG\n",
    "        APPEND LOG-POSTERIOR-PROBABILITY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram the log ratio of the \"real\" network's posterior probability to the posterior probabilities of each of the random networks. Does it appear that the published network is a reasonable approximation to the maximum *a posteriori* (MAP) estimate?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAYaUlEQVR4nO3de7QlZX3m8e8joICigHSwh0salAXBiSK2BBZeAkYF8R5HZRzFy9gZo+OYGCNoVtQ1yQSjiZd4QRQVjRoVARnQKIKXMUvBBhpoBAShVbCFxqCAEhD8zR/1HtycPqd7n+5TZ3d3fT9r7XWq3qra9VbtfZ5d+91Vb6WqkCQNx30mXQFJ0sIy+CVpYAx+SRoYg1+SBsbgl6SB2XrSFRjHLrvsUkuWLJl0NSRps3LBBRfcVFWLppdvFsG/ZMkSli9fPulqSNJmJckPZyq3qUeSBsbgl6SBMfglaWAMfkkaGINfkgbG4Jekgekt+JNsm+T8JBcnuSzJW1v5x5Jcm2RFexzQVx0kSWvr8zz+O4DDq+q2JNsA30rypTbt9VV1So/rliTNorfgr66j/9va6DbtYef/kjRhvV65m2Qr4ALgYcD7quq8JK8E/jbJXwPnAMdW1R0zLLsMWAaw55579llNbQGWHHvWRNa76vijJrJeaWP0+uNuVd1dVQcAuwMHJfnPwHHAfsBjgJ2BN8yy7IlVtbSqli5atFZXE5KkDbQgZ/VU1c+BrwFHVNXq6twBfBQ4aCHqIEnq9HlWz6IkO7bh7YAnAVckWdzKAjwLWNlXHSRJa+uzjX8xcHJr578P8NmqOjPJuUkWAQFWAP+jxzpIkqbp86yeS4BHzVB+eF/rlCSt32bRH7+ktXkmkzaUXTZI0sAY/JI0MAa/JA2MwS9JA2PwS9LAeFaPtBEmdWaNtDE84pekgTH4JWlgDH5JGhiDX5IGxuCXpIEx+CVpYAx+SRoYg1+SBsbgl6SBMfglaWAMfkkaGINfkgbG4Jekgekt+JNsm+T8JBcnuSzJW1v5XknOS3J1ks8kuW9fdZAkra3PI/47gMOr6pHAAcARSQ4G3ga8s6oeBtwMvLzHOkiSpukt+KtzWxvdpj0KOBw4pZWfDDyrrzpIktbWaxt/kq2SrABuBM4GfgD8vKruarNcB+w2y7LLkixPsnzNmjV9VlOSBqXX4K+qu6vqAGB34CBgvzkse2JVLa2qpYsWLeqtjpI0NAtyVk9V/Rz4GnAIsGOSqVs+7g5cvxB1kCR1+jyrZ1GSHdvwdsCTgMvpPgCe22Y7BvhCX3WQJK2tz5utLwZOTrIV3QfMZ6vqzCTfA/4lyd8AFwEn9VgHSdI0vQV/VV0CPGqG8mvo2vslSRPglbuSNDAGvyQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvyQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sAY/JI0ML0Ff5I9knwtyfeSXJbkf7XytyS5PsmK9nhqX3WQJK1t6x6f+y7gdVV1YZIdgAuSnN2mvbOq3tHjuiVJs+gt+KtqNbC6Dd+a5HJgt77WJ0kaz4K08SdZAjwKOK8VvTrJJUk+kmSnWZZZlmR5kuVr1qxZiGpK0iD0HvxJHgB8HnhtVd0CfAB4KHAA3TeCf5hpuao6saqWVtXSRYsW9V1NSRqMXoM/yTZ0of/JqjoVoKpuqKq7q+o3wIeAg/qsgyTp3vo8qyfAScDlVfWPI+WLR2Z7NrCyrzpIktbW51k9hwIvAi5NsqKVvRE4OskBQAGrgD/psQ6SpGn6PKvnW0BmmPTFvtYpSVo/r9yVpIEx+CVpYAx+SRoYg1+SBsbgl6SBMfglaWAMfkkaGINfkgbG4JekgRkr+JP8ft8VkSQtjHGP+N+f5Pwkf5rkQb3WSJLUq7GCv6oeB7wQ2IPuFoqfSvKkXmsmSerF2G38VXUV8FfAG4AnAO9JckWS5/RVOUnS/Bu3jf8RSd4JXA4cDjy9qn6vDb+zx/pJkubZuN0y/xPwYeCNVXX7VGFV/STJX/VSM0lSL8YN/qOA26vqboAk9wG2rapfVdUnequdJGnejdvG/1Vgu5Hx7VuZJGkzM27wb1tVt02NtOHt+6mSJKlP4wb/L5McODWS5NHA7euYX5K0iRq3jf+1wOeS/ITuProPAZ7fW60kSb0ZK/ir6rtJ9gP2bUVXVtWv17VMkj2AjwO7AgWcWFXvTrIz8BlgCbAKeF5V3bxh1ZckzdVcOml7DPAI4EDg6CQvXs/8dwGvq6r9gYOBVyXZHzgWOKeq9gHOaeOSpAUy1hF/kk8ADwVWAHe34qI7op9RVa0GVrfhW5NcDuwGPBP4wzbbycDX6a4GliQtgHHb+JcC+1dVbchKkiwBHgWcB+zaPhQAfkrXFCRJWiDjNvWspPtBd86SPAD4PPDaqrpldFr7IJnxwyTJsiTLkyxfs2bNhqxakjSDcY/4dwG+l+R84I6pwqp6xroWSrINXeh/sqpObcU3JFlcVauTLAZunGnZqjoROBFg6dKlG/RNQ5K0tnGD/y1zfeIkAU4CLq+qfxyZdAZwDHB8+/uFuT63JGnDjXs65zeS/C6wT1V9Ncn2wFbrWexQ4EXApUlWtLI30gX+Z5O8HPgh8LwNq7okaUOMe1bPK4BlwM50Z/fsBpwAPHG2ZarqW3QXe81k1uUkSf0at6nnVcBBdGflUFVXJfmd3mqljbLk2LMmst5Vxx81kfVKmptxz+q5o6runBpJsjWznI0jSdq0jRv830jyRmC7dq/dzwH/t79qSZL6Mm7wHwusAS4F/gT4It39dyVJm5lxz+r5DfCh9pAkbcbGPavnWmZo06+qvee9RpKkXs2lr54p2wL/he7UTknSZmasNv6q+tnI4/qqehfdDdglSZuZcZt6DhwZvQ/dN4Bxvy1IkjYh44b3P4wM30W7c9a810aS1Ltxz+o5rO+KSJIWxrhNPX++runTet+UJG3C5nJWz2PoulQGeDpwPnBVH5XS5mlSfQRJmptxg3934MCquhUgyVuAs6rqv/VVMUlSP8btsmFX4M6R8TvxXrmStFka94j/48D5SU5r488CTu6nSpKkPo17Vs/fJvkS8LhW9NKquqi/akmS+jJuUw/A9sAtVfVu4Loke/VUJ0lSj8YK/iRvBt4AHNeKtgH+ua9KSZL6M+4R/7OBZwC/BKiqnwA79FUpSVJ/xg3+O6uqaF0zJ7l/f1WSJPVp3OD/bJIPAjsmeQXwVdZzU5YkH0lyY5KVI2VvSXJ9khXt8dQNr7okaUOMe1bPO9q9dm8B9gX+uqrOXs9iHwPeS3cq6Kh3VtU75lpRSdL8WG/wJ9kK+GrrqG19YX+PqvpmkiUbXjVJUh/W29RTVXcDv0nyoHla56uTXNKagnaabaYky5IsT7J8zZo187RqSdK4bfy3AZcmOSnJe6YeG7C+DwAPBQ4AVnPvfv7vpapOrKqlVbV00aJFG7AqSdJMxu2y4dT22ChVdcPUcJIPAWdu7HNKkuZmncGfZM+q+lFVzUu/PEkWV9XqNvpsYOW65pckzb/1NfWcPjWQ5PNzeeIknwa+Deyb5LokLwf+PsmlSS4BDgP+bK4VliRtnPU19WRkeO+5PHFVHT1D8UlzeQ5J0vxb3xF/zTIsSdpMre+I/5FJbqE78t+uDdPGq6oe2GvtJEnzbp3BX1VbLVRFJEkLYy798UuStgAGvyQNjMEvSQMz7pW7kgTAkmPPmti6Vx1/1MTWvSXxiF+SBsbgl6SBMfglaWAMfkkaGINfkgbG4JekgTH4JWlgDH5JGhiDX5IGxuCXpIEx+CVpYAx+SRoYg1+SBqa34E/ykSQ3Jlk5UrZzkrOTXNX+7tTX+iVJM+vziP9jwBHTyo4FzqmqfYBz2rgkaQH1FvxV9U3g36cVPxM4uQ2fDDyrr/VLkma20G38u1bV6jb8U2DXBV6/JA3exO7AVVWVpGabnmQZsAxgzz33XLB6zadJ3qlIkmaz0Ef8NyRZDND+3jjbjFV1YlUtraqlixYtWrAKStKWbqGD/wzgmDZ8DPCFBV6/JA1en6dzfhr4NrBvkuuSvBw4HnhSkquAP2rjkqQF1Fsbf1UdPcukJ/a1TknS+nnlriQNzMTO6pGkuZrUmXKrjj9qIuvti0f8kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvyQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sBs8TdimdSNGyRpU+URvyQNjMEvSQMzkaaeJKuAW4G7gbuqaukk6iFJQzTJNv7DquqmCa5fkgbJph5JGphJBX8BX0lyQZJlM82QZFmS5UmWr1mzZoGrJ0lbrkkF/2Or6kDgSOBVSR4/fYaqOrGqllbV0kWLFi18DSVpCzWR4K+q69vfG4HTgIMmUQ9JGqIFD/4k90+yw9Qw8GRg5ULXQ5KGahJn9ewKnJZkav2fqqp/nUA9JGmQFjz4q+oa4JELvV5JUmeL76tHkjbWJPv8WnX8UfP+nJ7HL0kDY/BL0sAY/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvyQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvyQNzESCP8kRSa5McnWSYydRB0kaqgUP/iRbAe8DjgT2B45Osv9C10OShmoSR/wHAVdX1TVVdSfwL8AzJ1APSRqkrSewzt2AH4+MXwf8wfSZkiwDlrXR25JcuQB1m7RdgJsmXYkJcvuHvf3gPlhr+/O2jXq+352pcBLBP5aqOhE4cdL1WEhJllfV0knXY1Lc/mFvP7gPFmr7J9HUcz2wx8j47q1MkrQAJhH83wX2SbJXkvsCLwDOmEA9JGmQFrypp6ruSvJq4MvAVsBHquqyha7HJmpQTVszcPs19H2wINufqlqI9UiSNhFeuStJA2PwS9LAGPw9S7JVkouSnNnG90pyXuuu4jPtB26S3K+NX92mLxl5juNa+ZVJnjKZLdkwSXZMckqSK5JcnuSQJDsnOTvJVe3vTm3eJHlP29ZLkhw48jzHtPmvSnLM5LZo7pL8WZLLkqxM8ukk227J74MkH0lyY5KVI2Xz9poneXSSS9sy70mShd3C9ZtlH7y9/R9ckuS0JDuOTJvxtZ2te5vZ3j9jqyofPT6APwc+BZzZxj8LvKANnwC8sg3/KXBCG34B8Jk2vD9wMXA/YC/gB8BWk96uOWz/ycB/b8P3BXYE/h44tpUdC7ytDT8V+BIQ4GDgvFa+M3BN+7tTG95p0ts25vbvBlwLbDfy+r9kS34fAI8HDgRWjpTN22sOnN/mTVv2yElv85j74MnA1m34bSP7YMbXtj1+AOzd/ncuBvYfeR+t9f4Zu36T3kFb8oPuGoVzgMOBM9sb9aaRF/8Q4Mtt+MvAIW146zZfgOOA40ae8575NvUH8KAWeplWfiWwuA0vBq5swx8Ejp4+H3A08MGR8nvNtyk/+O2V6ju31/VM4Clb+vsAWDIt9OblNW/Trhgpv9d8m9Jj+j6YNu3ZwCfb8Iyv7ej7YnS+deXIuA+bevr1LuAvgd+08QcDP6+qu9r4dXTBACNdWbTpv2jzz9TFxW5sHvYC1gAfbc1dH05yf2DXqlrd5vkpsGsbnm1bN9t9UFXXA+8AfgSspntdL2BY7wOYv9d8tzY8vXxz8zK6bysw932wrhwZi8HfkyRPA26sqgsmXZcJ2pru6+4HqupRwC/pvubfo7pDli32nOLWlv1Mug/B/wTcHzhiopWasC39NV+fJG8C7gI+Oak6GPz9ORR4RpJVdD2QHg68G9gxydSFc6PdVdzTlUWb/iDgZ2zeXVxcB1xXVee18VPoPghuSLIYoP29sU2fbVs3533wR8C1VbWmqn4NnEr33hjS+wDm7zW/vg1PL98sJHkJ8DTghe0DEOa+D37G7O+fsRj8Pamq46pq96paQvcj3blV9ULga8Bz22zHAF9ow2e0cdr0c9sb4wzgBe1sj72Afeh+3NrkVdVPgR8n2bcVPRH4Hvfe1un74MXtTI+DgV+05oEvA09OslM7gn5yK9sc/Ag4OMn27eyTqX0wmPdBMy+veZt2S5KD2/588chzbdKSHEHX9PuMqvrVyKTZXtsZu7dp74fZ3j/jmfQPIEN4AH/Ib8/q2bu9qFcDnwPu18q3beNXt+l7jyz/Jrpf969kEzyDYT3bfgCwHLgEOJ3uDI0H0/3ofRXwVWDnNm/obtLzA+BSYOnI87ys7ZurgZdOervmuA/eClwBrAQ+QXf2xhb7PgA+Tfd7xq/pvvW9fD5fc2Bp25c/AN7LtJMHNoXHLPvgaro2+xXtccL6Xlu6s56+36a9aaR8xvfPuA+7bJCkgbGpR5IGxuCXpIEx+CVpYAx+SRoYg1+SBsbgl6SBMfi1Xklu28jlT0mydxte1brUvSTJV5I8pJUnyblJHtjGd03yqSTXJLkgybeTPHvjt2ad9VyVZJc2/PWpLpGTfCzJtUkuTvL9JB9PsvvIcl+d6mZ4pOyEJIeOLLsiyYVJDhmZ511JHt+Gt07yf1oXxCva4019bu9087H9Se6b5JsjV5VqE2Twq1dJHk7XffA1I8WHVdUj6C7semMreypwcVXd0q7IPB34ZlXtXVWPprtqcfRS/XWts4/QeX1VPRLYF7gIOHekD/RP0HWnPOpg4Dsjyx5A10/RB1sdHwwcXFXfbPP8DV1fPr/f5n0csM04FWsfmn3/L4+1/VV1J92FWs/vuT7aCAa/xtYC5u3pbihyaZLnt/L7JHl/uptMnJ3ki0mmLid/IbNfTv5N4GEzzHc4cGdVnTA1Y1X9sKr+aR11e0mSM5KcSxc8JHl9ku+2bxdvHZn39PYt4rIky+ayD6rzTroeJo9sxWfQdQ889fy/B3y/qu5ex/b+MfCvbf7tgVcA/7Oq/qOt59aqess6tndJuht0fJzuKtY9knwgyfK2XaPbuyrJW9s3jkuT7NfKH9y+dV2W5MN0V9Fu9PbTfWi/cH3Ppckx+DUXz6HrguGRdJ2PvT1dh1vPoet7fH/gRXT9g085lK4b4pk8je4y/enzPRy4cAPqdyDw3Kp6QpIn0/V5clCr86OnmlWAl7VvEUuB17Sj77m6ENgPoKpuBu438jxH0kJ9mqcz8/Y+DPhRVd06xzrsA7y/qh5eVT+ku6R/KfAI4AlJHjEy701VdSDwAeAvWtmbgW9V1cOB04A957DudW3/SuAxc9wWLSCDX3PxWODTVXV3Vd0AfIPuH/yxwOeq6jfVdcz2tZFlFtP1yT/qa0lWAA8E/q6V7Txb8CV5X2tf/u566nd2Vf17G35ye1zEb0NqnzbtNUkupmuK2WOkfC6mHx3fSNdUA92NVkaD/+1te5fR9dkCM++X7omTl7Y2/h8n2WOmeZofVtV3Rsafl+RCum1+ON0H8ZRT298L6D6kobtL1D8DVNVZwM3rWNda1Zw2fs/2t286dybZYQ7PpwXkDzDq2+10HY+NOqyqbppWdleS+1TVb4DL6JpCAKiqV7UfHZevZ12/HBkO8HdV9cHRGZL8Id23lUOq6ldJvj5D/cbxKFqTUrMtcHtrttmxqn4yMu31VXXKtOVH98vVwJ5JdmhNPB+lu3nNSrrb783mnu1N16vjXwCPqaqbk3xs2nbd0f7ezfz838+4/SPj9wP+Yx7Wox54xK+5+H/A89PdQH4R3RHj+cC/AX/c2vp3peuNdMrl/LZde12upOtxEOBcYNskrxyZvv0c6/pl4GVJHgCQZLckv0PXv/3NLfT3o/sRdmztd47X0B2xT7XRB3gIsAo4jHt/45nNPfului56TwLem2Tb9pxb0d1ndVwPpPsg+EV7DY5cz/zQ/ebwX9v6jqTrOXWdxtj+qR+ub6ru/gPaBBn8movT6LpXvpgunP+yNe18nq7r2e/RNR1cSHfLQICzuPcHwWzuma+6LmOfRddOfW2S8+lu2v6GcStaVV+hu8n9t5NcSncTmB3owmrrJJcDx/PbM2/W5+2teej7dM1bh7UzWAAeDXynulvhzda+P930/fImum58Vya5iO5D9mTgJ2svuraqupiuiecKuu3+tzEWeyvw+CSX0f1O86N1zDvu9kP34XfWOPXWZNgts+ZFkgdU1W3taO984NCq+mmS7eiOgA+d4SyX0eUXAx+vqictUJXXqTUBvaSqVo0x77vpbpBxTmtj/4NxjnaTfAt4WlX9fGPrO982dPvb+KnAsVX1/V4rqQ1mG7/my5lJdqRrnvjf7ZsAVXV7kjfT3Qx61iPKqlqd5ENJHlhVtyxMlefNyqnQa2fOjOt1dGfSbHLBP0crR0L/vsDphv6mzSN+bVaSPAV427Tia6tqXq/qTXdv1NMnfTTevkGdM8OkJ1bVz3pc70vYBLZf/TD4JWlg/HFXkgbG4JekgTH4JWlgDH5JGpj/D/8rTzH3LISAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot\n",
    "matplotlib.pyplot.hist(lp_real - lprobs_rand)\n",
    "matplotlib.pyplot.xlabel(\"log(P(G_real|D)/P(G_rand|D))\")\n",
    "matplotlib.pyplot.ylabel(\"Frequency\")\n",
    "matplotlib.pyplot.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}