{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats import bernoulli\n",
    "import matplotlib.pyplot as plt; plt.rcdefaults()\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sympy import latex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def generate_data(n):\n",
    "    \"\"\"\n",
    "    Generate n datapoints for 3 variables A, H, B\n",
    "    with the model\n",
    "    A -> H -> B\n",
    "    \"\"\"\n",
    "\n",
    "    data = []\n",
    "    for i in range(n):\n",
    "        N_A = bernoulli.rvs(1./20)\n",
    "        N_H = bernoulli.rvs(1./3)\n",
    "        N_B = bernoulli.rvs(1./2)\n",
    "        A = N_A\n",
    "        H = (A + N_H) % 2\n",
    "        B = (H + N_B) % 2\n",
    "        data.append((A, H, B))\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEZCAYAAABvpam5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu0nXWd3/H3JxjpjEGPWRlDEjBJp54zGY9W0DmhKnQj\nBBAMQmdAmQjpTJFpO9OZUlcq06ySnGKHi21pp8EuRRtjBWfiKJcDFROMGxlXzcFRlDAkEZeJuRFE\nQMM4MwnJt3/s3w47O/vyPDv7cvY5n9dae+W5/J7f7/s8O+d8z3P7/RQRmJmZ5TGt1wGYmVn/cfIw\nM7PcnDzMzCw3Jw8zM8vNycPMzHJz8jAzs9ycPMxykPRZSTe1uO0ySV9td0xmveDkYZOCpB2SfiHp\nQMXnTzvQVKRP/g0j7oqIC9scj1lPvKrXAZi1SQDvi4hNjQpJOikiDp9gWzrB7c36ns88bFKT9M8l\nfVPSf5P0HLBK0kxJY5J+Jmlc0sckPVqxza9J2ijpp5K2SrqiqtpZkjZI+rmkoqQ3Vmx7RNLvSdou\n6QVJa6pieTRNS9LtkvanOL4v6c1p3cWSnkz175b0kYo63ifp8VT3NyW9pVPHzqwRJw+bTOqdEYwA\nPwTeAPwJ8AngADAbWA5cQ7oUJek1wEbg88CvAB8EPiFpUUUby4D/BMwCHgfuqmrvEuAdwFuBKyXV\nulR1AXA28KaIeB1wBfDTtO4zwHUR8VrgzcCmFNsZad2HgZnAJ4H7Jb26yXExazsnD5ssBNyb/iIv\nf65N6/ZGxB0RcQQ4BPwzYFVE/F1EPAWs45XE8z7gRxGxLiKORMTjwJcp/XIveyAi/jIiDgIrgX8i\naV7F+lsi4ucRsQv4OvC2GvEeAk4BFkmaFhHbIuKZtO4g8GZJr42In0XEd9Py64BPRsRjUfI54O+B\ns1o9aGatcvKwySKA90fE6ys+n07rdlWU+xVK9/oql+2umJ4PLK5MQsBvUzpLKbdztHxE/A3wPDC3\noo5nKqZ/AbzmuGBL92bWAHcA+yV9UtIpafVvAhcDO9JlsXJymA98pCq204A5DY6LWUc4edhUUPl0\n1E+Al4HTK5ZVTv8YeKQqCZ0SEb9fq7ykGZQuIe3NHVTE/4yIdwC/DgwCK9Lyb0fEZZQS3b3A+orY\n/nNVbDMi4s/ztm12opw8bDJp+hRUetLqy8BqSb8k6deAq3klwTwIDEr6kKTp6fMbqVy5jYslvSvd\na7gJ+H8RsadBTMfFJekdkhZLmk7p7OTvgMOpvWWSXpdiPQCUnw67E/iXkkbSDffXSLokJTCzrnLy\nsMlkrOo9jy9T+72MPwBeR+ny0jrgC5TuMxARByjdzP4gsAfYB9wMlG9KB6Ub5Kso3eA+A/hQRd3V\nbVW2Xzn9WuBTlC557QCeAz6e1n0I+JGkn1G6z7EsxfZXlG6Wr0nb/YDSzX6zrpMHg7KpTtKtwBsi\n4nd6HYtZv/CZh005koYkvTVd+hkBfhe4p9dxmfUTv2FuU9EplC5VzQX2A/8lIu7vbUhm/cWXrczM\nLDdftjIzs9z66rKVJJ8mmZnlFBFt78yz7848IqIvP6tWrep5DI6/93E4/v789HP8ndJ3ycPMzHrP\nycPMzHJz8uiSQqHQ6xBOiOPvLcffW/0efyf01aO6kqKf4jUz6zVJhG+Ym5nZRODkYWZmuTl5mJlZ\nbk4eZmaWm5OHmZnl5uRhZma5OXmYmVluTh5mZpZbX/WqC7By5ada3nZgAFasuK6N0ZiZTU19lzzm\nz2/9l//Ona0nHjMze4UvW5mZWW5OHmZmlpuTh5mZ5ebkYWZmuTVNHpJOlvSIpGmSHpL0gqSxBuWv\nkPSkpMOSzmxQrmZdktZLWphvN8zMrJuynHksAx6IiCPAbcDVTco/AVwOfKNJuXp13QlcnyGu3MbH\nhztRrZnZlJMleVwF3AcQEZuAlxoVjoitEbG9WaUN6ioCF2eIK7dt2xZ0olozsymnYfKQdBIwnCUZ\ntEtEHAL2SFrUrTbNzCyfZi8JzgIOdCOQKnuBBcBT1SvGxlYfnR4cLDA0VOhWTGZmE16xWKRYLHa8\nnSxvmFePfdvOQcTr1SXgSK0VS5eubmPzZmaTS6FQoFAoHJ0fHR3tSDvN7nk8B8yoWnbcQOqSbpZ0\nWY3tVVFmnqSHm9WVzAF2NonNzMx6pGHyiIjDwBZJQwCSHgXWA+dJ2iVpSSo6DOxLZS6XtAs4C3hQ\n0ldSmTnAy+W669UlaTpwWkRsbddOlg0N7Wh3lWZmU1KWy1Z3AZcBt0bE2XXKTI+IzQARcQ9wT40y\ni4E15ZkGdZ0LPJAhrtxGRrYA7+xE1WZmU0qWR3XvBi6RVO8SExFxUbNKIuKOiMiSFK4Fbs9QzszM\neqTpmUdEHATO6UIs5fau7FZbZmbWGvdtZWZmuTl5mJlZbk4eZmaWW98NQ3siQ8kODLQxEDOzKUwR\n7XxhvLMkRT/Fa2bWa5KIiLpPy7bKl63MzCw3Jw8zM8vNycPMzHLruxvmK1ee2A3zFSuua2M0ZmZT\nU98lj/nzW//lfyJPapmZ2St82crMzHJz8jAzs9ycPMzMLDcnDzMzy61p8pB0sqRHJE2T9JCkFySN\nNSh/haQnJR2WdGaDcsslbU+fayqWr5e0MP+umJlZt2Q581gGPBARR4DbgKublH8CuBz4Rr0CkmYC\nNwIj6bNKUrnnqTuB6zPEldv4+HAnqjUzm3KyJI+rgPsAImIT8FKjwhGxNSK2N6nzQmBDRLwYES8C\nG4HyaIRF4OIMceW2bduCTlRrZjblNEwekk4ChjMkg7zmArsr5ncD8wAi4hCwR9KiNrdpZmZt0uwl\nwVnAgW4EUmUvsAB4qnrF2Njqo9ODgwWGhgrdisnMbMIrFosUi8WOt5PlDfPqrnzb0Sf6HqBQMX86\nsKmqzSO1Nly6dHUbmjczm5wKhQKFQuHo/OjoaEfaaXbP4zlgRtWy4/qFl3SzpMtqbK+KMvMkPZxm\nNwAXSBqQ9HpgCfDViu3mADubBW9mZr3RMHlExGFgi6QhAEmPAuuB8yTtkrQkFR0G9qUyl0vaBZwF\nPCjpK6nMHODlVO/zwE3AY8A4MJpunCNpOnBaRGxt326WDA3taHeVZmZTUpbLVncBlwG3RsTZdcpM\nj4jNABFxD3BPjTKLgTXlmYhYC6ytUe5c4IEMceU2MrIFeGcnqjYzm1KyPKp7N3CJpLrDGEbERfXW\nVZS5IyKyJIVrgdszlDMzsx5peuYREQeBc7oQS7m9K7vVlpmZtcZ9W5mZWW5OHmZmlpuTh5mZ5dZ3\nw9CeyFCyAwPNy5iZWXOKaMcL490hKfopXjOzXpNERNR9WrZVvmxlZma5OXmYmVluTh5mZpZb390w\nX7my9g3zgQFYseK6LkdjZjY19V3ymD+/doI4kaewzMwsH1+2MjOz3Jw8zMwsNycPMzPLzcnDzMxy\na5o8JJ0s6RFJ0yQ9JOkFSWMNys+UtFHSdkkbJNXsFKReXZLWS1qYf1fMzKxbspx5LAMeiIgjwG3A\n1U3K3wBsjIhB4GtpvpZ6dd0JXJ8hrmPs3DmHYjHvVmZm1oosyeMq4D6AiNgEvNSk/KXAujS9jtIQ\ntsdpUFcRuDhDXMfYuXOuk4eZWZc0TB6STgKGI2J7jjpnR8T+NL0fmJ0noIg4BOyRtCjPdmZm1j3N\nXhKcBRxotfKICEmtdIO7F1gAPFW9Ymxs9dHpwcECQ0MFwJetzMwAisUixS78Mszyhnl1V77NksF+\nSadGxDOS5gDPNihbry4BR2qtWLp0dc0N5s/fR6Ewt0loZmaTW6FQoFAoHJ0fHR3tSDvN7nk8B8yo\nWnZcv/CSbpZUvrdxP7A8TS8H7k1l5kl6uFldyRxgZ5PYzMysRxomj4g4DGyRNAQg6VFgPXCepF2S\nlqSiw8C+NH0LsETSduA9aR5KCeHlct316pI0HTgtIrbm2ZH58/dSkWzNzKyDsly2uovSE1O3RsTZ\ndcpMj4jNABHxPHB+jTKLgTXlmQZ1nQs8kCGuY5QuW+XdyszMWpHlUd27gUsk1R3GMCIualZJRNwR\nEVmSwrXA7RnKmZlZjzQ984iIg8A5XYil3N6V3WrLzMxa476tzMwsNycPMzPLzcnDzMxy67thaOsN\nNztQs+9eMzPrBEW00ntIb0iKforXzKzXJBERdZ+WbZUvW5mZWW5OHmZmlpuTh5mZ5dZ3N8xXrqx/\nw3zFiuu6HI2Z2dTUd8lj/vzaCaLeU1hmZtZ+vmxlZma5OXmYmVluTh5mZpZb0+Qh6WRJj0iaJukh\nSS9IGmtQfqakjZK2S9ogqea735KWpzLbJV1TsXy9pIWt7Y6ZmXVDljOPZcADEXEEuA24ukn5G4CN\nETEIfC3NH0PSTOBGYCR9VlUkmTuB67OFb2ZmvZAleVwF3AcQEZuAl5qUvxRYl6bXURqFsNqFwIaI\neDEiXgQ2AuUBpYrAxRniOsb4+HDeTczMrEUNk4ekk4DhiNieo87ZEbE/Te8HZtcoMxfYXTG/G5gH\nEBGHgD2SFuVok23bFuQpbmZmJ6DZmccs4ECrladeDFvpyXAvsKDVds3MrLOyvCRY3Rtjs2SwX9Kp\nEfGMpDnAszXK7AEKFfOnA5uq2jxSq/KxsdVHpwcHCwwNFWoVMzObkorFIsVisePtNOySPV222h0R\ncyqWFYCPRMTSimU3A5sj4l5JtwE/jYhbJd0ADETEDZLmAesi4vx0w/zbwJmUEsVfAWem+x9IegT4\nvYjYWhVPfPKTteP92Mf28uMfz23hEJiZTV496ZI9Ig4DWyQNpSAeBdYD50naJWlJKjoM7EvTtwBL\nJG0H3pPmAeYAL6d6nwduAh4DxoHRisQxHTitOnGYmdnEkeWy1V2Unpi6NSLOrlNmekRshqOJ4fwa\nZRYDa8ozEbEWWFuj3LnAAxniOsbQ0A5K9+HNzKzTsjyqezdwiaS6pz0RcVG9dRVl7oiILEnhWuD2\nDOWOMTKyJe8mZmbWoqZnHhFxEDinC7GU27uyW22ZmVlr3LeVmZnl5uRhZma5OXmYmVluTh5mZpZb\n3w1DW2+42YGaHb+bmVknNHzDfKKRFP0Ur5lZr/XkDXMzM7NanDzMzCw3Jw8zM8ut726Yr1xZ+4a5\nWTsMDMCKFdf1OgyzCa/vksf8+f7Bts6p9zSfmR3Ll63MzCw3Jw8zM8vNycPMzHJrmjwknSzpEUnT\nJC2XtD19rqlT/gpJT0o6LOnMBvU+JOkFSWNVy9dLWph/V8zMrFuynHksozSy3wBwIzCSPqsk1eoU\n5AngcuAbTeq9Dbi6xvI7geszxGVmZj2SJXlcBdwHXAhsiIgX03jjG4HjRhCMiK0Rsb1ZpRGxCXip\nxqoicHGGuMzabnx8uNchmPWFhslD0knAcEoG84DdFat3p2VtFRGHgD2SFrW7brNmtm1b0OsQzPpC\nszOPWcCBNN3NHgn3Agu62J6ZmeWQ5SXBcm+Me4BCxfLTgU0n2H69hCTgSK0VY2Orj04PDhYYGirU\nKmZmNiUVi0WKxWLH22mWPJ4DZqTpDcCfpJvkApYAHwWQdDOwOSLurdr+aDfAkuYB6yLi/Frrq8wB\ndtZasXTp6iYhm5lNXYVCgUKhcHR+dHS0I+00vGwVEYeBLZKGIuJ54CbgMWAcGE03zgGGgX0Aki6X\ntAs4C3hQ0ldSmTnAy+W6JT0KrAfOk7RL0pK0fDpwWkRsbddOmplZe2W5bHUXcBlwa0SsBdbWKDM9\nIjYDRMQ9wD01yiwG1pRnIuLsOu2dS+nRYLOuGxraAcztcRRmE1+WR3XvBi6RVHckqog47pHdGmXu\niIgsSeFa4PYM5czabmRkS69DMOsLTc88IuIgcE4XYim3d2W32jIzs9a4byszM8vNycPMzHJz8jAz\ns9ycPMzMLLe+G4bWw4RaJw3U6ifazI6jiG52WXViJEU/xWtm1muSiIi6r1q0ypetzMwsNycPMzPL\nre/ueaxc6Xse1l4DA7BixXW9DsOsr/Rd8pg/3z/k1l5+CMMsP1+2MjOz3Jw8zMwsNycPMzPLzcnD\nzMxya5o8JJ0s6RFJ0yQtl7Q9fa6pU/4KSU9KOizpzAb11qxL0npJC1vbHTMz64YsZx7LKI3sNwDc\nCIykz6o0nnm1J4DLgW/Uq1DSzAZ13Qlcn3UHzE7U+Phwr0Mw6ztZksdVwH3AhcCGiHgxjV2+EThu\nBMGI2BoR25vU2aiuInBxxvjNTti2bQt6HYJZ32mYPCSdBAynZDAP2F2xenda1oq59eqKiEPAHkmL\nWqzbzMw6rNlLgrOAA2m6mz0S7gUWAE9VrxgbW310enCwwNBQoVsxmZlNeMVikWKx2PF2srxhXu6N\ncQ9QqFh+OrCpxXab1SXgSK0Nly5d3WKTZmaTX6FQoFAoHJ0fHR3tSDvN7nk8B8xI0xuACyQNSHo9\nsAT4KoCkmyVdVmP7o90AS5on6eFmdSVzgJ2598bMzLqiYfKIiMPAFklDEfE8cBPwGDAOjKab3QDD\nwD4ASZdL2gWcBTwo6SupzBzg5VRv3bokTQdOi4it7dtNs/qGhnb0OgSzvpPlstVdwGXArRGxFlhb\no8z0iNgMEBH3APfUKLMYWFOeaVDXuZQeDTbripGRLcA7ex2GWV/J8qju3cAlkuqORBURxz2yW6PM\nHRGRJSlcC9yeoZyZmfVI0zOPiDgInNOFWMrtXdmttszMrDXu28rMzHJz8jAzs9ycPMzMLLe+G4bW\nQ4Zauw3U6t7TzBpSRDd7HTkxkqKf4jUz6zVJRETdp2Vb5ctWZmaWm5OHmZnl5uRhZma59d0N85Ur\nfcPczPrLwACsWHFdr8Noq75LHvPnT64vwMwmv8n4lKgvW5mZWW5OHmZmlpuTh5mZ5ebkYWZmuTVN\nHpJOlvSIpGmSlkvanj7X1Ck/U9LGVGaDpJqdP0h6SNILksaqlq+XtLC13TEzs27IcuaxjNLIfgPA\njcBI+qyqkxhuADZGxCDwtTRfy23A1TWW3wlcnyEuM7O+MD4+3OsQ2i5L8rgKuA+4ENgQES+m8cY3\nArVGELwUWJem11EawvY4EbEJeKnGqiJwcYa4zMz6wrZtC3odQts1TB6STgKGI2I7MA/YXbF6d1pW\nbXZE7E/T+4HZeQKKiEPAHkmL8mxnZmbd0+wlwVnAgTSduzvbiAhJrXSDuxdYADxVvWJsbPXR6cHB\nAkNDhRaqNzObnIrFIsVisePtZHnDvNyV7x6gULH8dGBTjfL7JZ0aEc9ImgM826DueolFwJFaK5Yu\nXd0wWDOzqaxQKFAoFI7Oj46OdqSdZvc8ngNmpOkNwAWSBiS9HlgCfBVA0s2Syvc27geWp+nlwL2p\nzDxJD1fVX6+P+TnAzsx7YWZmXdUweUTEYWCLpKGIeB64CXgMGAdG041zgGFgX5q+BVgiaTvwnjQP\npYTwcrluSY8C64HzJO2StCQtnw6cFhFb27GDZma9NjS0o9chtF2Wy1Z3UXpi6taIWAusrVFmekRs\nBkhJ5vwaZRYDa8ozEXF2nfbOpfRosJnZpDAysgV4Z6/DaKssj+reDVwiqe4whhFR65Hd6jJ3RESW\npHAtcHuGcmZm1iNNzzwi4iBwThdiKbd3ZbfaMjOz1rhvKzMzy83Jw8zMcnPyMDOz3PpuGNrJOJyj\nmU1uAzX7Fu9vimil95DekBT9FK+ZWa9JIiLqPi3bKl+2MjOz3Jw8zMwsNycPMzPLre9umK9c6Rvm\nZtY/BgZgxYrreh1G2/Vd8pg/f/J9CWY2eU3WJ0R92crMzHJz8jAzs9ycPMzMLDcnDzMzy61p8pB0\nsqRHJE2TtFzS9vS5pk75mZI2pjIbJNV8Mb9eXZLWS1rY+i6ZmVmnZTnzWEZpZL8B4EZgJH1W1UkM\nNwAbI2IQ+FqaP4akmQ3quhO4Pud+mJlNSOPjw70OoSOyJI+rgPuAC4ENEfFiGrt8I1BrBMFLgXVp\neh2lIWyrNaqrCFyceQ/MzCawbdsW9DqEjmiYPCSdBAxHxHZgHrC7YvXutKza7IjYn6b3A7NrlJlb\nr66IOATskbQo0x6YmVnXNXtJcBZwIE3n7s42IkJSK93g7gUWAE9VrxgbW310enCwwNBQoYXqzcwm\np2KxSLFY7Hg7Wd4wL3fluwcoVCw/HdhUo/x+SadGxDOS5gDP1ijTrC4BR2oFs3Tp6gwhm5lNTYVC\ngUKhcHR+dHS0I+00u+fxHDAjTW8ALpA0IOn1wBLgqwCSbpZUvrdxP7A8TS8H7k1l5kl6uFldyRxg\nZ+u7ZWZmndQweUTEYWCLpKGIeB64CXgMGAdG081ugGFgX5q+BVgiaTvwnjQPpYTwcqq3bl2SpgOn\nRcTW9uyimVnvDA3t6HUIHZHlstVdlJ6YujUi1gJra5SZHhGb4WhiOL9GmcXAmvJMg7rOpfRosJlZ\n3xsZ2QK8s9dhtF2WR3XvBi6RVHcYw4io9chudZk7IiJLUrgWuD1DOTMz65GmZx4RcRA4pwuxlNu7\nslttmZlZa9y3lZmZ5ebkYWZmuTl5mJlZbn03DO1kHdLRzCangZr9ivc/RbTSe0hvSIp+itfMrNck\nERF1n5ZtlS9bmZlZbk4eZmaWm5OHmZnl5uRhZma5OXmYmVluTh5mZpabk4eZmeXm5GFmZrk5eZiZ\nWW5OHl3SjQHpO8nx95bj761+j78TnDy6pN//8zn+3nL8vdXv8XeCk4eZmeXm5GFmZrn1Xa+6vY7B\nzKzfdKJX3b5KHmZmNjH4spWZmeXm5GFmZrn1RfKQdJGkrZJ+IOmjvY6nkqQdkr4v6buSxtOymZI2\nStouaYOkgYryf5z2Y6ukCyqWv13SE2nd/+hQrP9b0n5JT1Qsa1uskk6W9Odp+bckze9C/Ksl7U7H\n/7uS3juB4z9d0tclPSlpi6Q/TMv74jtoEP+E/w4k/QNJmyU9LumvJd2clvfLsa8Xf++OfURM6A9w\nEvA0sACYDjwOLOp1XBXx/QiYWbXsNuDfp+mPArek6V9P8U9P+/M0r9x3GgdG0vT/BS7qQKxnA2cA\nT3QiVuBfA59I0x8A/qwL8a8C/l2NshMx/lOBt6XpGcA2YFG/fAcN4u+L7wD45fTvq4BvAe/ul2Pf\nIP6eHft+OPMYAZ6OiB0RcQj4M+D9PY6pWvWTDJcC69L0OuCyNP1+4AsRcSgidlD6QhdLmgOcEhHj\nqdznKrZpm4h4FHihg7FW1vUl4LwuxA/HH3+YmPE/ExGPp+mXgKeAefTJd9AgfuiD7yAifpEmX03p\nj9IX6JNj3yB+6NGx74fkMQ/YVTG/m1f+w04EATws6duSPpyWzY6I/Wl6PzA7Tc+lFH9ZeV+ql++h\ne/vYzliPflcR8TLwM0kzOxR3pX8j6XuSPlNx2WFCxy9pAaWzqM304XdQEf+30qIJ/x1ImibpcUrH\n+OsR8SR9dOzrxA89Ovb9kDwm+rPE74qIM4D3Ar8v6ezKlVE6B5zo+wD0V6wV/hewEHgbsA/4r70N\npzlJMyj9ZfdHEXGgcl0/fAcp/r+gFP9L9Ml3EBFHIuJtwGnAOZLOrVo/oY99jfgL9PDY90Py2AOc\nXjF/Osdmzp6KiH3p358A91C6zLZf0qkA6TTx2VS8el9Oo7Qve9J05fI9nY38qHbEurtimzemul4F\nvC4inu9c6BARz0YCfJrS8S/HMuHilzSdUuL4PxFxb1rcN99BRfyfL8ffb99BRPwMeBB4O3107GvE\n/45eHvt+SB7fBt4kaYGkV1O6kXN/j2MCQNIvSzolTb8GuAB4glJ8y1Ox5UD5l8T9wAclvVrSQuBN\nwHhEPAP8XNJiSQKurtim09oR63016vot4GudDj79wJddTun4T8j4U3ufAf46Iv57xaq++A7qxd8P\n34GkWeVLOpJ+CVgCfJf+OfY14y8nvqS7x77ZHf6J8KF0SWgbpZs+f9zreCriWkjpiYbHgS3l2ICZ\nwMPAdmADMFCxzX9I+7EVuLBi+dvTF/808KcdivcLwF7gIKVrm7/TzliBk4H1wA8oXQtf0OH4f5fS\nDb/vA9+j9IM/ewLH/27gSPr/8t30uahfvoM68b+3H74D4C3Ad1Ls3wdWtPtntcPHvl78PTv27p7E\nzMxy64fLVmZmNsE4eZiZWW5OHmZmlpuTh5mZ5ebkYWZmuTl5mJlZbk4eNuFIukzSEUlDHaj7s5J+\ns8byBaro6r1q3cdV6oL81hNs++HyS6VpPtd+qtT9/8yK+YKksSbbLJD0typ11/24pG9KGkzr/rGk\nz7S6Pza1OXnYRHQV8ED6t91aebHpw8BbIiLTWDKpa4fqZe8BtsWxfVnl3c/q2LPuy9MRcUaU+kVa\nR+nlMSLie8CvSnpDxnrMjnLysAkldbq3GPgDSl3RlJcXJBUlfVHSU5I+n5a/Q68MhPOEpCNp+Ycl\njae/tv8idelQdk76C/yHtc5CquK5n9LYFd+RdKWktZXbSHqpIr5HJd0HPFmjqt/mlW4g6u5nlkNU\nZzqr1wGV/RV9BbiihXpsinPysInm/cBDEfFj4CeSzqxY9zbgjygNdPMPJb0rIr6d/qo+g9Ivwo+n\nsl+KiJH01/ZTwL9IywWcGhHvAt4H3NIomIi4FPjb1Mb6WkUqps8A/jAial2Gehelftqy7Gc9Ar5e\nTpbAnWQ7+/jVtM3TwL8Fbq9YNw6ck6EOs2M4edhEcxXwxTT9RY69pDMeEXuj1KfO45RGSANA0geA\nM4Eb0qK3pDOB7wPLKCUcKP2yLfcG+xSvjN/QDuMRsbPOurlxbA+ljfazngAKFcnyWrKdffwwbfOP\ngOuBT1Ws20fFcTTL6rhrs2a9km4GnwsMSwpKo6UFsCIV+fuK4odJ/38lDVMajvPseKWzts8Cl0bE\nE5KWA4WKbQ9WNpszzJdJf3RJmkZpVLeyv8lSQYb9zKqVy1ZjwNqqOtzBneXmMw+bSH4L+FxELIiI\nhRHxRuBHqhpgq1LqpvoLwNUR8dOKVTOAZ1Qaf+JDtO8X5A5KvZJCadjO6Rm321vxpFTD/ZS0NW9Q\nkkYkrWtekndT6k21bA5Q72zJrC4nD5tIPkhpQK1KX6J0SafWKG9B6Rf4G4FPp+v630nr/iOlIV7/\nktI9j+piAzOOAAAArUlEQVTtmk3XK38n8E9VGg70LOClDNuT4viNNF1vPz8oaVaDOmrtf3nZG4Ff\nUFv5nsfjwMcoXe4qGwG+0aBNs5rcJbtZF6g0ZOgHIuJfNSl3CbAwItbkrP82SmczW3JuVwSujIhn\nm5U1q+TkYdYlkh4GLq9616NnJL2V0tNh1zYtbFbFycPMzHLzPQ8zM8vNycPMzHJz8jAzs9ycPMzM\nLDcnDzMzy83Jw8zMcvv/FeE+zT68A9wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc41c826cd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p(H=0|A=0) = 0.67\n",
      "p(H=0|A=1) = 0.32\n",
      "p(H=0|B=0) = 0.65\n",
      "p(H=0|B=1) = 0.65\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "\n",
    "data = generate_data(100000)\n",
    "\n",
    "# Example data\n",
    "people = ('(0,0,0)', '(0,0,1)', '(0,1,0)', '(0,1,1)', '(1,0,0)', '(1,0,1)', '(1,1,0)', '(1,1,1)')\n",
    "get_bin = lambda x, n: x >= 0 and str(bin(x))[2:].zfill(n) or \"-\" + str(bin(x))[3:].zfill(n)\n",
    "y_pos = np.arange(len(people))\n",
    "# ('(0,0,0)', '(0,0,1)', '(0,1,0)', '(0,1,1)', '(1,0,0)', '(1,0,1)', '(1,1,0)', '(1,1,1)')\n",
    "count = [sum([1 for d in data if map(int, list(get_bin(i, 3))) == list(d)]) for i in range(8)]\n",
    "\n",
    "error = np.random.rand(len(people))\n",
    "\n",
    "plt.barh(y_pos, count, xerr=error, align='center', alpha=0.4)\n",
    "plt.yticks(y_pos, people)\n",
    "plt.xlabel('Anzahl fur (A, H, B)')\n",
    "plt.title('Ergebnisse')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "from sympy.interactive import printing\n",
    "printing.init_printing(use_latex=True)\n",
    "p_1 = (count[0] + count[1])  # h=0 and A=0\n",
    "p_2 = (count[4] + count[5])  # h=0 and A=1\n",
    "p_3 = (count[2] + count[3])  # h=1 and A=0\n",
    "p_4 = (count[6] + count[7])  # h=1 and A=1\n",
    "print(latex(\"p(H=0|A=0) = %0.2f\" % (p_1/float(p_1+p_3))))\n",
    "print(latex(\"p(H=0|A=1) = %0.2f\" % (p_2/float(p_2+p_4))))\n",
    "\n",
    "p_1 = (count[0] + count[4])  # h=0 and B=0\n",
    "p_2 = (count[1] + count[5])  # h=0 and B=1\n",
    "p_3 = (count[2] + count[6])  # h=1 and B=0\n",
    "p_4 = (count[3] + count[7])  # h=1 and B=1\n",
    "print(latex(\"p(H=0|B=0) = %0.2f\" % (p_1/float(p_1+p_3))))\n",
    "print(latex(\"p(H=0|B=1) = %0.2f\" % (p_2/float(p_2+p_4))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
