From af1eb5d57f2966851d9b72a56930c98609ced2e2 Mon Sep 17 00:00:00 2001 From: astrokathi Date: Tue, 11 Apr 2023 15:29:20 +0400 Subject: [PATCH] fix(Fixed-the-Sprinkler-notebook): added pos argument to the nx.draw method, inline comments added wherever necessary --- Sprinkler.ipynb | 375 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 374 insertions(+), 1 deletion(-) diff --git a/Sprinkler.ipynb b/Sprinkler.ipynb index 54425a7..7f6230b 100644 --- a/Sprinkler.ipynb +++ b/Sprinkler.ipynb @@ -1 +1,374 @@ -{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"Sprinkler.ipynb","provenance":[],"collapsed_sections":[]},"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.8.5"}},"cells":[{"cell_type":"markdown","metadata":{"id":"uYCyEcvDyvaG"},"source":["# Inference in Bayesian Networks\n","\n","This notebook illustrates the Bayesian Network for the \"Wet Grass\" example discussed in the course book \"Inference & Causality\", unit 1.2, based\n","on: Murphy, K. (2001). An introduction to graphical models.Rap. tech,96,1–19"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"8D7_xzjnyvaI","executionInfo":{"status":"ok","timestamp":1611834585233,"user_tz":-60,"elapsed":9618,"user":{"displayName":"Ulrich Kerzel","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64","userId":"10582137377195463323"}},"outputId":"09460824-4c8c-47f6-ec77-24ec6b5cd4ef"},"source":["!pip install pgmpy\n","\n","from pgmpy.models import BayesianModel\n","from pgmpy.factors.discrete import TabularCPD\n","import networkx as nx\n","import pylab as plt"],"execution_count":2,"outputs":[{"output_type":"stream","text":["Collecting pgmpy\n","\u001b[?25l Downloading https://files.pythonhosted.org/packages/06/19/d508949e8ac7b32e639f15e854a5f5ed710a4118e4f6692bddaccc390d88/pgmpy-0.1.13-py3-none-any.whl (324kB)\n","\r\u001b[K |█ | 10kB 18.1MB/s eta 0:00:01\r\u001b[K |██ | 20kB 21.4MB/s eta 0:00:01\r\u001b[K |███ | 30kB 25.5MB/s eta 0:00:01\r\u001b[K |████ | 40kB 18.6MB/s eta 0:00:01\r\u001b[K |█████ | 51kB 17.8MB/s eta 0:00:01\r\u001b[K |██████ | 61kB 16.3MB/s eta 0:00:01\r\u001b[K |███████ | 71kB 14.0MB/s eta 0:00:01\r\u001b[K |████████ | 81kB 13.6MB/s eta 0:00:01\r\u001b[K |█████████ | 92kB 13.5MB/s eta 0:00:01\r\u001b[K |██████████ | 102kB 12.5MB/s eta 0:00:01\r\u001b[K |███████████ | 112kB 12.5MB/s eta 0:00:01\r\u001b[K |████████████▏ | 122kB 12.5MB/s eta 0:00:01\r\u001b[K |█████████████▏ | 133kB 12.5MB/s eta 0:00:01\r\u001b[K |██████████████▏ | 143kB 12.5MB/s eta 0:00:01\r\u001b[K |███████████████▏ | 153kB 12.5MB/s eta 0:00:01\r\u001b[K |████████████████▏ | 163kB 12.5MB/s eta 0:00:01\r\u001b[K |█████████████████▏ | 174kB 12.5MB/s eta 0:00:01\r\u001b[K |██████████████████▏ | 184kB 12.5MB/s eta 0:00:01\r\u001b[K |███████████████████▏ | 194kB 12.5MB/s eta 0:00:01\r\u001b[K |████████████████████▏ | 204kB 12.5MB/s eta 0:00:01\r\u001b[K |█████████████████████▏ | 215kB 12.5MB/s eta 0:00:01\r\u001b[K |██████████████████████▏ | 225kB 12.5MB/s eta 0:00:01\r\u001b[K |███████████████████████▎ | 235kB 12.5MB/s eta 0:00:01\r\u001b[K |████████████████████████▎ | 245kB 12.5MB/s eta 0:00:01\r\u001b[K |█████████████████████████▎ | 256kB 12.5MB/s eta 0:00:01\r\u001b[K |██████████████████████████▎ | 266kB 12.5MB/s eta 0:00:01\r\u001b[K |███████████████████████████▎ | 276kB 12.5MB/s eta 0:00:01\r\u001b[K |████████████████████████████▎ | 286kB 12.5MB/s eta 0:00:01\r\u001b[K |█████████████████████████████▎ | 296kB 12.5MB/s eta 0:00:01\r\u001b[K |██████████████████████████████▎ | 307kB 12.5MB/s eta 0:00:01\r\u001b[K |███████████████████████████████▎| 317kB 12.5MB/s eta 0:00:01\r\u001b[K |████████████████████████████████| 327kB 12.5MB/s \n","\u001b[?25hRequirement already satisfied: networkx in /usr/local/lib/python3.6/dist-packages (from pgmpy) (2.5)\n","Requirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.19.5)\n","Requirement already satisfied: scipy in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.4.1)\n","Requirement already satisfied: statsmodels in /usr/local/lib/python3.6/dist-packages (from pgmpy) (0.10.2)\n","Requirement already satisfied: joblib in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.0.0)\n","Requirement already satisfied: torch in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.7.0+cu101)\n","Requirement already satisfied: tqdm in /usr/local/lib/python3.6/dist-packages (from pgmpy) (4.41.1)\n","Requirement already satisfied: pyparsing in /usr/local/lib/python3.6/dist-packages (from pgmpy) (2.4.7)\n","Requirement already satisfied: scikit-learn in /usr/local/lib/python3.6/dist-packages (from pgmpy) (0.22.2.post1)\n","Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (from pgmpy) (1.1.5)\n","Requirement already satisfied: decorator>=4.3.0 in /usr/local/lib/python3.6/dist-packages (from networkx->pgmpy) (4.4.2)\n","Requirement already satisfied: patsy>=0.4.0 in /usr/local/lib/python3.6/dist-packages (from statsmodels->pgmpy) (0.5.1)\n","Requirement already satisfied: dataclasses in /usr/local/lib/python3.6/dist-packages (from torch->pgmpy) (0.8)\n","Requirement already satisfied: future in /usr/local/lib/python3.6/dist-packages (from torch->pgmpy) (0.16.0)\n","Requirement already satisfied: typing-extensions in /usr/local/lib/python3.6/dist-packages (from torch->pgmpy) (3.7.4.3)\n","Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.6/dist-packages (from pandas->pgmpy) (2.8.1)\n","Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.6/dist-packages (from pandas->pgmpy) (2018.9)\n","Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from patsy>=0.4.0->statsmodels->pgmpy) (1.15.0)\n","Installing collected packages: pgmpy\n","Successfully installed pgmpy-0.1.13\n"],"name":"stdout"}]},{"cell_type":"markdown","metadata":{"id":"eS15uTMlzcKT"},"source":["First we define the Bayesian network by specifying which nodes are connected to which.\n","The following shortcuts are used:\n","* C: Cloudy\n","* S: Sprinkler\n","* R: Rain\n","* W: Wet grass\n","\n","We use the module ```networkx``` do visualize the network. The output is perhaps a bit counter-intuitive - but careful investigation shows that the arrows point indeed into the correct nodes:\n","\n","* *Cloudy* is a common cause to both *Sprinkler* and *Rain*\n","* *Sprinkler* is a cause to *Wet grass*\n","* *Rain* is a cause to *Wet grass*."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":319},"id":"8r5E1lCGyvaM","executionInfo":{"status":"ok","timestamp":1611834661134,"user_tz":-60,"elapsed":748,"user":{"displayName":"Ulrich Kerzel","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64","userId":"10582137377195463323"}},"outputId":"5a3fd0bf-b9cc-4351-b945-856343200b2a"},"source":["#define model\n","model = BayesianModel([('C','S'),('C','R'),('S','W'),('R','W')])\n","nx.draw(model, with_labels=True)\n","plt.show()\n"],"execution_count":3,"outputs":[{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVyU1eIG8GdYZIhFxd3AjUnJFHNHZCbtZhuiuKCmKGSmhWmbmUlZmZa/XMpc0xYrzSUyNTUrU2RXXHINFRQURAWVnQFm5v394WVuhOY2M2dm3uf7120YZh4/F+bhnPec8yokSZJAREQkEw6iAxAREVkSi4+IiGSFxUdERLLC4iMiIllh8RERkayw+IiISFZYfEREJCssPiIikhUWHxERyQqLj4iIZIXFR0REssLiIyIiWWHxERGRrLD4iIhIVlh8REQkKyw+IiKSFRYfERHJCouPiIhkhcVHRESywuIjIiJZYfEREZGssPiIiEhWnEQHILJn+SUViDmQjbSLRSjS6uCpdIJfU0+EdfVGA3cX0fGIZEkhSZIkOgSRvTl8vgBLYtOx51QeAKBCZzB+TenkAAlAn3aNEPWICp186glKSSRPLD4iE1udkonZ29Og1enxb79dCgWgdHJE9NN+CA9oZbF8RHLHqU4iE7peen+hvMpwy+dKElBepcfs7X8BAMuPyEK4uIXIRA6fL8Ds7Wm4GLcOlza8W+NrOZ8/f8PHSk/sQXmVAbO3p+FIdoEl4xLJFouPyESWxKZDq9PDxacDKnL+gmTQAwB0JVch6fWovJRR4zHdtVy4+HQAAGh1eiyNTReWnUhOWHxEJpBfUoE9p/IgSYBLswcAvR6Vl88CACrOH4eypT+cvbxrPOZUrxmcPBoAuD7tuftkHq6UVAj7NxDJBYuPyARiDmQb/7fC0Rl1mrdFxbljAICK88eg9G4PF+/2NR5z8XmoxmsoAMQczAYRmReLj8gE0i4W1dyy4NMB2vPXS057/jhcfB6C0uehGo8pW3Ss8RpanQFpucWWC00kUyw+IhMo0upq/LeLTwdUZJ+AvrwYhvIiOHvdD5f7H0RFThr05cWoyj9Xa8R3/XWqLBWZSLZYfEQm4KmsuTPI5X4/GCrKUHL4V7jc/yAAwMHlPji6e6Hk8K9wdPeCc72mN3gdZ4vkJZIzFh+RCfg19YSL0/9+nRycXeDSTIWifZtqjOyU3u1rPWb8mpMD/Jp5WCQvkZyx+IhMYGhX71qPufh0gKGsAC7e7f/22EMwlBVA+d9tDH8nARjapfbrEJFp8cgyIhMZ/91+/P7XpX89puxmFArgifZNsDy8m+mDEVENHPERmcjEPioonRzv6nuVTo6I6qMycSIiuhEWH5GJtPQAhrd1hqvznf1auTo7IPppP/h78y4NRJbA4iO6R6dOncKECRPQsGFDLJ8yCtFPPwhXZ0coFP/+fQoF4OrsiOinH+QB1UQWxLszEN2l3NxchIWF4cCBA6isrITBYEBkZCTCA1rB37selsamY/fJPChwfXN6ter78fVt1whRfVQc6RFZGIuP6C45Ojri/Pnz0Ol0MBgM8PDwwFNPPQUA8Peuh+Xh3XClpAIxB7ORlluM32LjoSsvxqtjh2NoF96BnUgUruokuge7d+/G448/DkdHR0iShMLCQiiVylrPkyQJnp6eKCsrw8aNGzFw4EABaYkI4DU+ort27tw5hIeHY926dRg6dCgCAgJuWHoAsHfvXuj1ehgMBkRERCAvL8/CaYmoGouP6C4UFRWhf//+eP311zFkyBCsXr0au3fvvunzV69ejYqK67ccKi0txZgxYywVlYj+gcVHdId0Oh1GjBiB3r1749VXXzU+7uBw81+nrVu3Gr/evHlz1KvHBS1EovAaH9EdkCQJkyZNwunTp7Ft2zY4Od3e+rDMzEy4urri8ccfx9dff40uXbqYOSkR3QxHfER34LPPPkNsbCw2bNhw26UHAK1atUKTJk2gUqmQnp5uxoREdCvczkB0m37++Wd8/PHHSEpKQt26de/qNVQqFTIyMkycjIjuBEd8RLfh0KFDGDt2LH766Se0bNnyrl+HIz4i8Vh8RLeQk5ODAQMGYNmyZejRo8c9vZavry+Lj0gwFh/RvygpKUFISAheeuklDB069J5fj1OdROJxVSfRTej1egwaNAiNGzfGypUrobjVqdO3+Zru7u64evUqXF1dTZCSiO4UR3xENzFlyhSUlpZi2bJlJik94Pr5nq1atcKZM2dM8npEdOdYfEQ3sHTpUuzYsQMxMTFwdnY26WtzgQuRWNzOQPQPO3bswAcffIDExETUr1/f5K/PBS5EYnHER/Q3R48exZgxY/Djjz+iTZs2ZnkPLnAhEovFR/Rfubm56N+/Pz777DMEBgaa7X041UkkFouPCEBZWRkGDBiA559/HiNGjDDre3Gqk0gsbmcg2TMYDAgLC4O7uztWrVplshWcN1NZWQlPT08UFRWhTp06Zn0vIqqNIz6SvWnTpiE/Px8rVqwwe+kBQJ06ddC8eXNkZWWZ/b2IqDYWH8naypUrsWnTJmzcuBEuLi4We18ucCESh9sZSLZ27tyJd955B/Hx8WjQoIFF35sLXIjEYfGRLJ04cQKjRo3CDz/8gAceeMDi788FLkTicKqTZOfSpUvo378/5s+fD41GIyQDpzqJxGHxkayUl5cjNDQUo0ePRnh4uLAcnOokEofbGUg2DAYDnnnmGTg6OmLNmjUWWcF5M+Xl5ahfvz5KS0vh6OgoLAeRHHHER7LxzjvvICcnB1999ZXQ0gMAV1dXNGzYEDk5OUJzEMkRi49kYdWqVVi3bh1++uknKJVK0XEAcLqTSBQWH9m92NhYvPnmm9i2bRsaNWokOo4RV3YSicHiI7t28uRJDB8+HGvXroWfn5/oODVwZSeRGCw+slv5+fkIDg7GnDlz8Oijj4qOUwtHfERisPjILlVUVGDQoEEICwvDs88+KzrODXHERyQGtzOQ3ZEkCeHh4aiqqsK6devg4GCdf98VFhbi/vvvR3FxsfBVpkRywiPLyO7MnDkTGRkZ2L17t9WWHgDUrVsXrq6uuHTpEpo2bSo6DpFsWO+nAtFdWLNmDVatWoXNmzfD1dVVdJxb4nQnkeWx+MhuxMfH49VXX8XWrVvRpEkT0XFuC/fyEVkei4/sQnp6OsLCwrBmzRo89NBDouPcNq7sJLI8Fh/ZvKtXryI4OBgzZ85Ev379RMe5I5zqJLI8Fh/ZtMrKSgwZMgQhISEYP3686Dh3jCM+IsvjdgayWZIk4dlnn0VhYSFiYmJs8i4HeXl58PPzw5UrV0RHIZINbmcgm/XRRx/h2LFj2LNnj02WHgA0bNgQOp0OV69ehZeXl+g4RLLA4iObtH79enz++edITk6Gm5ub6Dh3TaFQwNfXFxkZGSw+IgvhNT6yOcnJyZg0aRJ+/vlnNG/eXHSce8YFLkSWxeIjm3L27FkMHjwYq1atgr+/v+g4JsG9fESWxeIjm1FQUIDg4GBER0fj6aefFh3HZLiyk8iyWHxkE6qqqjB06FD069cPL730kug4JsWpTiLLYvGR1ZMkCVFRUVAqlViwYIHoOCbHqU4iy+KqTrJ68+bNQ2pqKuLj421228K/adasGQoLC1FaWmrTK1SJbAVHfGTVNm7ciM8++wxbt26Fh4eH6Dhm4eDggDZt2nC6k8hCWHxktVJTUzFhwgRs3rwZ3t7eouOYFRe4EFkOi4+s0rlz5xAaGoovvvgCXbp0ER3H7LjAhchyWHxkdYqKihAcHIwpU6Zg4MCBouNYBBe4EFkOi4+sik6nw/Dhw6FWq/HKK6+IjmMxnOokshwWH1kNSZLw8ssvAwA+++wzKBQKwYksh1OdRJbD7QxkNRYuXIi4uDgkJibCyUleP5otWrTAxYsXUVFRARcXF9FxiOwaR3xkFX7++WfMnTsXW7duhaenp+g4Fufk5AQfHx9kZmaKjkJk91h8JNyhQ4fw3HPPYdOmTWjZsqXoOMJwgQuRZbD4SKjs7GwMGDAAy5YtQ/fu3UXHEYoLXIgsg8VHwpSUlCAkJASTJk3CkCFDRMcRjgtciCyDxUdC6PV6jBw5Et26dcMbb7whOo5V4FQnkWWw+EiI119/HWVlZVi6dKmsti38G051ElmGvNaMk1VYsmQJfvvtNyQlJcHZ2Vl0HKvRunVrnDt3DjqdTnbbOYgsiSM+sqjt27dj1qxZ2LZtG+rVqyc6jlVRKpVo0qQJzp8/LzoKkV1j8ZHFHDlyBJGRkdi4cSNat24tOo5V4nQnkfmx+MgicnNzERISgkWLFqFXr16i41gtruwkMj8WH5ldaWkpQkJCMH78eAwfPlx0HKvGER+R+bH4yKwMBgNGjx6NDh06YPr06aLjWD2O+IjMj0vHyKzefPNNXL16FevWreO2hdvAvXxE5sfiI7NZsWIFNm/ejJSUFNSpU0d0HJvQpk0bZGRkQJIk/qFAZCac6iSz+P333zFjxgxs27YNXl5eouPYDA8PD3h6eiI3N1d0FCK7xeIjkzt+/DhGjRqFH374AQ888IDoODaH051E5sXiI5O6dOkS+vfvjwULFkCtVouOY5O4spPIvFh8ZDLl5eUYOHAgIiIiEB4eLjqOzeLKTiLzYvGRSRgMBkRERMDX1xfvvvuu6Dg2jVOdRObFVZ1kEu+88w4uXLiAnTt3cjXiPfL19eWIj8iMWHx0z77++musX78eKSkpUCqVouPYvOoRH7c0EJkHpzrpnuzevRvTpk3D1q1b0bBhQ9Fx7IKXlxcUCgWuXLkiOgqRXWLx0V07efIkRowYgXXr1sHPz090HLvCBS5E5sPio7uSn5+P4OBgzJkzB3379hUdx+5wgQuR+bD46I5ptVqEhoZi2LBhePbZZ0XHsUvcy0dkPiw+uiOSJOG5555D8+bNMWvWLNFx7BanOonMh8VHd+T9999HRkYGvvnmGzg48MfHXDjVSWQ+3M5At2316tX45ptvkJKSAldXV9Fx7Br38hGZj0KSJEl0CLJ+8fHxGDJkCGJjY9G+fXvRceyeJEnw8PDAhQsX4OnpKToOkV3hXBXdUnp6OsLCwrBmzRqWnoUoFArjvfmIyLRYfPSvrl69iuDgYMycORP9+vUTHUdWuMCFyDxYfHRTlZWVGDx4MAYMGIDx48eLjiM7XOBCZB4sProhSZIwfvx41K9fH//3f/8nOo4scS8fkXmw+OiGPvzwQxw7dgyrV6/mtgVBONVJZB7czkC1rF+/HitWrEBKSgrc3NxEx5EtTnUSmQe3M1ANycnJGDhwIHbu3Al/f3/RcWRNr9fDzc0N165d475JIhPiHBYZnTlzBoMHD8Y333zD0rMCjo6OaNWqFc6ePSs6CpFdYfERAKCgoADBwcF4++238dRTT4mOQ//FBS5EpsfiI1RVVWHo0KF44oknMHHiRNFx6G+4wIXI9Fh8MidJEqKiouDq6or58+eLjkP/wAUuRKbH4pO5uXPnYv/+/Vi7di0cHR1Fx6F/4FQnkelxO4OMbdy4EYsWLUJycjLc3d1Fx6Eb4FQnkelxO4NMpaamIjg4GL/++is6d+4sOg7dRGVlJTw9PVFcXAxnZ2fRcYjsAqc6ZSgrKwuhoaH44osvWHpWrk6dOmjWrBmysrJERyGyG5zqlJnCwkL0798fb7zxBgYMGCA6Dt2G6ulOlUolOgqRUX5JBWIOZCPtYhGKtDp4Kp3g19QTYV290cDdRXS8f8WpThnR6XTo378/fH19sXjxYigUCtGR6Da8+OKL6NChA7eakFU4fL4AS2LTsedUHgCgQmcwfk3p5AAJQJ92jRD1iAqdfOoJSvnvOOKTCUmSMHnyZCgUCixcuJClZ0N8fX25wIWswuqUTMzengatTo8bDZm0/y3B305cQtypfEQ/7YfwgFaWDXkbWHwy8emnnyIhIQEJCQlwcuL/7bZEpVIhLi5OdAySueul9xfKqwy3fK4kAeVVesze/hcAWF358RNQBrZs2YJ58+YhOTkZnp6eouPQHeJePhLt8PkCzN6eVqv0So/Hoih1E6quZMOhjiucm7RB3V7DoPR5CABQXmXA7O1p8PeuB39v65n2ZPHZuYMHD2LcuHHYtm0bWrRoIToO3YU2bdrg7NmzMBgMvDciCbEkNh1anb7GY0X7fkJhSgwaPDERytZdoHB0QvmZAyg/vddYfACg1emxNDYdy8O7WTr2TfG3yI5lZ2dj4MCBWL58Obp37y46Dt0lNzc3eHl5IScnR3QUkqH8kgrsOZVX45qeQVuKgvg18Hr8RdzXLhAOdZRQODrhvgd6ov6jY2t8vyQBu0/m4UpJhYWT3xyLz06VlJSgf//+mDx5MgYPHiw6Dt0jTneSKDEHsms9VnEhDZKuEve17XVbr6EAEHOw9uuIwuKzQ3q9Hs888wy6d++OKVOmiI5DJsCjy0iUtItFNbYsAIC+vAgO93lC4XB75/tqdQak5RabI95dYfHZoddeew1arRZLly7ltgU7wbs0kChFWl2txxxdPWEoK4Jk0N/gO272OlWmjHVPWHx2ZvHixdi5cyd++OEHnu1oR7iXj0TxVNZeA+nS3A8KJ2eUnUq+g9exns8jruq0I9u3b8eHH36IxMRE1KtnPUuH6d6pVCqcPn1adAySgby8PMTHxyMuLg7x8fHw6j0MLnX9a0x3OijdUC9oFK7+thwKB0coW3eGwsEJ2sw/oT13BPX71lzgonRygF8zD0v/U26KR5bZicOHD6Nfv37YvHkzevW6vQvOZDsKCgrg4+ODoqIiTl+TSZ07d85YdHFxccjNzUVgYCA0Gg3UajVa+/mj76cJta7zAUDJ8d0oTt2MqivnoajjCpemKnj2Gg6l94M1nufi5ICkNx+1mjM8OeKzA7m5uRgwYAAWL17M0rNT9erVg4uLC/Ly8tC4cWPRcchGSZKEU6dOGUsuPj4eZWVlxpJ74YUX4O/vX+um1I+0bYTf/7pU65gy94f6wv2hvv/6ngoF0LddI6spPYDFZ/NKS0sREhKCCRMmYNiwYaLjkBlVL3Bh8dHt0uv1OHLkSI2pS1dXV6jVamg0GkRHR6Ndu3a3nEV4CNnY41AHFfo7nyBUOjkiqo913VmExWfD9Ho9wsPD0aFDB7z11lui45CZVe/lCwwMFB2FrFRlZSX2799vHNElJSWhWbNm0Gg0CA0Nxfz589GyZcvbfj2dToepU6diy5YteOnT9fh83+XbOquzmquzA6Kf9rOq48oAFp9NmzZtGgoKCrB+/Xpe95EB7uWjfyotLUVycrJxNJeamop27dpBrVZj3LhxWLVq1V3PEOTn52P48OFwcnLCvn374OXlhYYN//3uDNUUiusjPd6dgUxqxYoV2LJlC5KTk1GnTh3RccgCVCoVduzYIToGCXT16lUkJiYaR3THjx/Hww8/DI1Gg6lTpyIwMBB169a95/f5888/MWjQIAwfPhyzZ882XvMLD2gFf+96WBqbjt0n86DA/25FBPzvfnx92zVCVB+V1Y30qnFVpw367bffMGbMGCQkJPCu3DKSlJSEV199FXv37hUdhSzkwoULNa7PZWZmIiAgwHiNrkePHnB1dTXpe65duxaTJ0/G4sWLMXz48Js+70pJBWIOZiMttxhF2ip4Kp3h18wDQ7vwDuxkYseOHcOjjz6KjRs3IigoSHQcsqDLly+jffv2yM/PFx2FzECSJJw5c8ZYcnFxcbh27RqCgoKMqy47d+5stoMpdDod3nrrLfz444/46aef0KlTJ7O8jzXgVKcNuXTpEkJCQvDJJ5+w9GSoUaNGqKioQEFBAQ8osAMGgwEnTpyosbUAgLHkXnvtNbRv394it6K6evUqRowYAUmSkJqaigYNGpj9PUVi8dmI8vJyDBw4EBERERg1apToOCSAQqEwLnDp2rWr6Dh0h6qqqnDo0CFjySUkJMDLywtqtRpPPvkkPvzwQ7Ru3driC9WOHDmCQYMGYfDgwfjoo4/g5GT/tcCpThtgMBgwfPhwuLi44LvvvuMKThkLCwvD0KFD//XaC1mH8vJy7Nu3zzii27t3L1q1amUc0anVajRv3lxoxg0bNmDixIlYuHAhRo4cKTSLJdl/tduBt99+GxcvXsTOnTtZejLH+/JZr8LCQiQlJRlHdIcOHULHjh2hVqsxefJk9O7dG15eXqJjAri+Bzg6Ohrr1q3Db7/9hs6dO4uOZFEsPiv39ddfY8OGDUhJSYGLi3WvlCLzU6lUSEpKEh2DcH2xUXx8vHEhyqlTp9C9e3doNBq89957CAgIgLu7u+iYtVy9ehUjR440bnZv2LCh6EgWx+KzYrt27cK0adMQFxcnyx9Oqk2lUuHbb78VHUOWsrKyamwtyM3NRe/evaFWq7F48WJ07drV6v84PXr0KAYNGoSQkBDMnTtXFtfzboTX+KxUWloaHnnkEaxfvx59+vQRHYesxPnz59GzZ09cuHBBdBS7JkkSTp48WWNrgVarNe6f02g06NixY63DnK1ZTEwMXnzxRSxYsACjR48WHUcoFp8VysvLQ69evfD2228jMjJSdByyIgaDAW5ubsjPz4ebm5voOHZDr9fj8OHDNUZ09913n3EhikajQdu2bW3yGrter8eMGTOwevVqbNy4kSuCwalOq6PVahEaGorhw4ez9KgWBwcHtG7dGmfOnEHHjh1Fx7FZFRUVxsOc4+PjkZSUhObNm0OtVmPw4MH45JNP0KJFC9Ex71lBQQFGjhyJsrIypKam8s4e/8URnxWRJAmjRo2CwWDA999/b5GNq2R7BgwYgLFjxyI0NFR0FJtRUlKC5ORk44hu//79aNeunXFEFxQUZHelcOLECYSGhuKpp57CvHnzzHbiiy3iiM+KvPfeezh79ix27drF0qObqr4vH93c1atXkZCQYBzRHT9+HJ07d4ZGo8G0adMQGBgIT09P0THN5qeffsL48eMxb948REREiI5jdVh8VuK7777Dd999h5SUFJMfOkv2xdfXF0ePHhUdw6pcuHChxkKUrKwsBAQEQKPRYO7cuejRoweUSqXomGZnMBjw3nvvYdWqVdi+fTu6d+8uOpJVYvFZgbi4OEyZMgW7d++2u+kWMj2VSoVNmzaJjiGMJEnIyMiosRDl2rVrxkUokZGR6Ny5s+yW6hcWFiI8PByFhYVITU1FkyZNREeyWvL6ybBCp0+fxrBhw7BmzRq0b99edByyAXKb6jQYDDh+/HiNEZ1CoTBuK5gyZQoefPBBWV8eSEtLQ2hoKB577DEsWLCA9+i8BS5uEejKlSvo1asX3njjDTz//POi45CNqKqqgru7O4qLi+3yA66qqgoHDx40llxiYiK8vLyMRadWq4Uc5myttmzZgnHjxmHOnDkYO3as6Dg2gcUnSGVlJfr164eePXvi448/Fh2HbIxKpcL27dvRtm1b0VHuWXl5Ofbu3Wsc0aWkpKBNmzY1DnNu1qyZ6JhWx2Aw4IMPPsAXX3yBmJgY9OzZU3Qkm8GpTgEkScLzzz+Phg0bYs6cOaLjkA2qPqzaFouvsLAQiYmJxhHdn3/+iY4dO0Kj0eDll1/Ghg0bUL9+fdExrVpRURFGjx6N/Px8pKamomnTpqIj2RQWnwCzZ8/GiRMnsGfPHllfl6C7Z0vX+aoPc64e0Z06dQo9evSARqPBzJkzERAQwFNo7sDJkycRGhqKPn364IcffrDL6W5zY/FZ2Nq1a7Fy5UqkpKTgvvvuEx2HbFT1DWmtUVZWVo2FKBcvXkTv3r2h0WiwZMkSdO3alR/Wd2nr1q0YO3YsZs+ezXUB94DFZ0FJSUl4+eWX8ccff/CaBd0TX19f7Ny5U3QMSJKEtLS0GiM6rVZrXIgSFRVlc4c5WyODwYAPP/wQy5Ytw6ZNmxAYGCg6kk1j8VnImTNnMGTIEHzzzTc8Y5HumagRX/VhztUlFx8fDzc3N2g0GvTp0wczZszAAw88wBWXJlRcXIyIiAjk5uYiNTVV+F3b7QFXdVrAtWvXEBgYiEmTJiEqKkp0HLIDWq0W9erVQ2lpqVlHUxUVFUhNTTWO6JKTk9G8efMaWwt8fHzM9v5yd/r0aYSGhiIwMBCLFy+2+vv92QoWn5lVVVXhySefhL+/Pz755BPRcciO+Pj4ID4+Hq1atTLZa1Yf5lw9otu/fz/8/PxqHObcqFEjk70f3dwvv/yCiIgIzJw5ExMmTOAo2oQ41WlGkiThxRdfhJubG+bNmyc6DtmZ6unOeym+K1euICEhwTiiO3HiBLp06QK1Wo233noLvXr1suvDnK2RJEmYM2cOFi1ahI0bNyIoKEh0JLvD4jOjjz/+GAcPHkRcXBwv7pPJVe/l+89//nPb35OTk1NjIUpWVhZ69eoFjUaDefPmyeYwZ2tVUlKCZ599FufOncO+ffvg7e0tOpJdYvGZSUxMDBYvXoyUlBS4u7uLjkN26FZ7+aoPc/771oLCwkLjaShjx47Fww8/LLvDnK1VRkYGQkND0b17d+zZs4d/gJgRf+LNYN++fYiKisKvv/6K+++/X3QcslMqlQr79u0z/rfBYMCxY8dqjOgcHBx4mLMN+PXXXzFmzBjMmDEDUVFRvJ5nZlzcYmJZWVkIDAzE8uXLERISIjoO2bF9+/Zh5MiRmDBhgvEw54YNGxoXomg0GrRq1YofolZMkiTMnTsXn3zyCdavXw+NRiM6kixwxGdChYWFCA4OxtSpU1l6ZHJlZWXYu3evcUSXkpKC8vJyZGVlYcyYMVi5ciXPbLQhpaWleO6555CRkYF9+/ZxW4gFccRnIjqdDv3794dKpcKiRYv4Vzbds+rDnOPi4hAXF4cjR44YD3NWq9Xo3bs3/Pz88Oeff/IkIBtz9uxZhIaGonPnzli2bBlcXV1FR5IVjvhMQJIkTJo0CQ4ODvj0009ZenRXLl26ZDwNJS4uDunp6ejRowfUajVmzZqFnj171jrMuXplJ4vPduzcuROjRo1CdHQ0Jk2axM8LAVh8JvDpp58iMTERCQkJXCFHt0WSJGRlZdVYiHLp0iXjYc5Lly69rcOcq/fyqdVqCyWnuyVJEubPn4/58+dj/fr16NOnj+hIssVP6Xu0efNmzNFdsoYAABX7SURBVJ8/H0lJSdzoSzdVfZhz9bRlfHw8KisrjdOWL730Ejp06HDH+z2rR3xk3crKyjBu3DikpaUhJSUFLVu2FB1J1lh89+DAgQMYN24cfvnlF7Ro0UJ0HLIiOp0Ohw8frjGi8/DwgFqtxqOPPor33nsPKpXqnqe5VCoVfv75ZxOlJnPIzMzEoEGD0KFDByQmJvJ6nhVg8d2l7OxsDBw4ECtWrEC3bt1ExyHBqg9zrh7RJScnw9vbGxqNBkOHDsXChQvNsmrPmu/LR8CuXbswcuRIvPnmm3jllVd4Pc9KcFXnXSguLoZarUZ4eDimTJkiOg4JUFxcXOMw5wMHDuDBBx807p8LCgpCw4YNzZ4jPz8fKpUK165d44eqFZEkCQsXLsScOXOwZs2aOzpWjsyPxXeHdDodQkND0bx5c3z++ef8sJGJ6sOcq0d0f/31F7p06WI8FaVXr17w8PCweC5JklC/fn1kZGSgQYMGFn9/qq28vBzjx4/H0aNHsWnTJpPePYNMg1Odd+i1115DZWUllixZwtKzY9nZ2TWuz507dw6BgYHQaDRYsGABunfvbhVnKSoUCuN0J4tPvHPnzmHQoEFo27YtkpKScN9994mORDfA4rsDixYtwh9//IHExEQ4OzuLjkMmIkkS0tPTaxzmXFRUZJy2HDduHDp16mS1W1WqV3b26NFDdBRZ27NnD0aMGIHXX38dr7/+Ov8wtmLW+ZtshbZt24aPPvoIiYmJqFevnug4dA+qD3P++9YCJycn47Tl1KlT4efnZzOHOXOBi1iSJGHx4sWYNWsWVq9ejX79+omORLfA4rsNhw8fxrPPPostW7agdevWouPQHaqqqsKBAweMJZeQkIDGjRtDrVajf//++Pjjj9GyZUub/Qvd19cXe/bsER1DlrRaLV544QUcPHgQycnJaNOmjehIdBtYfLdw4cIFhISEYPHixQgICBAdh25DWVkZUlJSjNOW+/btg6+vLzQaDSIiIuzuMGeVSoUvv/xSdAzZyc7OxuDBg9GqVSskJyfXOk6OrBdXdf6L0tJSaDQaDBkyBNOnTxcdh26ioKDAeJhzfHw8Dh8+jE6dOtU4zNmep6cvXLiALl264OLFi6KjyEZ8fDyGDx+Ol19+GVOnTrXZ2QK5YvHdhF6vx5AhQ1C/fn189dVX/MG2IhcvXqxxmHNGRgZ69OhhvEbXs2dPWa2mkyQJbm5uuHTpkpAtFXIiSRKWLVuG999/H9988w2efPJJ0ZHoLnCq8ybefPNNFBYWYsOGDSw9gSRJQmZmZo2tBZcvX0ZQUBA0Gg2WL1+OLl263PIwZ3umUCjg6+uLjIwMPPzww6Lj2C2tVouJEydi7969SExMhEqlEh2J7hKL7wY+//xzbN26FUlJSbL+QBVBkiT89ddfNbYW6HQ647Tl5MmT0aFDB5tZcWkpLD7zysnJweDBg+Ht7Y3k5GSOrG0ci+8ffv31V7z33ntISEiAl5eX6Dh2T6fT4c8//zSWXEJCAjw8PKDRaPCf//zHZIc52zuVSsW7NJhJYmIihg0bhqioKEyfPp0/i3aAxfc3x44dw+jRo7Fx40b4+vqKjmOXtFptjcOcU1JS4OPjA41Gg2HDhmHRokXw9vYWHdPmqFQqHDhwQHQMu/P555/jnXfewddff43g4GDRcchEWHz/denSJYSEhODTTz9FUFCQ6Dh2o7i4GElJScYR3cGDB9G+fXuo1WpERUXh+++/51FbJuDr64sNGzaIjmE3KioqMGnSJCQkJCAxMREPPPCA6EhkQiw+XN/3NWDAAERGRmLkyJGi49i0/Pz8Goc5p6WloWvXrtBoNHj77beFHeZs7zjVaTq5ubkYMmQImjRpgr179/Ln1Q7JfjuDwWDA8OHDoVQq8e2333L+/g6dP3++xorL7OxsBAYGGs+57N69O1xcXETHtHs6nQ7u7u4oKCiwisOzbVVycjLCwsIwYcIEREdHcxGVnZL9iC86OhoXL17Ezp07WXq3IEkSTp8+bSy6uLg4lJSUGEtu/Pjx8Pf3t9rDnO2Zk5MTWrRogbNnz+LBBx8UHccmffHFF5g+fTq+/PJLhISEiI5DZiTrT6ivvvoKMTExSE5O5qjkBvR6fa3DnOvUqWPcWjBt2jT4+fnxDwYrUT3dyeK7M5WVlXj55Zexe/duxMfHo127dqIjkZnJtvh27dqF6dOnIy4uziJ3yrYFlZWVNQ5zTkxMRJMmTaBWqzFgwADMmzcPLVu2FB2TbqJ6Lx/dvosXL2Lo0KHw8vLC3r17UbduXdGRyAJkWXxpaWl45plnsGHDBrRt21Z0HGGqD3OuHtGlpqZCpVJBo9EgMjISX375JZo0aSI6Jt0mlUqF06dPi45hM/bt24chQ4bgueeew4wZM3g9T0ZkV3x5eXkIDg7Gxx9/jEceeUR0HIu6du1ajcOcjxw5gocffhhqtRpTpkxBYGCgXR/mbO9UKhV++eUX0TFswtdff42pU6di5cqVCA0NFR2HLExWxafVahEaGopnnnkGERERouOYXW5ubo3DnM+cOYOePXtCo9Hgo48+Qo8ePWR1mLO941TnrVVVVeHVV1/Fb7/9hj179qB9+/aiI5EAstnOIEkSRo4cCUmS8P3339vdtEb1Yc5/X4iSn5+PoKAg46rLLl26wNnZWXRUMpOKigp4enqitLSUK2tv4PLlywgLC4OHhwdWr17N2Q0Zk81vx7vvvovMzEzs2rXLLkrPYDDUOszZYDAYS+6VV17BQw89ZBf/Vro9Li4uaNasGc6dO8c7gf/D/v37MXjwYEREROD999/n74XMyaL4vv32W6xevRopKSlwdXUVHeeu6HQ6HDp0qMZhznXr1oVGo0G/fv0wc+ZM+Pr6cmuBzPn6+iI9PZ3F9zfffvstXn/9dXz++ecYPHiw6DhkBey++OLi4jBlyhTExsaicePGouPcNq1Wi3379hlHdMnJyWjZsiXUajVGjBiBJUuW4P777xcdk6xM9V6+xx9/XHQU4aqqqvDGG29g69at2L17Nzp06CA6ElkJuy6+06dPY9iwYfj++++t/iJ2UVFRjcOcDx06hPbt20Oj0WDixIk8zJluCxe4XJeXl4dhw4ZBqVQiNTUV9evXFx2JrIjdFt+VK1cQHByMWbNm4bHHHhMdp5a8vDzjYc7x8fFIS0tDt27doNFoMGPGDPTq1Qvu7u6iY5KNUalUSExMFB1DqIMHD2Lw4MEYOXIkPvjgAzg6OoqORFbGLouvoqICgwYNwqBBgzBu3DjRcQBcP8z57wtRcnJyEBgYCI1Gg4ULF6Jbt248No3umdzv0rBmzRq88sorWLp0KcLCwkTHIStld9sZJElCREQESktL8cMPPwhZvSVJEk6dOlXjrgUlJSXQaDTGcy47derEv0TJ5EpKStC4cWOUlJTIauWiTqfDm2++iU2bNmHTpk3o2LGj6EhkxexuxDdr1iykpaUhNjbWYr/4er0eR48erTGic3FxMRbd9OnT0a5dO664JLNzd3dH3bp1ceHCBdncyT4/Px8jRoyAo6MjUlNT4eXlJToSWTmbKb78kgrEHMhG2sUiFGl18FQ6wa+pJ8K6eqOB+/UpwrVr1+LLL79ESkqKWU8kqaysxP79+40ll5iYiKZNm0Kj0WDgwIE8zJmEqp7ulEPxHT58GIMGDUJYWBg+/PBDzqLQbbH6qc7D5wuwJDYde07lAQAqdAbj15RODpAA9GnXCEFeZXg9cij++OMPk09zlJaW1jrMuW3btsZpy6CgIB7mTFYjMjISarUazz33nOgoZrVu3TpMmjQJixYtwogRI0THIRti1SO+1SmZmL09DVqdHjeqZ+1/S/C345ewQ1eB5+asMknpXbt2DQkJCcYR3dGjR/Hwww9Do9Fg6tSpCAwM5O1LyGrZ+wIXvV6Pt956CzExMdi5cyc6deokOhLZGKstvuul9xfKqwy3fK4EQOHkgk2ZCjyUkonwgFY1vp6dnf2v0z7VhzlXj+gyMzONhznPmTMHPXv2tNkTX0h+fH198dNPP4mOYRZXr17FiBEjYDAYkJqayr2tdFessvgOny/A7O1pNUove+lYGMoKAIUDFHWUcG3dFV6PvwCHOv8rpPIqA2ZvT4O/dz34e9eDJEmYM2cOoqOjcfjwYXTs2BGSJOHs2bM1DnO+cuUK1Go11Go1xowZg86dO/MwZ7JZ9jriO3LkiHGb0pw5c3gQN901q7zGN/67/fj9r0s1pjezl45Fg6cnw7XVw9CXXMOl9e/AVdUD9R8ZU+N7FQrgifZN8Nkwf0RGRmLz5s2oqqpCaGgoHB0dERcXB0mSamwt4GHOZE+uXbuGli1borCw0G5WEm/YsAETJ07EwoULMXLkSNFxyMZZ3Z9M+SUV2HMq74bX9Ko5uteHa5suqLp8ptbXJAnYlXYZvu074cLZUzAYro8ak5KS8MEHH2DWrFlo06aN3XwgEP1T/fr14eTkhPz8fDRq1Eh0nHui1+sRHR2NdevW4ddff0WXLl1ERyI7YHXFF3Mg+5bP0RXlo/zMAShb+N/4CZIEqVUPOJ7LgKurK6qqqlBcXIzIyEgWHslC9XSnLRfftWvX8Mwzz6CyshKpqak2/W8h62J1xZd2sajGloW/y/txFqBQQKosh7KlP+qpR93weZUGIOz5VzF3x5c4evQoEhIScPLkSXPGJrIq1cXXq1cv0VHuyrFjxxAaGoqQkBDMnTuX1/PIpKzup6lIq7vp1xoNeRuurR6G9txR5G+ZC315ERyUNz7IuUhbBScnJ3Tu3BmdO3c2V1wiq2TLd2n48ccf8cILL2DBggUYPXq06Dhkh6xuRYen8tZdrGzREW4dH8O1XV/9y+twVSbJly2u7NTr9Xj77bfx2muvYceOHSw9MhurKz6/pp5wcbp1LM/uA6HNPITKS7UXuCidHODXzMMc8Yhsgq2N+AoKCjBgwAAkJCQgNTUVXbt2FR2J7JjVFd/Qrrd3vqDjfXXh1uFRFCauq/U1CcDQLvZ/TiHRzdjSiO/EiRPo0aMHfH198fvvv6Nx48aiI5Gds5l9fLereh/f8vBupg9GZCMkSYKHhwdycnKs+ni9TZs24fnnn8fcuXMRGRkpOg7JhNWN+ABgYh8VlE53d8q60skRUX1UJk5EZFsUCoVVT3caDAa8++67mDx5MrZv387SI4uyyuLr5FMP0U/7wdX5zuK5Ojsg+mk/+HvXM1MyItthrdOdhYWFCA0Nxa5du5Camoru3buLjkQyY5XFBwDhAa0Q/fSDcHV2xK32nCsUgKuzI6KffrDWAdVEcmWNxZeWloaePXvCx8cHf/zxB2/nRUJYbfEB18tv/fgAPNG+CVycHKD8x2pPpZMDXJwc8ET7Jlg/PoClR/Q31jbV+fPPP0Oj0eCNN97AkiVLUKdOHdGRSKasbgP7P/l718Py8G64UlKBmIPZSMstRpG2Cp5KZ/g188DQLv+7AzsR/Y9KpcKaNWtEx4DBYMCsWbOwcuVKbNmyBQEBAaIjkcxZ5apOIrp3WVlZCAwMRE5OjrAMRUVFiIiIwOXLlxETE4NmzZoJy0JUzaqnOono7nl7e+PKlSsoKysT8v6nTp1CQEAAmjRpgt27d7P0yGqw+IjslKOjI1q3bo0zZ2qfbmRu27ZtQ1BQEF555RUsX76c1/PIqrD4iOyYpRe4SJKE2bNnY/z48di0aRPGjx9vsfcmul1Wv7iFiO6eJbc0lJSUIDIyEjk5OUhNTUXz5s0t8r5Ed4ojPiI7ZqniS09PR0BAAOrXr4/Y2FiWHlk1Fh+RHbPEVOeOHTvQu3dvvPTSS1ixYgVcXLi9iKwbpzqJ7Jg5R3ySJOH//u//sGjRIvz4448ICgoyy/sQmRr38RHZscrKSnh4eKC4uNikKytLS0sxduxYZGZm4scff4S3N28DRraDU51EdqxOnTq4//77kZWVZbLXPHPmDHr16gU3Nzfs2bOHpUc2h8VHZOdMOd35+++/IzAwEBMmTMCXX34JpVJpktclsiRe4yOyc6ZY4CJJEubPn48FCxZgw4YN0Gg0JkpHZHksPiI7d68jvrKyMowbNw6nTp3C3r174ePjY8J0RJbHqU4iO3cvxZeZmYnevXvD2dkZ8fHxLD2yCyw+Ijt3J1Od+/fvR3BwMCoqKvDHH38gICAAkZGRWLVqFVxdXc2clMgyuJ2ByM6VlZXBy8sLpaWlcHR0/NfnPvbYY4iNjUWnTp2Qk5ODtWvXom/fvhZKSmQZLD4iO5dfUoFnZy5H685qlFTq4al0gl9TT4R1rXkT57S0NHTu3BlarRYODg6YMWMG3n33XYHJicyDxUdkpw6fL8CS2HTsOZUHAKjQGYxfUzo5QALQp10jRD2iQiefeggNDcXmzZsBAAqFAgqFAunp6WjdurWI+ERmw+IjskOrUzIxe3satDo9/u03XKEAlE6OeOOxNnjukXZwdHSEv78/+vbtC7VajZCQkFtOjxLZGhYfkZ25Xnp/obzKcOsn/5erswMiO9XFlNCeLDqyeyw+Ijty+HwBRqxMQXmVvsbj2vPHURD7NSrzz0GhcIBzAx/Uf+x5uDRra3yOq7Mj1o8PgL93PUvHJrIobmAnsiNLYtOh1dUsPUNFGS7HzESDJ6Jwn18QJL0OFdnHoXB0rvE8rU6PpbHpWB7ezZKRiSyOxUdkJ/JLKrDnVF6ta3pVV3MAAG7tHwEAKBwc4dq6S63vlyRg98k8XCmpqLHak8jecAM7kZ2IOZB9w8edve6HQuGA/K0LUJ6xH3ptyU1fQwEg5uCNX4fIXrD4iOxE2sWiGlsWqjm43Iem4R8DUODKjkXIXjgSl2NmQl96rdZztToD0nKLLZCWSBxOdRLZiSKt7qZfc27og4b9XwUAVF05j/yf5+PqzpVoNHDqDV6nymwZiawBR3xEdsJTeXt/xzo38IFbx/+gKv/GN6f1VDrf8HEie8HiI7ITfk094eJU+1e66sp5FO3dCF1RPgBAV5SHshNxcGnertZzlU4O8GvmYfasRCJxqpPITgzt6o1Pdp6q9biijisqck+hKHUTDBWlcHBxg6uqB+r3HVvruRKAoV28LZCWSBwWH5GdaOjugkfaNsLvf12qsaXByaMhGoVOu+X3KxRA33aNuJWB7B6nOonsyMQ+Kiid7u7IMaWTI6L6qEyciMj6sPiI7Egnn3qIftoPrs539qvt6uyA6Kf9eFwZyQKnOonsTHhAKwC4o7szRD/tZ/w+InvHQ6qJ7NSR7AIsjU3H7pN5UOD65vRq1ffj69uuEaL6qDjSI1lh8RHZuSslFYg5mI203GIUaavgqXSGXzMPDO3izYUsJEssPiIikhUubiEiIllh8RERkayw+IiISFZYfEREJCssPiIikhUWHxERyQqLj4iIZIXFR0REssLiIyIiWWHxERGRrLD4iIhIVlh8REQkKyw+IiKSFRYfERHJCouPiIhkhcVHRESywuIjIiJZYfEREZGssPiIiEhWWHxERCQrLD4iIpKV/wcWnaFhYgKpSgAAAABJRU5ErkJggg==\n","text/plain":["
"]},"metadata":{"tags":[]}}]},{"cell_type":"markdown","metadata":{"id":"Ho61X-iE0fET"},"source":["Next, we define the conditional probability tables for all nodes, in this case: C,S,R,W\n","\n","We then add them to the Bayesian network and cross-check that the netowrk is well defined."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"6W0wgisEyvaS","executionInfo":{"status":"ok","timestamp":1611834676467,"user_tz":-60,"elapsed":858,"user":{"displayName":"Ulrich Kerzel","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64","userId":"10582137377195463323"}},"outputId":"2a247281-010d-4bd6-8a6f-4efd25313917"},"source":["#define CPT\n","\n","#cloudy\n","CPT_C = TabularCPD(variable='C', variable_card=2, \n"," values=[[0.5], [0.5]])\n","\n","#sprinkler\n","CPT_S = TabularCPD(variable='S', variable_card=2,\n"," values=[ [0.5, 0.9],\n"," [0.5, 0.1]\n"," ],\n"," evidence=['C'], evidence_card=[2]\n"," )\n","\n","#rain\n","CPT_R = TabularCPD(variable='R', variable_card=2,\n"," values=[ [0.8, 0.2],\n"," [0.2, 0.8]\n"," ],\n"," evidence=['C'], evidence_card=[2]\n"," )\n","\n","# wet grass\n","CPT_W = TabularCPD(variable='W', variable_card=2,\n"," values=[ [1.0, 0.1, 0.1, 0.001],\n"," [0.0, 0.9, 0.9, 0.99]\n"," ],\n"," evidence=['S','R'], evidence_card=[2,2]\n"," )\n","model.add_cpds(CPT_C, CPT_S, CPT_R, CPT_W)\n","model.check_model()"],"execution_count":5,"outputs":[{"output_type":"execute_result","data":{"text/plain":["True"]},"metadata":{"tags":[]},"execution_count":5}]},{"cell_type":"markdown","metadata":{"id":"Q4QIQeXO02BA"},"source":["Once the Bayesian network is fully defined we can use it for inference.\n","\n","In the example below, we with to estimate the probability for rain (R) given that we observe that the grass is wet (W=1).\n","\n","Exerise: Experiment and infer the values of different quantities given some observations."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"HTeS6XwEyvaV","executionInfo":{"status":"ok","timestamp":1611834713348,"user_tz":-60,"elapsed":693,"user":{"displayName":"Ulrich Kerzel","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64","userId":"10582137377195463323"}},"outputId":"d757422b-0004-4540-e866-3db0bc59fe13"},"source":["# Infering the posterior probability \n","from pgmpy.inference import VariableElimination\n","\n","infer = VariableElimination(model)\n","posterior_p = infer.query(variables=['R'], evidence={'W': 1})\n","print(posterior_p)"],"execution_count":8,"outputs":[{"output_type":"stream","text":["Finding Elimination Order: : 100%|██████████| 2/2 [00:00<00:00, 238.02it/s]\n","Eliminating: C: 100%|██████████| 2/2 [00:00<00:00, 376.69it/s]"],"name":"stderr"},{"output_type":"stream","text":["+------+----------+\n","| R | phi(R) |\n","+======+==========+\n","| R(0) | 0.2921 |\n","+------+----------+\n","| R(1) | 0.7079 |\n","+------+----------+\n"],"name":"stdout"},{"output_type":"stream","text":["\n"],"name":"stderr"}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"SlMEg9fAyvaZ","executionInfo":{"status":"ok","timestamp":1611834830851,"user_tz":-60,"elapsed":675,"user":{"displayName":"Ulrich Kerzel","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64","userId":"10582137377195463323"}},"outputId":"c3a7a572-c3ad-402d-ebb1-ea9c5a94d8c6"},"source":["# independence conditions\r\n","model.get_independencies()"],"execution_count":9,"outputs":[{"output_type":"execute_result","data":{"text/plain":["(C _|_ W | S, R)\n","(S _|_ R | C)\n","(R _|_ S | C)\n","(W _|_ C | S, R)"]},"metadata":{"tags":[]},"execution_count":9}]}]} \ No newline at end of file +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "uYCyEcvDyvaG" + }, + "source": [ + "# Inference in Bayesian Networks\n", + "\n", + "This notebook illustrates the Bayesian Network for the \"Wet Grass\" example discussed in the course book \"Inference & Causality\", unit 1.2, based\n", + "on: Murphy, K. (2001). An introduction to graphical models.Rap. tech,96,1–19" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 9618, + "status": "ok", + "timestamp": 1611834585233, + "user": { + "displayName": "Ulrich Kerzel", + "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64", + "userId": "10582137377195463323" + }, + "user_tz": -60 + }, + "id": "8D7_xzjnyvaI", + "outputId": "09460824-4c8c-47f6-ec77-24ec6b5cd4ef" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Requirement already satisfied: pgmpy in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (0.1.22)\n", + "Requirement already satisfied: scikit-learn in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (1.0.2)\n", + "Requirement already satisfied: joblib in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (1.1.0)\n", + "Requirement already satisfied: statsmodels in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (0.13.2)\n", + "Requirement already satisfied: tqdm in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (4.64.1)\n", + "Requirement already satisfied: scipy in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (1.9.1)\n", + "Requirement already satisfied: pandas in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (1.4.4)\n", + "Requirement already satisfied: pyparsing in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (3.0.9)\n", + "Requirement already satisfied: opt-einsum in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (3.3.0)\n", + "Requirement already satisfied: networkx in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (2.8.4)\n", + "Requirement already satisfied: torch in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (2.0.0)\n", + "Requirement already satisfied: numpy in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pgmpy) (1.21.5)\n", + "Requirement already satisfied: python-dateutil>=2.8.1 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pandas->pgmpy) (2.8.2)\n", + "Requirement already satisfied: pytz>=2020.1 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from pandas->pgmpy) (2022.1)\n", + "Requirement already satisfied: threadpoolctl>=2.0.0 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from scikit-learn->pgmpy) (2.2.0)\n", + "Requirement already satisfied: packaging>=21.3 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from statsmodels->pgmpy) (21.3)\n", + "Requirement already satisfied: patsy>=0.5.2 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from statsmodels->pgmpy) (0.5.2)\n", + "Requirement already satisfied: sympy in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from torch->pgmpy) (1.10.1)\n", + "Requirement already satisfied: jinja2 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from torch->pgmpy) (2.11.3)\n", + "Requirement already satisfied: filelock in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from torch->pgmpy) (3.6.0)\n", + "Requirement already satisfied: typing-extensions in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from torch->pgmpy) (4.3.0)\n", + "Requirement already satisfied: six in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from patsy>=0.5.2->statsmodels->pgmpy) (1.16.0)\n", + "Requirement already satisfied: MarkupSafe>=0.23 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from jinja2->torch->pgmpy) (2.0.1)\n", + "Requirement already satisfied: mpmath>=0.19 in /Users/kathi/opt/anaconda3/lib/python3.9/site-packages (from sympy->torch->pgmpy) (1.2.1)\n" + ] + } + ], + "source": [ + "!pip install pgmpy\n", + "\n", + "from pgmpy.models import BayesianModel\n", + "from pgmpy.factors.discrete import TabularCPD\n", + "import networkx as nx\n", + "import pylab as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "eS15uTMlzcKT" + }, + "source": [ + "First we define the Bayesian network by specifying which nodes are connected to which.\n", + "The following shortcuts are used:\n", + "* C: Cloudy\n", + "* S: Sprinkler\n", + "* R: Rain\n", + "* W: Wet grass\n", + "\n", + "We use the module ```networkx``` do visualize the network. The output is perhaps a bit counter-intuitive - but careful investigation shows that the arrows point indeed into the correct nodes:\n", + "\n", + "* *Cloudy* is a common cause to both *Sprinkler* and *Rain*\n", + "* *Sprinkler* is a cause to *Wet grass*\n", + "* *Rain* is a cause to *Wet grass*." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 319 + }, + "executionInfo": { + "elapsed": 748, + "status": "ok", + "timestamp": 1611834661134, + "user": { + "displayName": "Ulrich Kerzel", + "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64", + "userId": "10582137377195463323" + }, + "user_tz": -60 + }, + "id": "8r5E1lCGyvaM", + "outputId": "5a3fd0bf-b9cc-4351-b945-856343200b2a" + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/kathi/opt/anaconda3/lib/python3.9/site-packages/pgmpy/models/BayesianModel.py:8: FutureWarning: BayesianModel has been renamed to BayesianNetwork. Please use BayesianNetwork class, BayesianModel will be removed in future.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAHzCAYAAACe1o1DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhc0lEQVR4nO3dfZTVhX3n8e+9d2CGASRhKDhRBtIAJhF0G2ox2ZrI1tqeNnFXE59WbR4sUVs3idvk1MP22JN6jiZttrhpTqM1PbqrHFGrtCbd7aZpqbE1UI5uIrpRhjUyGBEiRBhmmIGZe/cPhPAwM8xwn34Pr9d/zH36jYczvvl9fvdOoVKpVAIAAE5RsdkHAABAuglKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqghKAACqIigBAKiKoAQAoCqCEgCAqrQ0+wAAkqpvcChe2dUXB4bKMbmlGPM7psbUVj82AY7nJyPAUbp39MbqDT2x7qWd0bO7PypH3VaIiK6Z7bH8rNlxzbKuWDhnerMOEyBRCpVKpXLyuwFk27bd/bFy7aZ4assbUSoWYrg8+o/Gw7dfsGBW3HHpkpg7s72BRwqQPIISyL01G3viD594IYbKlTFD8nilYiFaioX44iVnx1XnddXxCAGSTVACufa1dd3xlW9vrvp5Pn/xorh5+cIaHBFA+niXN5Bbazb21CQmIyK+8u3N8fDGnpo8F0DaOEMJ5NK23f1x0aonY3CoPOLtg6+9FHu+92gc2PH/Yrjvp1FsnRYtb5sTrWe8J2b+ym+P+JjWlmJ855YPuaYSyB1nKIFcWrl2UwyNcr1k/5aN8foDX4jKgf54+4WfiDlX3h4zL1oRbWe+N/pffGrU5xwqV2Ll2k31OmSAxHKGEsid7h298at3fXfU219ffWsM9+6Kd3z67igUS8fcVqmUo1AY+9/i37nlg7Fgto8UAvLDGUogd1Zv6IlSsTDq7eX9e6PYftoJMRkRJ43JUrEQD653LSWQL4ISyJ11L+0c8+OBWs94dxx47aXY/ff3xOBrL0VleGjczz1crsS6zTtrcZgAqeE35QC5sm9wKHp29495n7dd+Ik4uOvV6H3mm9H7zDcjii3R2rkwpiz4pZi+9MNRnDxlzMf37OqPvsEhv6YRyA3XUAK58sJre+I3/+yfx3Xfwe3dMbD1B3Fge3cM9GyK8v69UZoxJzo//qdRap8x5mP/9j/9cpz9jrHvA5AV/vkM5MqBUT4maCStnQujtfPQh5VXhofip/90X/Ru/JvYu+GxePvyT9XsdQDSzjWUQK5Mbjm1H3uFUku87d/+x4iIOPCTrXV7HYA08hMPyJX5HVNj9Pd3HzK0b/eIXz+4a1tERLRM6xjz8YW3XgcgL0zeQK5MbW2JrpntsXWMN+bsfPi2KE3viPYFy2JSx5lRqZTj4I6XY+/Gv47C5Ckx/RcvGfM1ujravSEHyBU/8YDcWX7W7Hhgw9ZRPzpoxgeujP7u9bF341/HcN9PozJ0MErT3h5t886NGe+/IibNmjvqc5eKhVi+aHa9Dh0gkbzLG8idk/2mnGr5TTlA3riGEsidhXOmxwULZo3523JORalYiAsWzBKTQO44Qwnk0rbd/XHRqidjsJYf7zN8MPb/1cr4jQ8ti/e9731xzjnnxDnnnBNvf/vba/caAAkkKIHcWrOxJ259fFPNnq9z2z/G+tV/GhERhUIhDv94Pf300+OP//iP47rrrqvZawEkickbyK2rzuuKz1+8qCbP9YWLz4rH/+TzMWnSpIiIOPrf6q+//nqUyz7oHMguQQnk2s3LF8aXLlsSrS3FCV9TWSoWorWlGF++bEn87vIF0dnZGTfddFOUSqUj9ykUCvHRj340fuu3fqvWhw6QGCZvgDh0TeXKtZviqS1vRKlYGPUjhSLiyO0XLJgVd1y6JObObD9y2/bt22P+/Plx4MCBKBaLUS6XY/HixfH444/HwoULG/GtADScM5QAETF3Zns8cP2y+PvPfTCuWzYv5nW0n/AbdQoRMa+jPa5bNi++c8sH44Hrlx0TkxERnZ2dceONN0ZERKlUijVr1sTg4GC8733vizVr1jTmmwFoMGcoAUbRNzgU/+G6FVEplOK/rfqvMb9j6rh+A8727dtj6dKlcfvtt8f1118fvb29ccMNN8RDDz0UN954Y6xatSra2toa8B0ANIagBBjDJZcc+jWLTzzxxIQeVy6Xo1j82QhUqVTi3nvvjc985jPx7ne/Ox599FETOJAZJm+AOjg6JiMOvTnn05/+dGzYsCH6+/tN4ECmCEqABjr33HPjmWeeiY985CNx9dVXx0033RQDAwPNPiyAqghKgAabPn16rF69Ou65556477774vzzz4/u7u5mHxbAKROUAE1gAgeyRFACNJEJHMgCQQnQZCZwIO0EJUACmMCBNBOUAAliAgfSSFACJIwJHEgbQQmQQCZwIE0EJUCCmcCBNBCUAAlnAgeSTlACpIAJHEgyQQmQIiZwIIkEJUDKmMCBpBGUAClkAgeSRFACpJgJHEgCQQmQciZwoNkEJUAGmMCBZhKUABliAgeaQVACZIwJHGg0QQmQQSZwoJEEJUCGmcCBRhCUABlnAgfqTVAC5IAJHKgnQQmQIyZwoB4EJUDOmMCBWhOUADlkAgdqSVAC5JgJHKgFQQmQcyZwoFqCEgATOFAVQQnAESZw4FQISgCOYQIHJkpQAnACEzgwEYISgFGZwIHxEJQAjMkEDpyMoATgpEzgwFgEJQDjZgIHRiIoAZiQkSbwzZs3N/uwgCYSlABM2PET+NKlS03gkGOCEoBTZgIHIgQlAFUygQOCEoCqmcAh3wQlADVjAod8EpQA1JQJHPJHUAJQcyZwyBdBCUDdmMAhHwQlAHVlAofsE5QA1J0JHLJNUALQMCZwyCZBCUBDmcAhewQlAA1nAodsEZQANI0JHLJBUALQVCZwSD9BCUDTmcAh3QQlAIlhAod0EpQAJIoJHNJHUAKQOCZwSBdBCUBimcAhHQQlAIlmAofkE5QAJJ4JHJJNUAKQGiZwSCZBCUCqmMAheQQlAKljAodkEZQApJYJHJJBUAKQaiZwaD5BCUDqmcChuQQlAJlhAofmEJQAZIoJHBpPUAKQOSZwaCxBCUBmmcChMQQlAJlmAof6E5QAZJ4JHOpLUAKQGyZwqA9BCUCumMCh9gQlALljAofaEpQA5JYJHGpDUAKQayZwqJ6gBCD3TOBQHUEJAG8xgcOpEZQAcBQTOEycoASA45jAYWIEJQCMwgQO4yMoAWAMJnA4OUEJACdhAoexCUoAGKfjJ/Abb7wx9u/f3+zDgqYTlAAwAUdP4Pfff3+8//3vN4GTe4ISACbIBA7HEpQAcIpM4HCIoASAKpjAQVACQNVM4OSdoASAGjGBk1eCEgBqyAROHglKAKgxEzh5IygBoE5M4OSFoASAOjKBkweCEgDqzARO1glKAGgQEzhZJSgBoIFM4GSRoASABjOBkzWCEgCaxAROVghKAGgiEzhZICgBoMlM4KSdoASAhDCBk1aCEgASxAROGglKAEgYEzhpIygBIKFM4KSFoASABDOBkwaCEgASzgRO0glKAEgJEzhJJSgBIEVM4CSRoASAlDGBkzSCEgBSygROUghKAEgxEzhJICgBIOVM4DSboASAjDCB0yyCEgAyxAROMwhKAMgYEziNJigBIKNM4DSKoASADDOB0wiCEgAyzgROvQlKAMgJEzj1IigBIEdM4NSDoASAnDGBU2uCEgByygROrQhKAMgxEzi1ICgBIOdM4FRLUAIAEWEC59QJSgDgCBM4p0JQAgDHMIEzUYISABiRCZzxEpQAwKhM4IyHoAQAxmQC52QEJQAwLiZwRiMoAYBxM4EzEkEJAEyICZzjCUoA4JSYwDlMUAIAp8wEToSgBACqZAJHUAIANWECzy9BCQDUjAk8nwQlAFBTJvD8EZQAQF2YwPNDUAIAdWMCzwdBCQDUlQk8+wQlANAQJvDsEpQAQMOYwLNJUAIADWUCzx5BCQA0hQk8OwQlANA0JvBsEJQAQFOZwNNPUAIAiWACTy9BCQAkhgk8nQQlAJAoJvD0EZQAQCKZwNNDUAIAiWUCTwdBCQAkmgk8+QQlAJAKJvDkEpQAQGqYwJNJUAIAqWICTx5BCQCkkgk8OQQlAJBaJvBkEJQAQKqZwJtPUAIAmWACbx5BCQBkhgm8OQQlAJApJvDGE5QAQCaZwBtHUAIAmWUCbwxBCQBkmgm8/gQlAJALJvD6aWn2ATRb3+BQvLKrLw4MlWNySzHmd0yNqa25/88CAJl0eAK/8MIL4zOf+UysX78+HnnkkVi0aNGEnkc/HCuX33n3jt5YvaEn1r20M3p290flqNsKEdE1sz2WnzU7rlnWFQvnTG/WYQIAdXB4Al+2bFlcfvnlsXTp0rj33nvjqquuGvNx+mF0hUqlUjn53bJh2+7+WLl2Uzy15Y0oFQsxXB79Wz98+wULZsUdly6JuTPbG3ikQFJccsklERHxxBNPNPlIgHro7e2NG264IR566KG44YYbYtWqVTFlypRj7qMfTi4311Cu2dgTF616Mp5+eVdExJh/GY6+/emXd8VFq56MNRt76n6MAEBjnexd4PphfHIRlF9b1x23Pr4pBofKJ/2LcLzhciUGh8px6+Ob4mvruut0hABAs4z2LnD9MH6ZD8o1G3viK9+uzedNfeXbm+PhnPxLAwDy5uh3ga+48z79MAGZvoZy2+7+uGjVkzE4VK7Zc7a2FOM7t3woN9dEQN65hhLyp2dXX/zKn/5THDxJPhzY+aPYu/FvYqBnUwzv2x2FYikmzTwj2t9zQUw799eiNOVnb8zJej9k+l3eK9duiqERTlHve+47set/3vWzLxSKUWqfEa1di+NtF1wbk2aeMepzDpUrsXLtpnjg+mV1OGIAoNn+y18/H+UoRMTo59x6v/93sfvbX49JM8+IGcsui0kdc6NSHo4Dr3fHvv/zv2Lwxy/G7I/+wZH7Z70fMhuU3Tt646ktb4x5n47f+FxM6jgzKkMHYvDHP4w9Tz8SA1s3xTs+fXeU2qaN+JjhciWe2vJGbNnZGwtm5+sjAQAg68bTD4M//mHs/t9/Hm3v/IWYfdkfRKFl0pHbprzzF+K0X7o09r/87DGPyXo/ZPYaytUbeqJULIx5n0k/Ny9az3h3tM07J2Z84Mo47fyPRbn/zdi/+XtjPq5ULMSD60+8FmLfvn3x4IMPxt69e6s6dgCgOcbTD3uefiSiUIiOX7/5mJg8rFCaFO0LTzwTOVo/ZEFmg3LdSzsn/I6s1s4FEREx3PfmmPcbLldi3eadR/68b9+++PKXvxxnnnlmXHfdda61AoCUOlk/VMrDMdDzXEw+fUG0nPZzE3ru4/shSzI5ee8bHIqe3f0TftzQmzsiIsa8hvKwnl39sXP3nrjv3rvjzjvvjL1798bh9zcNDw9P+LUBgOYaTz+U9++NysHBaJkx55Reo2dXf/QNDmXu1zRm67t5y9ZdfWNcRnuUSjkq5eFD11C++n9jz9MPR+vcxTFlhNPUJzw0IuYvPi/2bz/xs6XuuuuueOyxxyZ83EDy/Ou//mtE/Ozd3kB27W/riMq7PlbX16hExCu7+uLsd8yo6+s0WiaD8sA4Pybo9f/xe8f8eVLH3Jj90T+IQrE0vhcqZfI/HwDkUqVw8v//F6ecFoVJrTG0Z8cpv854OyVNMllEk1vGd2lox4f/86G3+R/YH30//G7s+/7fxU+e+JOYc8UXx/X47677x/iHx/77CZP35z73ufj4xz9+yscPJIfPoYT8eOG1PfGbf/bPY96nUCxF27xzY//Lz8TQ3jei5bRZE36d8XZKmmTvO4qI+R1TY+z3Zx0yqWNutHYujLZ550THr98c0869OAZefib6Xhz7L1NERCEi3jN3Vvz+7/9+vPrqq3HnnXfGjBmHTl+XSuM8wwkAJMZ4+2HG+y+PqFRi19/9WVSGD55we2V4KPq7N4z42MJbr5M1mQzKqa0t0XUKn0T/tuWfimLbtNjz1OqoVMY+Hd3V0X7kgtpp06YdCcsHH3zQtVYAkELj7YfWM94TM3/td2Lgle/H9vs/F73P/m0M9GyK/a98P/ZseCxe+8ZNse+5vx/xsUf3Q5Zk7zt6y/KzZscDG7ZO6KODSm3T4rT3Xx5vrrsv+l54MqYtXj7y/YqFWL5o9glfnzZtWlxzzTWnfMwAQHONtx+m/5tfj9bORbF349/EnvV/FcN9P41CsSUmzTwjpr73wpi+9MMnPGa0fsiCzAblNcu64v7vvTLhx5229CPR+8y3Ys+/PBRT3/vBEd+gM1yuxLXnd9XgKAGAJJlIP0ye8/Mx68O3jPu5s9wPmQ3KhXOmxwULZsXTL+864V8Z0865KKadc9GIjyu0TI4zf+e+UZ+3VCzEB36+I5O/NgkA8m6sfqhG1vshk9dQHnbHpUui5SS/PmkiKpVKVIYOxsqL31mz5wQAkqXW/RAR0VIsxB2XLqnpcyZJpoNy7sz2+OIlZ9fs+QqFQuxd95fx73/ll+MHP/hBzZ4XAEiOWvdDRMQfXXJ2zD2FNwynRaaDMiLiqvO64vMXL6rJc33h4rPi6Qe/Eu3t7bFs2bL4i7/4iyOfPQkAZEet++HK87J57eRhmQ/KiIibly+ML122JFpbilGa4CnsUrEQrS3F+PJlS+J3ly+IRYsWxfr16+OTn/xk3HDDDXHNNddEb29vnY4cAGiWWvZD1hUqOTrFtm13f6xcuyme2vJGlIqFMS+2PXz7BQtmxR2XLhnxNPWaNWtixYoV0dnZGY8++mice+659Tx8oAn8phyg1v2QRbkKysO6d/TG6g09sW7zzujZ1R9H/wcoxKEPHV2+aHZce37XSd+NtXnz5rjiiivixRdfjK9+9auxYsWKKBRqeyEv0DyCEjislv2QNbkMyqP1DQ7FK7v64sBQOSa3FGN+x9QJf4L9wMBA3HLLLXH33XfH1VdfHffcc09Mn56vv0iQVYISGEkt+iFLch+UtWQCh+wRlAAnl4s35TTKVVddFc8884x3gQMAuSIoa8y7wAGAvBGUddDW1hZf//rX46GHHopvfvObsXTpUh+EDgBklqCsIxM4AJAHgrLOTOAAQNYJygYwgQMAWSYoG8gEDgBkkaBsMBM4AJA1grIJTOAAQJYIyiYygQMAWSAom8wEDgCknaBMABM4AJBmgjJBTOAAQBoJyoQxgQMAaSMoE8gEDgCkiaBMMBM4AJAGgjLhTOAAQNIJyhQwgQMASSYoU8QEDgAkkaBMGRM4AJA0gjKFTOAAQJIIyhQzgQMASSAoU84EDgA0m6DMABM4ANBMgjJDTOAAQDMIyowxgQMAjSYoM8gEDgA0kqDMMBM4ANAIgjLjTOAAQL0JyhwwgQMA9SQoc8QEDgDUg6DMGRM4AFBrgjKHTOAAQC0JyhwzgQMAtSAoc84EDgBUS1BiAgcAqiIoOcIEDgCcCkHJMUzgAMBECUpOYAIHACZCUDIqEzgAMB6CkjGZwAGAkxGUnJQJHAAYi6Bk3EzgAMBIBCUTYgIHAI4nKJkwEzgAcDRBySkzgQMAEYKSKpnAAQBBSdVM4ACQb4KSmjGBA0A+CUpqygQOAPkjKKk5EzgA5IugpG5M4ACQD4KSujKBA0D2CUrqzgQOANkmKGkYEzgAZJOgpKFM4ACQPYKShjOBA0C2CEqaxgQOANkgKGkqEzgApJ+gpOlM4ACQboKSxDCBA0A6CUoSxQQOAOkjKEkcEzgApIugJLFM4ACQDoKSRDOBA0DyCUoSzwQOAMkmKEkNEzgAJJOgJFVM4ACQPIKS1DGBA0CyCEpSywQOAMkgKEk1EzgANJ+gJPVM4ADQXIKSzDCBA0BzCEoyxQQOAI0nKMkcEzgANJagJLNM4ADQGIKSTDOBA0D9CUoyzwQOAPUlKMkNEzgA1IegJFdM4ABQe4KS3DGBA0BtCUpyywQOALUhKMk1EzgAVE9QknsmcACojqCEt5jAAeDUCEo4igkcACZOUMJxTOAAMDGCEkZhAgeA8RGUMAYTOACcnKCEkzCBA8DYBCWMkwkcAEYmKGECTOAAcCJBCRNkAgeAYwlKOEUmcAA4RFBCFUzgACAooWomcADyTlBCjZjAAcgrQQk1ZAIHII8EJdSYCRyAvBGUUCcmcADyQlBCHZnAAcgDQQl1ZgIHIOsEJTSICRyArBKU0EAmcACySFBCg5nAAcgaQQlNYgIHICsEJTSRCRyALBCU0GQmcADSTlBCQpjAAUgrQQkJYgIHII0EJSSMCRyAtBGUkFAmcADSQlBCgpnAAUgDQQkJZwIHIOkEJaSECRyApBKUkCImcACSSFBCypjAAUgaQQkpZQIHICkEJaSYCRyAJBCUkHImcACaTVBCRpjAAWgWQQkZYgIHoBkEJWSMCRyARhOUkFEmcAAaRVBChpnAAWgEQQkZZwIHoN4EJeSECRyAehGUkCMmcADqQVBCzpjAAag1QQk5ZQIHoFYEJeSYCRyAWhCUkHMmcACqJSiBiDg0gT/77LMmcAAmTFACRyxcuNAEDsCECUrgGCZwACZKUAIjMoEDMF6CEhiVCRyA8RCUwJhM4ACcjKAExsUEDsBoBCUwbiZwAEYiKIEJMYEDcDxBCZwSEzgAhwlK4JSZwAGIEJRAlUzgAAhKoCZM4AD5JSiBmjGBA+SToARqygQOkD+CEqgLEzhAfghKoG5M4AD5ICiBujKBA2SfoAQawgQOkF2CEmgYEzhANglKoKFM4ADZIyiBpjCBA2SHoASaxgQOkA2CEmgqEzhA+glKIBFM4ADpJSiBxDCBA6SToAQSxQQOkD6CEkgkEzhAeghKILFM4ADpICiBRDOBAySfoARSwQQOkFyFip/IQIoMDAzELbfcEnfffXdcffXVcc8998T06dNr+hrbt2+P+++/P4aHh+PBBx+MiIhrr702SqVSfOITn4jOzs6avh5A2glKIJXWrFkTK1asiM7Oznj00Ufj3HPPrdlzP/LII3HllVdGqVSKcrkcERHFYjGGh4fj4YcfjiuuuKJmrwWQBYISSK3u7u64/PLL48UXX4yvfvWrsWLFiigUClU/7+DgYMyfPz9ef/31Y77e2dkZP/rRj6K1tbXq1wDIEtdQAqlVr3eBt7a2xm233XbC12+77TYxCTACZyiBTKj1BH78WUpnJwFG5wwlkAnjeRf44eshx+P4s5TOTgKMTlACmTHWBP6Nb3wjzjzzzNi+ffu4n+9Tn/pUtLa2Rmtra3zyk5+s12EDpJ7JG8ikoyfw22+/Pa677ro4ePBgfPazn4277rprXM/RNzgU9z/2tzFUjrjo310Y8zumxtTWlvoeOEAKCUogs7q7u+Oyyy6L559/PorFYpTL5Zg8eXK88soro36WZPeO3li9oSfWvbQzenb3x9E/IAsR0TWzPZafNTuuWdYVC+fU9vMvAdJKUAKZValU4mMf+1isXbv2yPWUpVIpbr755hPOUm7b3R8r126Kp7a8EaViIYbLo/9oPHz7BQtmxR2XLom5M9vr+W0AJJ6gBDLr/vvvH/Hax0mTJsXWrVuPnKVcs7En/vCJF2KoXBkzJI9XKhaipViIL15ydlx1XlfNjhsgbbwpB8isUql0zLR9+EPPDx48GNdff31ERHxtXXfc+vimGBwqTygmIyKGy5UYHCrHrY9viq+t667dgQOkjDOUQOa9+eab8dxzz8Vzzz0Xzz77bHzrW9+KefPmxe/9+WNx6+ObavY6X75sSVzpTCWQQ4ISyKVtu/vjolVPxuDQ+D+b8mRaW4rxnVs+5JpKIHdM3kAurVy7KYbemrj7Xvzn2PqlD0ffD797wv1e+8ubY+uXPhz7X37mhNt+fPdvx/b7Pnvkz0PlSqxcW7szngBpISiB3One0RtPbXnjyDWTbV1LIqIQA1ufO+Z+w/t74+BPtkZhUlsM9BwbikN734ihN1+P1q4lP7t/uRJPbXkjtuys/veJA6SJoARyZ/WGnigVC0f+XGqfEZN+bl4M9Dx/zP0GezZFFEsx7ZxfPSE2B3oO/blt3jnHfL1ULMSD63vqdOQAySQogdxZ99LOE97R3da1JIZ2vxpD+3Yf+dpAz6Zo7VwYU971i3Hg9S1RHuw/5rYoFKPtzLOPeZ7hciXWbd5Z328AIGEEJZAr+waHomd3/wlfP3ymcfCoaXugZ1O0di2O1jPfG1EoxOCrLxy5bXDrczH59HdFsW3qCc/Vs6s/+gaH6nD0AMkkKIFc2bqrL0b6aIvWriURheKRayWH9++Ngz/ZGm1zF0dx8pSYPOddMbD10G1De38SQ3t2RFvXOSM8U0QlIl7Z1Ven7wAgeQQlkCsHRvmYoFLbtJg8+51HgnKw5/mIYvHQ2cmIaOtafOS6ycPXUx5//eR4XgcgiwQlkCuTW0b/sdfatSSGdv84hnp3xcDW52Ly6QuiOHnKodvmLokDO16O8kDfoegslo7E5kRfByBr/MQDcmV+x9QojHLb0ddRDmzbFG1zF//strficWDb80ferHM4No9XeOt1APJCUAK5MrW1JbpG+U02bXMXRxSK0ffSv8TBn/S89fmUhxTbpsbk2e+Mvuf/IYb37IjWUa6fjIjo6miPqa0tNT92gKQSlEDuLD9r9jGfQ3lYsbU9Jp/+rti/eX1EoXDCpN3atTj6X/peRIx+/WSpWIjli2bX/qABEkxQArlzzbKuEz6H8rBD79yuxOQ574pi67FnMtvmLomISkSpJVrPePeIjx8uV+La87tqfMQAyVaoVCoj/1QFyLDr/nJDPP3yrlHD8lSUioX4wM93xAPXL6vZcwKkgTOUQC7dcemSaBlh9q5GS7EQd1y65OR3BMgYQQnk0tyZ7fHFS84++R0n4I8uOTvmjvKGH4AsE5RAbl11Xld8/uJFNXmuL1x8Vlx5nmsngXxyDSWQe2s29sQfPvFCDJUrE7qmslQsREuxEH90ydliEsg1QQkQEdt298fKtZviqS1vRKlYGDMsD99+wYJZccelS8zcQO4JSoCjdO/ojdUbemLd5p3Rs6s/jv4BWYhDH1q+fNHsuPb8rlgwe3qzDhMgUQQlwCj6BofilV19cWCoHJNbijG/Y6rfgAMwAkEJAEBVvMsbAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICqCEoAAKoiKAEAqIqgBACgKoISAICq/H8M9EAlqgaqBgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#define model\n", + "model = BayesianModel([('C','S'),('C','R'),('S','W'),('R','W')])\n", + "\n", + "\n", + "# Need to calculate position and pass it as a parameter to the nx.draw method\n", + "# The below line is referred from the stackoverflow\n", + "# https://stackoverflow.com/questions/71607514/stopiteration-error-while-drawing-a-pgmpy-networkx-graph\n", + "pos = nx.circular_layout(model)\n", + "\n", + "#Pass pos as pos argument\n", + "nx.draw(model, pos=pos, with_labels=True)\n", + "\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Ho61X-iE0fET" + }, + "source": [ + "Next, we define the conditional probability tables for all nodes, in this case: C,S,R,W\n", + "\n", + "We then add them to the Bayesian network and cross-check that the netowrk is well defined." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 858, + "status": "ok", + "timestamp": 1611834676467, + "user": { + "displayName": "Ulrich Kerzel", + "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64", + "userId": "10582137377195463323" + }, + "user_tz": -60 + }, + "id": "6W0wgisEyvaS", + "outputId": "2a247281-010d-4bd6-8a6f-4efd25313917" + }, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#define CPT\n", + "\n", + "#cloudy\n", + "CPT_C = TabularCPD(variable='C', variable_card=2, \n", + " values=[[0.5], [0.5]])\n", + "\n", + "#sprinkler\n", + "CPT_S = TabularCPD(variable='S', variable_card=2,\n", + " values=[ [0.5, 0.9],\n", + " [0.5, 0.1]\n", + " ],\n", + " evidence=['C'], evidence_card=[2]\n", + " )\n", + "\n", + "#rain\n", + "CPT_R = TabularCPD(variable='R', variable_card=2,\n", + " values=[ [0.8, 0.2],\n", + " [0.2, 0.8]\n", + " ],\n", + " evidence=['C'], evidence_card=[2]\n", + " )\n", + "\n", + "# wet grass\n", + "CPT_W = TabularCPD(variable='W', variable_card=2,\n", + " values=[ [1.0, 0.1, 0.1, 0.001],\n", + " [0.0, 0.9, 0.9, 0.99]\n", + " ],\n", + " evidence=['S','R'], evidence_card=[2,2]\n", + " )\n", + "model.add_cpds(CPT_C, CPT_S, CPT_R, CPT_W)\n", + "model.check_model()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Q4QIQeXO02BA" + }, + "source": [ + "Once the Bayesian network is fully defined we can use it for inference.\n", + "\n", + "In the example below, we with to estimate the probability for rain (R) given that we observe that the grass is wet (W=1).\n", + "\n", + "Exerise: Experiment and infer the values of different quantities given some observations." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 693, + "status": "ok", + "timestamp": 1611834713348, + "user": { + "displayName": "Ulrich Kerzel", + "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64", + "userId": "10582137377195463323" + }, + "user_tz": -60 + }, + "id": "HTeS6XwEyvaV", + "outputId": "d757422b-0004-4540-e866-3db0bc59fe13" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "+------+----------+\n", + "| R | phi(R) |\n", + "+======+==========+\n", + "| R(0) | 0.2921 |\n", + "+------+----------+\n", + "| R(1) | 0.7079 |\n", + "+------+----------+\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/kathi/opt/anaconda3/lib/python3.9/site-packages/pgmpy/models/BayesianModel.py:8: FutureWarning: BayesianModel has been renamed to BayesianNetwork. Please use BayesianNetwork class, BayesianModel will be removed in future.\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "# Infering the posterior probability \n", + "from pgmpy.inference import VariableElimination\n", + "\n", + "infer = VariableElimination(model)\n", + "posterior_p = infer.query(variables=['R'], evidence={'W': 1})\n", + "print(posterior_p)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "executionInfo": { + "elapsed": 675, + "status": "ok", + "timestamp": 1611834830851, + "user": { + "displayName": "Ulrich Kerzel", + "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14GgUe7aETi8_2hp0GY6lXuHo8CbZgw1Wjq6SJtE-ow=s64", + "userId": "10582137377195463323" + }, + "user_tz": -60 + }, + "id": "SlMEg9fAyvaZ", + "outputId": "c3a7a572-c3ad-402d-ebb1-ea9c5a94d8c6" + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(C ⟂ W | R, S)\n", + "(R ⟂ S | C)\n", + "(S ⟂ R | C)\n", + "(W ⟂ C | R, S)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# independence conditions\n", + "model.get_independencies()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "colab": { + "collapsed_sections": [], + "name": "Sprinkler.ipynb", + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.9.13" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}