{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Logistic regression" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Logistic regression can be used to learn the probability of an event being true or false as a function of one or more features.\n", "\n", "Here we present as example a simple 'linear' logistic regression." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We import the following packages, classes, and functions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# for handling data:\n", "import numpy as np\n", "import pandas as pd\n", "\n", "# for plotting:\n", "import matplotlib.pyplot as plt\n", "\n", "import halerium.core as hal\n", "\n", "# for graphs:\n", "from halerium.core import Graph, Entity, Variable, StaticVariable\n", "from halerium.core.regression import linear_regression, polynomial_regression, connect_via_regression\n", "from halerium.core.distribution import BernoulliDistribution\n", "\n", "# for models:\n", "from halerium.core import DataLinker, get_data_linker\n", "from halerium.core.model import MAPModel, ForwardModel, Trainer\n", "from halerium.core.model import get_posterior_model\n", "\n", "# for predictions:\n", "from halerium import Predictor\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create example data, we simply build a forward model of logistic regression:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "n_data = 100\n", "\n", "x_scatter = 10\n", "\n", "with Graph(\"graph\") as graph:\n", "\n", " x = Variable(\"x\", shape=(2,), mean=0, variance=x_scatter**2)\n", " y = Variable(\"y\", shape=(), distribution=BernoulliDistribution)\n", "\n", " connect_via_regression(\n", " name_prefix=\"parameters\",\n", " inputs=x,\n", " outputs=y,\n", " order=1,\n", " )\n", "\n", "slope = graph.parameters_y.location.slope\n", "intercept = graph.parameters_y.location.intercept\n", "\n", "model = ForwardModel(graph, data=DataLinker(n_data))\n", "x_data, y_data, slope_data, intercept_data = model.get_example((x, y, slope, intercept))\n", "\n", "data = pd.DataFrame()\n", "data[\"x_1\"] = x_data[:,0]\n", "data[\"x_2\"] = x_data[:,1]\n", "data[\"y\"] = y_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's visualize the generated data:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEICAYAAABMGMOEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt+ElEQVR4nO3de3hU1bn48e8iQLirQEQEAyogNwExCvkFCnLJREVQqIVWzwPVHqCU9tRbQasNcCxgRVs9hVYrtSg8KgdBeaxM5BYFCVgi5CBE5H6TAiGi4U6S9/fHTGICIdmZ277k/TzPPDOZPbNnrRnY717vWnstIyIopZRSValldwGUUkq5gwYMpZRSlmjAUEopZYkGDKWUUpZowFBKKWWJBgyllFKW1La7AJHUvHlzadu2rd3FUEopV8nOzs4TkYSqXuepgNG2bVs2btxodzGUUspVjDH7rLxOU1JKKaUs0YChlFLKEg0YSimlLNGAoZRSyhINGEoppSzRgBFFWVkwY0bgXiml3M5Tw2qdJCsLBg6E8+ehbl1YuRKSk+0ulVJKhU5bGFGSmRkIFkVFgfvMTLtLpJRS4dGAESX9+wdaFnFxgfv+/e0ukVJKhcf2lJQxph7wCRBPoDyLRCTdGNMUeAdoC+wFfiQi39hVzupKTg6koTIzA8FC01FKKbezPWAA54ABInLSGFMHWGuMWQYMB1aKyExjzGRgMjDJzoJWV3KyBgqllHfYnpKSgJPBP+sEbwIMA+YFn58H3Bv70imllCphe8AAMMbEGWM2A0eB5SKyAWghIocBgvdX21hEpZSq8RwRMESkSER6AK2B240xXa2+1xgz1hiz0Riz8dixY1Ero1JK1XSOCBglROQEkAmkAUeMMS0BgvdHL/OeV0UkSUSSEhKqnM5dKaVUiGwPGMaYBGPMlcHH9YFBwJfAUmB08GWjgfdtKaBSSinAGaOkWgLzjDFxBALYQhH5wBiTBSw0xjwM7Afut7OQSilV09keMETk/4BbKnj+ODAw9iVSSilVEdtTUkoppdxBA4ZSSilLNGAopZSyRAOGUkopSzRgKKWUskQDhlJKKUs0YCillLJEA4ZSSilLNGAopZSyRAOGUkopSzRgKKWUskQDhlJKKUs0YCillLJEA4ZSSilLNGAopZSyRAOGUkopSzRgKKWUskQDhlJKKUs0YCillLJEA4ZSSilLNGAopZSyRAOGUkopSzRgKKWUskQDhlJKKUs0YCilVARlZcGMGYF7r6ltdwGUUsorsrJg4EA4fx7q1oWVKyE52e5SRY7tLQxjzHXGmNXGmFxjzFZjzH8Fn29qjFlujNkRvL/K7rIqpVRlMjMDwaKoKHCfmWl3iSLL9oABFAKPiUgnoDfwC2NMZ2AysFJE2gMrg38rpWzi5VRLpPTvH2hZxMUF7vv3t7tEkWV7SkpEDgOHg48LjDG5QCtgGNA/+LJ5QCYwyYYiKlXjeT3VEinJyYHvJjMzECy89h3ZHjDKMsa0BW4BNgAtgsEEETlsjLn6Mu8ZC4wFSExMjFFJlapZKkq1eO1gGCnJyd79bpyQkgLAGNMIeBf4tYh8Z/V9IvKqiCSJSFJCQkL0CqhUNXkpheP1VIuyxhEtDGNMHQLBYoGILA4+fcQY0zLYumgJHLWvhEpVj9dSOF5PtShrbA8YxhgDzAVyReTFMpuWAqOBmcH7920onlIh8WIKx8upFmWN7QEDSAH+A9hijNkcfO4pAoFioTHmYWA/cL89xVOq+kpSOCUtDE3hKC+wPWCIyFrAXGbzwFiWRalIiVUKJytL00QqdmwPGEp5VbRTOF7rJ1HO55hRUkqp6vH6VcXKeTRgKOVSbhvq6qVhxjWVpqSUcik3DXXV9Jk3aMBQysXcMtTVi8OMayJNSSmlos5t6TNVMW1hKGWTmjQk1k3pM3V5GjBUSGrSwS4aamJO3y3pM3V5GjBUtdXEg12kaU5fuZH2Yahq0/H/4avpOX0dYutO2sJQ1abzJIWvJuf0tYXqXhowVLXV5INdJNXUnL6m49xLA4YKSU092KnwaQvVvTRgKKViSluo7qUBQykVc9pCdScdJaWUUsoSDRg1gBeHMHqxTm6jv0HNoykpj/PiEEYv1slt9DeombSF4XFevMjOi3VyG/0NaiYNGB7nxSuKvVgnt9HfoGbSlJTHeXEIoxfr5Db6G9RMRkTsLkPEJCUlycaNG+0uhnI5nYlX1TTGmGwRSarqddrCUID7DpLRKq925ip1eRowlOsOktEsr85zpNzi+PHjrFixAr/fz+OPP06XLl2i/pna6a1cN+KlbHnPnoU33ojcviPVmavXKKhoyM/PJz09nd69e5OQkMCoUaN4//332blzZ0w+3/Y+DGPM34EhwFER6Rp8rinwDtAW2Av8SES+qWpf2ocRGje2MO64A86dC/xdt25kWwLhprvc9n0q5zp48CAZGRlcddVVDB8+nIKCAq6++mp69OhBWloaaWlpJCUlERcXF9bnuKkP4x/An4Gy54mTgZUiMtMYMzn49yQbylYjuG3ES3Iy/PSn8MorIBJoaUQyYIQ7z1G00lpu62eyyqv1CtXq1av54IMPyMjIYOvWrQAMGzaM4cOH07hxY44dO0ajRo3sKZyI2H4j0JL4oszf24GWwcctge1W9nPrrbeKqhnWrROpX18kLi5wv26d3SX6XjTK5uT6hsOr9bKquLhYtm/fLgsWLCh9zufzSd26dWXQoEHy/PPPy5YtW6S4uDiq5QA2ioVjrBNaGBVpISKHAUTksDHm6su90BgzFhgLkJiYGKPiqVi7+CzUya2iysoW6tm0VzvjvVqvynz33XesWrWKjIwM/H4/e/fuxRiDz+ejWbNmvPLKKzRv3pyGDRvaXdRLWYkq0b5xaQvjxEXbv7GyH21heJNXzkLDqYdXvoOLebVeZRUVFUl2drbk5+eLiMjs2bMFkEaNGsmwYcPkL3/5i+zevdvWMuLyFsYRY0xLCbQuWgJH7S6Qso9XzkLDqYeTW1Th8Gq9jh49yvLly/H7/Xz00UccPXqU119/nTFjxjBixAi6dOlCcnIydevWtbuo1eLUgLEUGA3MDN6/b29xlJ28sqRnuPXw6qJDXqjXhQsXyM/Pp0WLFhw5coSWLVsiIjRv3pzU1FR8Ph933nknAC1atKBFixY2lzg0ThhW+xbQH2gOHAHSgfeAhUAisB+4X0Tyq9qXDqv1Lq+MpPFKPRTs3bu3tB9i5cqVDBo0iMWLFwMwe/ZsevXqRc+ePalVy/mXu1kdVmt7wIgkDRhKeZfdwbawsJDatQNJmREjRpQGh8TERNLS0hg2bBh33XVX7AsWAW66DkMppSplx8WQIsK2bdtKWxHZ2dl8/fXXxMfHM3ToUPr27UtaWho33XQTxpjoFsYhNGAopRwv1gMf3n//fSZOnMjBgwcB6Ny5M2PGjOHUqVPEx8czevTo6H24gzk/uaaUSzh5/ignl82KaC3YVFRUxGeffca0adNISUnB7/cD0KpVK3r16sXf/vY39u/fz9atW3nhhRdo2rRpZD7YpbSFoVQEOHn+KCeXzapID789ceIE48ePZ/ny5eTn52OMISkpieLiYgCSkpJYtGhR2OX2Gg0YSkWA1ZRJNDpuq9qnV65jCXX47blz5/j0009LJ/GbPHkyTZo0ITc3lyFDhuDz+UhNTaV58+aRL7THaMBQNVqkDuBWrrGIxpm+lX165TqW6lqwYAFvv/02q1ev5tSpU9SpU4eRI0cCUKtWLXJycmwuoftowFA1ViQP4FZSJtE407eyT6dcTR3NYbEnT55k9erVrF27lpkzZ2KMYdWqVWzbto3Ro0fj8/m44447aNy4cWQ/uIbRgKEcJZZj7SN9AK8qZRKNM32r+7T7aupotK727dvHO++8g9/vZ+3atVy4cIEGDRowYcIE2rRpw+zZs6lXr15kKqAADRjKQWLdORvrVE00zvSd0nqoSiSC8/Hjx1m+fDlJSUm0a9eOzZs3M2nSJLp168YjjzyCz+cjJSWF+Ph4AA0WUaABQ8WM0zpn7TjYRuNM3+7WgxWhBOeioiI2bNhQeuHcv/71L0SEmTNnMmnSJFJTUzl06BDXXntttIuvgnRqEBUTVloPXhj+qS7PSrrx4MGDHD16lJ49e3L69GmaNm3KhQsX6NWrF2lpafh8vogsSarK06lBlKO4qXNWRUdFLaGzZ8/yySeflLYitm3bxu23386GDRto0KABGRkZdOvWjauuusqeQqtyNGComHBL56wqr7JWQSgDFESEPXv2cMMNNwDwk5/8hCVLllC3bl369evHQw89RFpaWunr+/XrF5GyqsjQlJSKGf0P7S6VpQirkz789ttvWbVqFX6/n4yMDPbt28fhw4e55pprWLNmDQUFBfTr1y+sJUk1nRkeTUkpx9HWQ+WcFlArSyNWtq24uJjCwkLq1q3L+++/z4gRIygqKqJx48YMHDiQyZMnU79+fQD69u0b9bLaYsqUwM1jdPJBpRyg5Az5mWcC906YJLCyCf8u3tat2xHefPNNHnzwQa655hrmzZsHwK233spvfvMbPv74Y44fP86SJUsYP348V1xxRczKaoupU20uQHRYbmEYYwYDPwJmi8hmY8xYEXk1ekVTqmJ2nIlH+zMdd4ZM5YMQSrYtX36Wt97qw5Ah2QAkJCSQmppKhw4dAGjdujXTp0+3tawqcqqTkpoA/BR42hjTFOgRlRIpVYlQc9XhHPBjkR936nxPZdOIe/bsKR3N1LBhQxYsWEBycj327evOf/zHcHw+H7fccottS5LanvKcMqV8y6JkUaX0dM+kp6oTMI6JyAngcWPMTOC26BRJqcsL5Uw83AN+LM7+I3GGHK1W0EsvvcScOXP46quvAGjTpg0jRowo3T537tzIfZible23MAY8NKCoRHUCxj9LHojIZGPML6NQHqUqFcqZeLgH/Fid/YdzhhyJVpCIsHXrVjIyMli5ciWLFi2iQYMGnD59mhtvvJEJEyaQlpZGhw4dasySpOoiIlLpDfgTweG3Tr/deuutouy1bp3I9OmBe6d8xrp1IvXri8TFBe5DKVss6hWO6dMD9YPA/fTp1t5XXFwsX3zxhTzUo4e0atVKAAGkc+fO8uWXX1a/IOnp1X+PF7nsewA2ioVjbJXXYRhjngW6AyNF5LQxJhVIF5GUqEayEOh1GPZy8lh4pw1ZjbRXX4WJEwOtqPj4y3/3RUVFbNy4Eb/fz4ABA+jbty8bN25k0G23MWjEiNLpN6677rrQCuLRVIzXRew6DBF52hjzE+BjY8w54BQwOQJlVB7jxJE+JWzvEI2irCz49a8D33utWvCnP5Wva2FhIfPnz8fv95dbkjQ+Pp6+ffvSs2dP8oDauiSpqkKVwxmMMQOB/yQQKBKAX4nImmgXTLmP48bCVyIrC2bMqPp6B6uvs1NJoC4uDpzcHzlyjpUrV/K///u/AMTFxfHUU0/x8ccfc8899/DWW29x7NgxJp89C8ZQKy4ucOZoTOBW3RE9U6Z8/95w9qMcz0pKahXwOxFZa4y5GXgTeFREVsWigNWhKSn7WUn92J0espo6c3KKraysLLjjjt2cP/8hxvipW3c1Z8+e5oYbbmDXrl0ApdOAX7azOlKpJE1JuZLVlFSVLQwRGSAia4OPtwB3As+GX8SqGWPSjDHbjTE7jTGaBnOB5GR48snKg4XdVzRXlDoL53V2KCgo4IMPPqC4uJjkZLjzzj8g8ktatvyShx/+KUuXLi23ZnWrVq10ZJMKW7XnkhKRw8E0VVQZY+KA2cBg4CDwL2PMUhHZFu3PVtHjhH4Oq8NknXQxnYiQk5NDRkYGCxf6ycn5lKKiC2RnZ9OzZ09mzXqC559/nHbt2oX2AenpkSlopPajHCmkyQdF5EykC1KB24GdIrIbwBjzNjAM0IDhYk44CCcnBzqG330XRoy4fMCye7qJvLw8RISEhAT8fj933XUXAMZ0Ax6hbl0fBQVdALjxxhvD+7BI9Tdov4WnOXm22lbAgTJ/HwR62VQWFSF2H4Th+1FF58/DmjVw882VB41YlbGwsJANGzbg9/vx+/1kZ2fzu9/9jilTpvCDH/yA119/ne3bU3n++WspKgq00tatg2osGaFUWJwcMCpKuF7Sm2aMGQuMBUhMTIx2mVQE2D3E1QlpsRKnTp2iYcOGFBUV0aZNG77++mtq1apF7969mTp1Kvfddx8ADRs2ZMyYMWRlwUsvOSNNpmoeJweMg0DZq4daA19f/CIJzJj7KgRGScWmaMrNyqbF4uJg//5AqyMWQePMmTPlliStV68en3/+OXFxcfzmN7/h2muvZdCgQZddkrSkhVb49BRqPzvFkaO2lIdZuRzcjhuBYLYbuB6oC+QAXSp7j04Noqxat05k/HiR+PjwpgypSnFxcenjadOmSb169QSQ+Ph4GTx4sLz44ovlXmMZhFcwl01d4Wge+C6xODWIYxdQEpFCYCKQAeQCC0Vkq72lii03XDRmhRPrkZwMiYlQWBjisNlKOne//fZbFi9ezLhx42jbti0HDx4EoFOnTowbN44PP/yQ/Px8PvroIx555BF7hrt6dIEfW9Sk79JKVHHLzUstjEhMmOcETq5HWGWr4Aw/JydH+vTpI3FxcQJIkyZN5L777pPc3NzIFDg9PfC5F99COcMNt4WivueB7xK3tzBqOidfNFYdTq5HSX/Af/93FVdxV9CaOAK8+eabPPDAA8yfPx+A5s2bc/r0aSZNmsQnn3xCXl4eixcvpmPHjlGrQ7XoFB6RU0O/yyqnBnETL00N4pZpKariiXoEp7uQ9HSenjaNZcCm4KYE4Cmfj1/7/baUybb3q+954LuM2Gy1yh5OuF4hEtxcjz179uD3+8kDngHM1Kms+ugjGtety/RPPuHQL7L58Y97kJJSzYZ62ZXZlHITK3krt9y81IfhRU5fhEhEZN26dfLLX/5S2jdt+v1iQiBFwf6CwmeeCZQfQu+XiUTOO9yROR4Y2eMYHvgusdiHoS0MFRNOTE2JCF988QUZGRmMHz+eLVsa8cwzy1m79jUGDryDiT4faWlptL/pJkww5RAHZM6AZGy+8C/cFoq2cCKnBn2X2umtYsIpnd8FBQUsXLiQhx56iNatW9OtWzeeeOIJXn99AwMHwurV/4Ux+Tz99D/51a9+FVi/uuTNwY7OJ58KPCMYCosMD+6cUvUH19BO0qjR780WGjBUTNi1uFJRURHr169n69bAJTy7d+9m5MiRLFmyhJSUFObOncuBAwc4eXJgcBGiK7hwoV75gFYyA+uUKd8PZgVmTBey1gnXzZ1y6QdffEC76L2lj/XAF5qadO2Dg+goKRUzsVo46dChQ3z00Uf4/X5WrFhBfn4+SUk/4+WX/0bv3sKGDRtISkqidu3vM7LVTplVNDKmbGd2ZSNnPDCqxnb6HUaU1VFStndUR/Kmnd4109mzZ2Xr1q0i6elSXFwsiYmJAkjLli3lrrvGSJ06b0utWnlVdlBXq1O+oo7Osp3ZlXVse6CT1BaRvHBRlYPFTm/bD/KRvGnAqBmKi4vlq6++kpdfflnuvvtuadCggbRo0SIwkklE/H6/5IAUFxfL9OmBK7khcD99ejU+qLoHoooOZjX5gBbNenvg6monsRowNCVVQ9i9jna4CgoKaNCgAXFxcTzzzDM8+2xgleB27drh8/nw+XzcNXQocSX/noMpi7BGZ1lJe0yZcvl8uof+b4UkmmkjTUlFVMTW9Fbu54R1tKuruLiYTZs2MXPmTPr370/Tpk0pORm45557mD17Njt37mTHAw/w59mzuWfoUOLgkpFIyf/PsP3HU6qe/iNUl+vMrs77VfXpUrD2sNIMcctNU1IVCystE0Ml03xv2bJFWrRoUXrhXPfu3WXSpEmya9euyndgtQ+hMuHkyct+ptV0jNdSK9rP4EpoSkqVcOJFcxBYknT9+vWliwndeeedTJs2jTNnzjB27FgGDx5Mamoq11xzTdU7K0kNXZSSCkt19xHKlB9eTq14uW4eo3NJqVJOnM/ppz/9KUuWLOHbb78lLi6O3r17c/311wNQv3593nzzzertcOrU8mkKO1IWVoPFxf0eJSm09HRNUSlH0xaGiqqSJUn9fj/79u1j8eLFAPz85z+nsLCQtLQ0Bg4cyJVXXhneB0XjbDYWkwR6+SxcJ1l0Db0OQ9nqgw8+kNTU1HJLkqampsrZs2cj9yFeyJd7rQ9DxV4E/r2jCyipWDlx4gTvvvsuY8eOZd++fQD8+9//Zv/+/eWWJM3IyCA+Pj5yH+yk6TZC/Uwd7aPCFcNpUjQlpUJy7NgxXnnlFfx+P+vXr6eoqIgmTZqwcOFCfD5fYERFLNeqtju1U/L5moZxFy/8XhH4t6/XYaiI+ve//80bb7xBRkYGEEhlpqenc+7cOSZPnsyaNWvIy8vD5/MBxDZYgHPO1HVSPHdx6+9l1+zHVvJWbrlpH0YYKsiDfvLJJzJ58mTp0aNH6TURo0aNKt2en58fwwI6kBf6UMqKZLnd8h14oQ8pAnVA55JS1QKya9cuee+990qfSk5Oltq1a0u/fv1kxowZ8vnnn0tRUZGNhXQorwSOSB48nXwg9srvVSKGAUP7MGqwU6dOkZmZid/vJ+PPf2YHULduXfLz82nYsCFfffUVLVu2pHHjxnYXNXzRzlWX5JHt7ksJRyTL7pbvwS3lrEwE/m1rH4a6hIiwZcsWTp48CcDLQ4cyZMgQ5v75z7QHXga2nD9Pgz/8AYAOHTp4I1hA9HPVTulDqa5I5sJ1VUF7xPD71RaGRW6d7TU/P5/ly5eTkZFBRkYGX3/9Ne+++y7Dhw9n//79fPXVV/Tp04d69esH3uChfw/lxOpM0s2jbmpiC8PNv1cEaQsjgtw022thYSHffPMNEFiONCEhgVGjRrFkyRL69OnD3Llz6dOnDwCJiYkMWrv2+2AB3jortOOM1wvfW02iv1f1WOnoiNYNuB/YChQDSRdtexLYCWwHfFb2F61Ob6fP9nrgwAF57bXX5P7775crr7xSxowZIyKB2V9nzpwp69atkwsXLlS+k5KOQK+yo25u60StiaOklIi4pNPbGNMpGCxeAR4XkY3B5zsDbwG3A9cCK4AOIlJU2f6ilZJy2myvxcXF1KoVaBzec889fPDBBwBce+21+Hw+RowYwd133139HbsljRAKO+rm5e9TeYorZqsVkVyo8CKvYcDbInIO2GOM2UkgeNiSDLJ7tlcRYceOHfj9fvx+P1u3bmX37t3ExcXh8/no168fPp+Prl27hnfBnFs7bq3wct2UihGnTm/eClhf5u+DwecuYYwZC4yFQE4+WpKT7WlVLFq0iCeeeIK9e/cC0L59e4YNG8apU6do0qQJEydOjNyHeTmfG6u6RWPqcrs6Zm3uEL5w4QIHDx7k7NmztpXBa+rVq0fr1q2pU6dOSO+PekrKGLMCqGgFnN+KyPvB12RSPiU1G8gSkfnBv+cCH4rIu5V9lpuvwyguLmbz5s2liwnNnDmT5ORkVq1axUsvvURaWho+n48bbrjB7qIqqyKVkrIrtWVzSm3Pnj00btyYZs2axX6qGQ8SEY4fP05BQUHp2jMlHJOSEpFBIbztIHBdmb9bA19HpkTOcuzYMR577DEyMjI4evQoAD169KCgoACAAQMGMGDAgNgWSocaKgc4e/Ysbdu21WARIcYYmjVrxrFjx0Leh1OH1S4FRhlj4o0x1wPtgc9sLlPYLly4wJo1a3j66aeZM2cOAFdccQVr1qxh0KBBzJs3j8OHD7Np0yZSU1PtK6hbJ2RzmnD6Tey6CM7Oi+8q+AwNFpEV7vdp9yip+4D/ARKAE8BmEfEFt/0WeAgoBH4tIsuq2p9TU1Lz589n8eLFrFy5ku+++464uDhGjx7N3LlzgUBT0VH/MXR0j7PUlJTURZ+Xm5tLp06dYvf5NURF36srLtwTkSUi0lpE4kWkRUmwCG77vYjcKCI3WQkWTnH69Gn8fj/PPvts6XPvvfce2dnZjBw5knfffZe8vLzSYAEOOYvSaR0U6O+tKmflYg233OyarXbv3r3ywgsvSGpqqsTHxwsg9erVk6NHj4qISEFBgRQXF9tStpA45QI+vfgrIJbfQ9nfPj09+p9dycyx27Zti+5n11AVfa/oEq3R880337Bo0SIOHToEwMcff8xjjz3GgQMHmDBhAn6/n/z8fBISEgBo1KhR7FsRXjhT1L6UALt+y4uHCEfrM0rCBNi7zG4EnThxorSf0ks0YFhQVFTEZ599xrRp00hJSaF58+bcf//9vPfeewDce++97Nu3j23btvHiiy/i8/moX3Z+JjuE8x9dL3KrWTQdGXGVBQwRobi4OMYligwNGJdx+PBhtm3bBgR+/N69ezNlyhQuXLjAU089xZo1axg3bhwATZo0iepFgzFn54FCD16xd/FZfskJQ8lJR6x+gwicqGRlwYwZkZsgdMuWLaSkpJT+/fnnn1sa5j558mR27dpFjx49Si+87dSpExMmTKBnz56sWbOGrl27lr5+1qxZTAl+v/Pnz+f222+nR48ejBs3jqKiSmdEii0reSu33MLpwzh37pysWrVKJk2aJN27dxdABg8eXLr9ww8/lGPHjoW8/5jQlcTcyUm/z8XfuY2/QXX7MNatE6lfPzBBaP36gb/DVVRUJC1atJDCwkIREenfv79kZ2eLiEifPn2ke/ful9yWL18ue/bskS5dupTuZ8+ePWKMkaysrNK/y25//vnnJT3YbzNkyBA5f/68iIj8/Oc/l3nz5oVfkTLC6cNw6tQgMTdkyBCWL19O7dq1SUlJYcaMGdx5552l28s+dqyyF9zp0NjQ2HHR4tSpzmlBuTgdmZkZmCC0qChwn5kZ/nQ+tWrVokuXLmzdupUdO3aQmJhIz549AVizZs1l31cylU9Zbdq0oXfv3pV+3sqVK8nOzua2224D4MyZM1x99dWhVyDCNGAEPfroo/ziF79gwIAB3lllzu3sOHg56eBth4vr7qIA0r9/YDbpklml+/ePzH579+7Np59+ypw5c/D7/aXP9+3bt3RGhrJmzZpFu3btLnm+YcOGpY9r165drh+jZL4sEWH06NHMmDEjMoWPNCvNELfc7BpW60hOSnO4SaxSMF5LH0ZBKMNq160LrFcTiXRUiaVLl0rTpk3lmWeesfyevLw8SUxMLP374hTU+fPnpVmzZpKXlydnz56VXr16SXp6umzdulXatWsnR44cERGR48ePy969eyNXGdFhtaoiNfksubrsWplPvDec1G7JyfDkk5GdWbpjx47Ex8czadIky+9p1qwZKSkpdO3alSeeeOKS7XXq1OF3v/sdvXr1YsiQIXTs2BGAzp078+yzz5Kamkq3bt0YPHgwhw8fjlhdwqVreitVli605BhOmRpk4sSJ3HbbbYwePdruokSEa6cGUUrhqn6CmmTXrl107NiRM2fOeCZYhEs7vZUqy46Dt6ahHOnGG2/kyy+/tLsYjqItDKXK0oO3UpelAUMppZQlGjBUZOiZuVKepwFDRYbOLKuU52nAUEopZYkGDBU6nVlWqRpFh9Wq0Olkh0rVKNrCUEopZYkGDBUZerWyUuW8/PLLdOrUiQceeKDS1zVq1ChGJQqfBgwVGdpvoZzCIf8W58yZw4cffsiCBQvsLkrEaMBQNYtDDiYqiiI4xDvUJVrHjx/P7t27GTp0KH/84x8BuPfee7n11lvp0qULr776arnXnzp1irvvvpvu3bvTtWtX3nnnndJtjlqy1coc6G656XoYqko1ZdlXDwhlPQwRiehvHOoSrSIibdq0Kbes8/Hjx0VE5PTp09KlSxfJy8sTEZGGDRvKokWL5Gc/+1npa0+cOCEiEpUlW3WJVqVUzTZlSvmWRclQ7/T0sFqVoS7RWpGXX36ZJUuWAHDgwAF27NhBs2bNALj55pt5/PHHmTRpEkOGDKFv376A85ZstTUlZYx53hjzpTHm/4wxS4wxV5bZ9qQxZqcxZrsxxmdjMZXb6fUi7hLK7xLFBalKlmidMmUK06dPL32+b9++9OjR45LbihUrLtlHZmYmK1asICsri5ycHG655ZbSZVkBOnToQHZ2NjfffDNPPvkk06ZNC1YjsGTr5s2b2bx5M9u3b2eKnf9urTRDonUDUoHawcfPAc8FH3cGcoB44HpgFxBX1f40JaWqpCkp5wv+Rk5ISYmEtkSrSPmU1HvvvSdDhgwREZHc3FyJj4+X1atXi0ggJXXo0CE5c+aMiIgsWbJEhg0bJiISlSVbXZuSEpGPyvy5Hvhh8PEw4G0ROQfsMcbsBG4HsmJcRKWU20R4iHcoS7ReLC0tjb/+9a9069aNm266id69e5fbvmXLFp544glq1apFnTp1+Mtf/gKUX7K1uLiYOnXqMHv2bNq0aRNWnULlpD6Mh4CSoQGtCASQEgeDz13CGDMWGAuQmJgYzfIpL9DrRZypoj6IZcvgiivg2murv68Ieumll5gxYwYNGzas1vv27t1b+jg+Pp5ly5ZV+LqTJ08C4PNVnHkfOXIkI0eOrNZnR0vU+zCMMSuMMV9UcBtW5jW/BQqBkgHLpoJdVTjvhIi8KiJJIpKUkJAQ+Qoob9F+C2eqqA+iTZvqB4sI0iVaLxX1FoaIDKpsuzFmNDAEGBjMpUGgRXFdmZe1Br6OTgmV8qCy83ypkOgSrZeye5RUGjAJGCoip8tsWgqMMsbEG2OuB9oDn9lRRqXKcctB2M3rk2ja0LHsvtL7z0BjYLkxZrMx5q8AIrIVWAhsA/zAL0TExssblQpy84HYLdwSlGsgWwOGiLQTketEpEfwNr7Mtt+LyI0icpOIVNxbpJT6nl5voqLM7haGUs7nlgNxFC9eUwqcNaxWKWfShaKUArSFoZQ3acexigINGEpVh1sOxJqGUlGgAUOp6tADsarBNGAopZTDPfTQQ1x99dV07dr1km3jxo3j448/5o477qBTp0506dKFl156KSrl0IChlFION2bMGPx+f4XbNmzYQLt27XjhhRfIzc1l/fr1zJ49m23btkW8HBowlFKqEqNGjWLkyJH06tWLNm3a8M9//jPmZfjBD35A06ZNL3k+NzeXDh060KpVq9KFnRo3bkynTp04dOhQxMuhw2qVUq7Qv3//S5770Y9+xIQJEzh9+jR33XXXJdvHjBnDmDFjyMvL44c//GG5bZmZmZY+Nycnh3vvvZd33nmHtWvX8uijj3L33XfzzTffcNVVV4VSFSCwAFNBQcElz8+aNYtBgyqdgq/UsmXLSEtLK/fc3r172bRpE7169Qq5bJejAUMppS7jzJkz5OXlkR4cHde5c2e++eYbAB555BH+8Y9/lHu9iGBMRZNtX6q6S7xWJCMjg9dff73075MnTzJixAj+9Kc/0aRJk7D3fzENGEopV6isRdCgQYNKtzdv3txyi6KsL774gvbt21OvXj0APv/8c7p3747f7+fLL79k1qxZPPjggwwfPpyhQ4eSkpLChg0bePzxx5kwYQIzZszgueee4/Tp05w/f545c+aU7jvcFsbp06c5ceIE1wangL9w4QIjRozggQceYPjw4dWuqxUaMJRS6jJycnLYv38/Z8+epaioiPT0dP7whz8QHx/Pgw8+yMSJE1m2bBmjRo3iV7/6FW+88Qbdu3cHAgf0+fPnc+bMGa688kp2795dbt/htjBWr17NHXfcAQRaNg8//DCdOnXi0UcfDWu/ldGAoZSX6boYYcnJyeGBBx6gf//+fPfddzz11FOkpKTw97//vTQwbN68mXvvvRcILLXq8/n47rvvMMawadMmZs+eTXx8fFjl+PGPf0xmZiZ5eXm0bt2aqVOnsmnTptJ+mU8//ZQ333yTm2++mR49egAwffr0Cvt1wqEBQykvmzpVA0YYcnJy+Nvf/sZzzz1X7vnmzZvz2muv0bx5c3bs2MFNN90EBNb/njVrFrVr16Zjx4507tyZMWPGcN111zFgwIBLOqiteuutty55rmfPnvzxj38EoE+fPkgM5jgzsfiQWElKSpKNGzfaXQylnMPFkyXm5ubSqVMnW8vQqlUrDhw4QK1a3rkCoaLv1RiTLSJJVb3XO9+CUirALdOxu8ChQ4c8FSzCpSkppbxGp2NXUaKhUymllCUaMJTyMrdMx34ZXupjdYJwv08NGEp5mYv7LerVq8fx48c1aESIiHD8+PHSixBDoX0YSilHat26NQcPHuTYsWN2F8Uz6tWrR+vWrUN+vwYMpZQj1alTh+uvv97uYqgyNCWllFLKEg0YSimlLNGAoZRSyhJPTQ1ijDkG7LO7HFHQHMizuxBR4uW6gdbPzbxcNyhfvzYiklDVGzwVMLzKGLPRyjwvbuTluoHWz828XDcIrX6aklJKKWWJBgyllFKWaMBwh1ftLkAUebluoPVzMy/XDUKon/ZhKKWUskRbGEoppSzRgKGUUsoSDRgOZYx53hjzpTHm/4wxS4wxV5bZ9qQxZqcxZrsxxmdjMUNmjLnfGLPVGFNsjEm6aJsX6pcWLP9OY8xku8sTLmPM340xR40xX5R5rqkxZrkxZkfw/io7yxgOY8x1xpjVxpjc4L/L/wo+7/o6GmPqGWM+M8bkBOs2Nfh8teumAcO5lgNdRaQb8BXwJIAxpjMwCugCpAFzjDFxtpUydF8Aw4FPyj7phfoFyzsbuBPoDPw4WC83+weB36OsycBKEWkPrAz+7VaFwGMi0gnoDfwi+Jt5oY7ngAEi0h3oAaQZY3oTQt00YDiUiHwkIoXBP9cDJXMSDwPeFpFzIrIH2AncbkcZwyEiuSKyvYJNXqjf7cBOEdktIueBtwnUy7VE5BMg/6KnhwHzgo/nAffGskyRJCKHReTz4OMCIBdohQfqKAEng3/WCd6EEOqmAcMdHgKWBR+3Ag6U2XYw+JxXeKF+XqiDFS1E5DAEDrjA1TaXJyKMMW2BW4ANeKSOxpg4Y8xm4CiwXERCqpuuh2EjY8wK4JoKNv1WRN4Pvua3BJrLC0reVsHrHTk22kr9KnpbBc85sn6V8EIdaiRjTCPgXeDXIvKdMRX9lO4jIkVAj2Bf6BJjTNdQ9qMBw0YiMqiy7caY0cAQYKB8f8HMQeC6Mi9rDXwdnRKGp6r6XYZr6lcJL9TBiiPGmJYictgY05LA2atrGWPqEAgWC0RkcfBpT9VRRE4YYzIJ9EdVu26aknIoY0waMAkYKiKny2xaCowyxsQbY64H2gOf2VHGKPFC/f4FtDfGXG+MqUugE3+pzWWKhqXA6ODj0cDlWo2OZwJNiblAroi8WGaT6+tojEkoGWVpjKkPDAK+JIS66ZXeDmWM2QnEA8eDT60XkfHBbb8l0K9RSKDpvKzivTiXMeY+4H+ABOAEsFlEfMFtXqjfXcCfgDjg7yLye3tLFB5jzFtAfwJTYh8B0oH3gIVAIrAfuF9ELu4YdwVjTB9gDbAFKA4+/RSBfgxX19EY041Ap3YcgUbCQhGZZoxpRjXrpgFDKaWUJZqSUkopZYkGDKWUUpZowFBKKWWJBgyllFKWaMBQSilliQYMpZRSlmjAUCpKgtNlDw4+ftYY87LdZVIqHDo1iFLRkw5MM8ZcTWAyu6E2l0epsOiFe0pFkTHmY6AR0F9ECowxNwC/Ba4QkR/aWzqlqkdTUkpFiTHmZqAlcC64xgLBNTIetrdkSoVGA4ZSURCc/XMBgUVqTrl1qVmlytKAoVSEGWMaAIsJLPmZC/w3MMXWQikVAdqHoVQMBWcI/T0wGHhNRGbYXCSlLNOAoZRSyhJNSSmllLJEA4ZSSilLNGAopZSyRAOGUkopSzRgKKWUskQDhlJKKUs0YCillLJEA4ZSSilLNGAopZSy5P8DZLpuA+ow/pYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x_true = x_data[y_data]\n", "x_false = x_data[~y_data]\n", "plt.plot(x_true[:,0], x_true[:,1], '.b')\n", "plt.plot(x_false[:,0], x_false[:,1], '+r')\n", "\n", "r = np.linspace((-1,-1), (1, 1)) * 3 * x_scatter / np.linalg.norm(slope_data)\n", "x_r = r * (np.array([[0,1],[-1,0]]) @ slope_data) - intercept_data / np.linalg.norm(slope_data)**2 * slope_data\n", "plt.plot(x_r[:,0], x_r[:,1], '--k');\n", "plt.xlabel(\"$x_1$\");\n", "plt.ylabel(\"$x_2$\");\n", "plt.legend([\"$y$=true\", \"$y$=false\", \"$p_{true}=1/2$\"]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Logistic regression model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let us bulid an train a logistic regression model." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with Graph(\"graph\") as graph:\n", " \n", " x = Variable(\"x\", shape=(2,), mean=0, variance=x_scatter**2)\n", " y = Variable(\"y\", shape=(), distribution=BernoulliDistribution)\n", "\n", " connect_via_regression(\n", " name_prefix=\"parameters\",\n", " inputs=x,\n", " outputs=y,\n", " order=1,\n", " ) \n", " \n", "trained_graph = Trainer(graph=graph, data = {graph.x: data[[\"x_1\", \"x_2\"]], graph.y: data[\"y\"]})()\n", "\n", "trained_graph.parameters_y.location.slope.mean" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use the trained graph to predict the probability of y=true for a given set of x-values:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x_1x_2p_y_pred
0-10.00.00.913
1-8.00.00.854
2-6.00.00.812
3-4.00.00.718
4-2.00.00.619
50.00.00.505
62.00.00.396
74.00.00.293
86.00.00.198
98.00.00.150
1010.00.00.094
\n", "
" ], "text/plain": [ " x_1 x_2 p_y_pred\n", "0 -10.0 0.0 0.913\n", "1 -8.0 0.0 0.854\n", "2 -6.0 0.0 0.812\n", "3 -4.0 0.0 0.718\n", "4 -2.0 0.0 0.619\n", "5 0.0 0.0 0.505\n", "6 2.0 0.0 0.396\n", "7 4.0 0.0 0.293\n", "8 6.0 0.0 0.198\n", "9 8.0 0.0 0.150\n", "10 10.0 0.0 0.094" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "prediction_x_data = np.linspace((-10, 0), (10, 0), 11)\n", "\n", "predictor = Predictor(graph=trained_graph, data={trained_graph.x: prediction_x_data}, n_samples=1000)\n", "prediction_y_data = predictor(trained_graph.y)\n", "\n", "prediction_data = pd.DataFrame()\n", "prediction_data[\"x_1\"] = prediction_x_data[:,0]\n", "prediction_data[\"x_2\"] = prediction_x_data[:,1]\n", "prediction_data[\"p_y_pred\"] = prediction_y_data\n", "display(prediction_data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also compare the predicted vs. true $p_{true}=1/2$-line:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEICAYAAACqMQjAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6q0lEQVR4nO3dd3hUZfrw8e+TRkJCEZIAEprUQAJI79LEQELVRVZQUFZ00d21A/pTQFHQxbIW9EUFdWERXRaF0BQEBQwoPWBAWuiQAKGFkDLzvH/MZEwgE5JMOTOT+3NduaacmTn3M0nOfc5TldYaIYQQoih+RgcghBDCc0mSEEIIYZckCSGEEHZJkhBCCGGXJAkhhBB2SZIQQghhV4DRAThTeHi4rl+/vtFhCCGEV9m6detZrXVEUdt8KknUr1+fLVu2GB2GEEJ4FaXUEXvbpLpJCCGEXZIkhBBC2CVJQgghhF0+1SZRlNzcXI4fP861a9eMDkUYIDg4mKioKAIDA40ORQiv5PNJ4vjx41SqVIn69eujlDI6HOFGWmvOnTvH8ePHadCggdHhCOGVfL666dq1a1SvXt2jE8SVK3DqlOVWOI9SiurVq8tVpBAO8PkrCcDjE8Tvv4PZDH5+0KQJhIUZHZXv8OTfvRDewOevJDzd5cuWBAGW28uXjY1HCCEKkiRhsEqVLFcQYLmtVMnYeIQQoiDDk4RSKlgp9YtSaqdSao9Saqr1+WpKqe+VUvutt7cYHasrhIVZqphq176xqunChQvMmjXLuOCEEOWe4UkCyAZ6a61bAa2BOKVUJ2AisEZr3RhYY33sk8LCoFatG9siiksSWmvM+fVUQgjhIoYnCW2R368n0PqjgcHA59bnPweGuCumpCSYPt1y6wzJycl07drV9njbtm307t37pu+bOHEiBw8epHXr1jz77LOkpqYSHR3N+PHjadOmDevXrycmJsb2+pkzZzJlyhQA5s2bR4cOHWjdujWPPPIIJpPJOYURQpQrHtG7SSnlD2wFGgEfaK03K6VqaK1PAWitTymlIt0RS1IS9OkDOTkQFARr1kDnzo59ZosWLTh48CAmkwl/f3+efvpp3nzzTQC6d+/O5SJaq2fOnMmMGTPYvXs3O3bsACA1NZV9+/Yxd+5cZs2aRWpqapH7S0lJYeHChWzcuJHAwEDGjx/P/PnzeeCBBxwriBCi3PGIJKG1NgGtlVJVgcVKqZibvMVGKTUOGAdQt25dh2NZt86SIEwmy+26dY4nCT8/P1q0aMGePXvYv38/devWpU2bNgCsX7/e7vuKSgL16tWjU6dOxe5vzZo1bN26lfbt2wOQlZVFZKRbcqwQwsd4RJLIp7W+oJRaB8QBZ5RStaxXEbWANDvvmQ3MBmjXrp12NIaePS1XEPlXEj17OvqJFp06dWLjxo3MmjWLlStX2p4v7kqiUaNGNzwfGhpqux8QEFCoXSJ/0JjWmtGjRzN9+nTnBC+EKLcMTxJKqQgg15ogQoC+wOvAEmA0MMN6+6074unc2VLFtG6dJUE4ehWRr1OnTowZM4bHHnuM2rVr254v7kri3LlzRSaQfDVq1CAtLY1z584RFhZGYmIicXFx9OnTh8GDB/Pkk08SGRnJ+fPnuXz5MvXq1XNOYYQQ5YbhSQKoBXxubZfwA77SWicqpZKAr5RSY4GjwJ/cFVDnzs5LDvmaNWtGhQoVmDBhQonfU716dbp27UpMTAz9+/fnscceK7Q9MDCQl156iY4dO9KgQQOaNWsGQPPmzZk2bRr9+vXDbDYTGBjIBx98IElCCFFqSmuHa2g8Rrt27fT1K9OlpKQQHR1tUER/ePzxx2nfvj2jR482OpRyx1P+BoTwVEqprVrrdkVtM7wLrK87ePAgzZo1IysrSxKEEMLreEJ1k09r2LAhe/fuNToMIYQoE7mSEEIIYZckCSGEEHZJkhBCCGGXJAkhhBB2SZIQQghhlyQJIYQQdkmSEEIIYZckCSGEEHZJkvBw7777LtHR0YwcObLY14Vdv6ydEEI4gSQJe6wrvBlt1qxZLF++nPnz5xsdihCiHJIkYc/UqU77qLIuX/roo49y6NAhBg0axNtvvw3AkCFDaNu2LS1atGD27NmFXp+ZmUl8fDytWrUiJiaGhQsX2rbJcqZCiLKQuZvcoKzLl3700UesXLmStWvXEh4eDsCcOXOoVq0aWVlZtG/fnrvvvpvq1asDsHLlSm699VaWLVsGwMWLFwFZzlQIUXaSJAqaMqXwFYRSltvJkx2qfirr8qVFeffdd1m8eDEAx44dY//+/bYkERsbyzPPPMOECRNISEige/fugCxnKoQoO0kSBU2Z8kcyUAqcuNZGWZYv7du3b6Hn1q1bx+rVq0lKSqJixYr07NnTtmQpQJMmTdi6dSvLly9n0qRJ9OvXj5deekmWMxVClJkkCTcpy/Kl17t48SK33HILFStWZO/evWzatKnQ9pMnT1KtWjVGjRpFWFgYn332GYAsZyqEKDNJEvZMnuzUjyvL8qXXi4uL46OPPqJly5Y0bdqUTp06FdqenJzMs88+i5+fH4GBgXz44YeALGcqhCg7Wb7UTWT5UuN4yt+AEJ5Kli81kCxfKoTwZlLd5GKyfKkQwpvJlYQQQgi7JEkIIYSwS5KEEEIIuyRJCCGEA5KSYPp0y60vkoZrIYQoo6Qk6NMHcnIgKAjWrIHOnY2OyrkMv5JQStVRSq1VSqUopfYopf5hfb6aUup7pdR+6+0tRscqhBAFrVtnSRAmk+V23TqjI3I+w5MEkAc8rbWOBjoBjymlmgMTgTVa68bAGutjIYRBfL1apSx69rRcQfj7W2579jQ6IuczvLpJa30KOGW9f1kplQLUBgYDPa0v+xxYB5R9TgshRJmVh2qVsujc2fJdrFtnSRC++J14wpWEjVKqPnA7sBmoYU0g+YmkyLmtlVLjlFJblFJb0tPT3RarN3rooYeIjIwkJibmhm2PPPIIP/74I7169SI6OpoWLVrwr3/9y4AohScqD9UqZdW5M0ya5JsJAjwoSSilwoBFwBNa60slfZ/WerbWup3Wul1ERITrAvQBY8aMKTRNeUGbN2+mUaNGvPnmm6SkpLBp0yY++OADfvvtNzdH6Zu8vaqmPFSriKIZXt0EoJQKxJIg5mut/2d9+oxSqpbW+pRSqhaQZlyEjhsxYgRaa1JTUzl9+jSzZs0iPj7erTH06NGD1NTUG55PSUmhSZMm1K5d2zaNeaVKlYiOjubEiRM0b97crXH6Gl+oqikP1SqiaIYnCaWUAj4FUrTWbxXYtAQYDcyw3n7rjP31LOIUaPjw4YwfP56rV68yYMCAG7aPGTOGMWPGcPbsWe65555C29aV8Lp7586dDBkyhIULF7Jhwwaeeuop4uPjycjI4JZbyt5xqzSLFtmzYsUK4uLiCj2XmprK9u3b6dixY5ljExZFVdV440G2c2fvjFs4xvAkAXQF7geSlVI7rM89jyU5fKWUGgscBf5kTHiOy8rK4uzZs0y2rlHRvHlzMjIyAHjyySdtiwPl01qj8pdOvYnSLn9alFWrVjF37lzb4ytXrnD33XfzzjvvULlyZYc/v7zLr6rJv5KQqhrhTQxPElrrDYC9I2IfZ++vuDP/ihUrFrs9PDy8xFcOBe3evZvGjRsTHBwMwLZt22jVqhUrV65k7969zJw5k1GjRjFs2DAGDRpE165d2bx5M8888wzjx49n+vTpvP7661y9epWcnBxmzZpl+2xHrySuXr3KhQsXuPXWWwHIzc3l7rvvZuTIkQwbNqzUZRU3cnVVTVKSVAMJ1zE8SZQHO3fu5OjRo1y7dg2TycTkyZN54403qFChAqNGjeLxxx9nxYoVjBgxgr///e988cUXtGrVCrAcxOfNm0dWVhZVq1bl0KFDhT7b0SuJtWvX0qtXL8ByBTN27Fiio6N56qmnHPpcUZirqmp8ob1DeDZJEm6wc+dORo4cSc+ePbl06RLPP/88Xbt2Zc6cObZksGPHDoYMGQJYliG96667uHTpEkoptm/fzgcffECFChUciuPPf/4z69at4+zZs0RFRTF16lS2b99ua2fZuHEj//73v4mNjaV169YAvPbaa0W20wjP4CvtHcJzSZJwg507d/Lxxx/z+uuvF3o+PDycTz75hPDwcPbv30/Tpk0By3rYM2fOJCAggGbNmtG8eXPGjBlDnTp16N279w2NzCW1YMGCG55r06YNb7/9NgDdunXDl5azLQ+8pb1DqsS8l6xx7Qa1a9fm2LFj+Pl5zLCUcsUT/gZcydMPwFIl5vmKW+NariTc4MSJE0aHIHyYp3dNlSox7yantkIIl5LR2t5NriSEMJinVxc5SkZrezdJEsKprlyBy5ehUiUICzM6Gs9XXurrPb1KTNgn1U3Caa5cgd9/hxMnLLdXrhgdkeeT2VWFp5MkIZzm8mUwmy33zWbLY1G88lZf7+2z4ZZHUt0knKZSJfDzsyQIPz/LY1G88lRfX16q1nyNJAnhNGFh0KSJtEmUVnmpr5eusN5JkoRwqrAwSQ6iaN4yOlwUJm0SbtKlS5ebvubdd98lOjqakSNHujyeKVOmMHPmTJfvpyyKW2YVZKlVb5VftfbKK1LV5E0kSbjJzz//fNPXzJo1i+XLlzN//vybvlZrjTm/lbiIx96suGVWQZZa9Wa+vh60L5Ik4SZhYWGkpqYSHR3Nww8/TIsWLejXrx9ZWVkAPProoxw6dIhBgwbx9ttvM2/ePDp06EDr1q155JFHMJlMtvePHz+eNm3asH79+kKPjx07BlDkewFeffVVmjZtSt++fdm3b1+RcY4YMYJ7772Xjh07Uq9ePZYtW+aeL6iAHj16UK1atSK3FVxqtU2bNkDhpVaFEM5V7token7W84bnhrcYzvj247mae5UB84tYvrT1GMa0HsPZq2e556vrli8ds65U+9+/fz8LFizg448/Zvjw4SxatIhRo0bx0UcfsXLlStauXUt6ejrPPfccGzduJDAwkPHjxzN//nx69OjBvn37mDt3LrNmzSI1NbXQY7AcRBcuXHjDe1u0aMGXX37J9u3bycvLo02bNrRt2xYoPADO3jKrgMcvtdqiRUdOnZJGc3fx9ZHiwqLcJQmjNWjQwLZWQ9u2bUlNTb3hNWvWrGHr1q20b98esCx/GhkZSY8ePahXrx6dOnWyvfb6x/bee/78eYYOHUrFihUBGDRoEPDHADizGXJyskhPL3qZVbhxqdXSLLMKrl1qdfr0dzh9urKt+22TJpIoXEm6s5Yf5S5JFHfmXzGwYrHbwyuGl/rK4XoFFw7y9/e3VTcVpLVm9OjRTJ8+vdDzqamphIaGFnru+sf23vvOO+8UeUAvOABu//7d1K9/4zKrgG2p1WnTprF8+fIil1mdOXMmgYGBTJ482e1LrfbuPYz82qb8gXySJFxHurOWH9Im4YH69OnDf//7X9LS0gA4f/48R44ccei9PXr0YPHixWRlZXH58mWWLl0K/DEADuDAgZ2cOmVZZjUzM5PJkyfz5JNPApYFkkaNGkXbtm0ZMWIEEydO5PDhw4WWWa1YsSKzZ8+2LbV65bp5OdavX8+OHTtu+ClpVVNxS60WLIcM5HO98jZSvDwrd1cS3qB58+ZMmzaNfv36YTabCQwM5IMPPqBmzZplfm+nTp249957ad26NfXq1aN79+5A4QFw6ek7GTXqxmVWAXbt2kWrVq3YsGGD3WVWAacstVrUMqtjx45lxYoVxS61+uKLr9GlywBpk3CD8jRSvLyTlemETY8ePfj4449ty6gWtGTJEhYtWsTRo0dZs2YNfn5+fPrpp+zdu5eAgACqVq3KhAkTWLp0Kf/5z38cXmq1KG3atGHz5s0EBgaW6n2O/g1IA63wdcWtTCdJQtgUt8yqN0wBbi9GR/4GpIFWlAfFJQlpkxA2J06csJsgPH0KcFfFKFN5C6OdPn2aOXPmMGzYMBYvXuz2/UuSEDflDVOAXx/juXPO+VxnNNDK9NiiLL788ks6dOhArVq1GDt2LFu2bOHixYtuj8Pwhmul1BwgAUjTWsdYn6sGLATqA6nAcK11hr3PEK7lDVOAV6oESkF+7enZs1C9uuNVY4420Ep1lSiJK1eusHr1apYtW8bLL79MrVq1uHTpEgEBAbz66qskJCQQGxtbqnFJzmJ4kgA+A94Hvijw3ERgjdZ6hlJqovXxBANiE3jHFOBhYRAeDunplsdaO2+shCNTeTt7PIGvNaL7WnlKIyMjg/nz55OYmMjatWvJycmhcuXKjBgxglq1ajFu3DjGjRtndJjGJwmt9U9KqfrXPT0Y6Gm9/zmwDgeSRGlHBosbecMU4NWrW6qZCl7xGN0xw5nTY/vaVYmvledm8vLySEpKIigoiI4dO5KZmcnf/vY3mjZtyuOPP05CQgLdunUrde89VzM8SdhRQ2t9CkBrfUopFWnvhUqpccA4gLp1696wPTg4mHPnzlG9enVJFD7o+h5NBa94QkM1586ds40gN0Jx1VWlPYv2tVHOvlaeomRkZLBy5UoSExNZsWIFGRkZDBkyhMWLFxMVFcXhw4epX7++0WEWy1OTRIlprWcDs8HSBfb67VFRURw/fpz0/HoI4TOys+HMGUvVklJQowbkj+G7cMFyGxwcTFRUlGExQtHVVWU5i/a1RXt8rTxguXI9fvw4UVFRKKXo168fW7ZsITw8nIEDBzJw4EDuvPNO2+s9PUGA5yaJM0qpWtariFpAWlk/KDAwkAYNGjgxNOEppk+HF1+0nIn6+1sWs5k0yeioSqYsZ9G+NsrZV8qTnZ3Njz/+SGJiIomJiZw5c4azZ88SEhLCjBkzCA0NpX379vj7+xsdapl4apJYAowGZlhvvzU2HOGJvPlMtKyx+9p62N5engULFvDwww+TmZlJcHAwffv25bnnnrMtANanTx+DI3Sc4UlCKbUASyN1uFLqODAZS3L4Sik1FjgK/Mm4CIWn8uYzUW+OvTwym81s376dxMREli1bxtSpU+nfvz/NmzfngQceICEhgV69ehESEmJ0qE7n89NyCCG8l9FdZC9fvsxTTz3FsmXLOHXqFEopOnXqxEsvveTUecmMVty0HIZfSQghRFGM6CJ7+PBhli1bhlKKxx57jNDQUNavX0/37t2Jj4+nf//+REREuDYIDyNJQgjhkdzVRfbXX39l0aJFJCYmsmfPHgB69erFY489hp+fHykpKR7ZfV5rzZGLR6gUVInqFau7bD8yd5MQLuCJ8zV5YkzFcdXCRhkZGXz99de2gZazZ8/mzTffpEaNGrz99tv8/vvv/PDDD7bXe0qCuJR9ifc2v8e4pePo8mkXqsyoQoN/NWDhnoUu3a+0SQjhZJ44ktgTYyoJZ7RJaK3Zu3evrYvqxo0bMZlMbNu2jdtvv52TJ08SGhpKlSpVnBl6mVzNvcqetD3sTttNcloyu9N2069hP57p8gyXsy9TeUZlqoVUIzYy1vJTI5beDXrTqFojh/YrbRJCuFFJqkmc3SB7s8/z1tHNZe0im52dTXZ2NpUrV2bVqlX0798fgFatWjFhwgQGDhxIy5YtAWxrprtTnjmPA+cPkHwmmQC/AIZGD0VrTdRbUWRcs8xlGhIQQvOI5gT5BwFQqUIlTj99msjQSLde3UiSEKIAZxy8bzYGwtln9SX5PG8eU1JSp0+fZvny5SQmJvLdd98xYcIEXnzxRbp3786HH35IfHw8derUcWtMWmvOZ523tRm8sOYFlh9YTkp6CtmmbADa3dqOodFDUUrxzzv/yS0htxAbGcttt9yGv1/hAXg1wmq4NX6QJCGEjbMO3jcbA+Hss/qSfJ7R4zJc2ZVVa03Pnj356aefAMtUPPfffz+9evUCIDQ0lEcffdS5O7Uj+UwyG45uIDkt2VZdFBwQzKmnTwFwMfsiNcNq0rdBX2JrxBITGUN0+B+rJo5tM9YtcZaGJAnhFdzRX96ZB+/iqkmcfVZf0s8zanSzM6+crly5wpo1a2zTXyxZsgSlFB07duSuu+5yy7oL1/KukZKeYms32JO+h/8N/x8VAiowZ/sc3tn8DlUqVCG2RiwjWowgtkYsZm3GT/nx/oD3XRaXq0iSEB7PXY2u7qqScfZZvdFXCTfjjOS7bNky3nvvPdatW2dra4iLiyMvL4+AgADeeOMNp8dtMps4fOEwyWeS6Vm/J7eE3MLHWz/m0WWPYtaWaTeC/IOIDo8mLTONOlXq8HSXp3mq81NEVY7ymF5RjpIkIQznKY2u7jzYOvus3pPnQCpt8s3Ly2PTpk0kJibyj3/8g1q1anHs2DEOHz7M+PHjGThwoFPXXdBaY9Zm/P382Xt2L29sfIPktGR+S/+Nq7lXAVh+33L6N+5P21vb8ny354mtYeld1KhaIwL9/4gjqrKxMw67gnSBFYYqyVWCt3bfFH+42YlAZmYmS5cuta27cP78eQICAli8eDEJCQmYTCanzKKaY8ph68mtljaDM3+0G8zoO4O/tPkLu87sot+/+1naCyJibMkgJjKGkEDfm5cpn3SBFR7LGxpdheOuv9LJH7uQm5tLy5YtycjI4M9//jPh4eEkJCQQHx/PXXfdZRu7UNoEkWPKYd/ZfbZk0KpmK0bEjOBy9mW6zOkCQFhQGDGRMQxtNpTG1RoDEBsZy+lnTjun0D5CkoQwlKc3upZnxZ39l6UjQXZ2Nj/99JNtUNuhQ4cYPHgw33zzDVFRUWzfvp3Y2NhSJYSNP5v5Zt0RWra7zP39WqK1ptOnndh2aht55jwAAvwCeLz944yIGUH1itVZft9ymoU3o17VevipwpNO+Eo7gjNJdZMwnNEzfYobFVfFV5rqv0uXLlG5cmUAunXrxsaNGwkODqZPnz62K4bSjl2Yt2seP6b+SNKhZPak7YEKV/A70ZkND/5M587w1KqnCA4IJiYyhtjIWJqGN7UNSBNFk+om4dHkKqEwT0iaxVUDFret4LoLiYmJpKSkkJ6eTkhICBMmTMDPz49evXpRsWJFu/vOzMlkT/oeks8k27qZXsq+xC8P/wLAf3/7LxuObqDytVjUzgfRZ2LgzO2sa2CJ46273nLlV3OjKVMsPz5KkoQQHsRTGumLqwa0t23JkiU8+uijhdZdeP7558nNzSUkJISBAwcW2keeOY/fz/3O7rTd7E7bzUt3vESAXwDPff8cs7bMAixTU7SIbEHLyJaYzCb8/fxZcPcCggOC2bRJ0WeGB4winzpVkgSAUupOYDjwgdZ6h1JqnNZ6tutCE6Lk3Hn27cp9ecocS8V1FujcGebNS+WzzxJJT1/G5ctPAv2oU6cOXbt2ZeDAgYXWXdBac+ziMcIrhhMSGMK3e7/lpXUvsffsXnJMOQD4K3/GtB7DbbfcxoO3P0jf2ywjkhtUbXDD1BT5vYykQ4N7lOZKYjzwIPB/SqlqQGuXRCREKZXl7LusB3pXn+l70hxL11cDXr16lZdffrnQuguNGzfm0qVLANx+++18/fXXnL5ymq9++4rkzcm2q4SL2RdZff9q+tzWh9CgUG6tdCt3NbzLNpNps/BmBAcEA5a5jNrdWmT1+E1jdJspUyxXEPnyG7wnT/a5q4rSJIl0rfUF4Bml1AygvWtCEqJ0Snv27ciB3tVn+o6eHTvzKicjI4NVq1aRlZXFgw8+SHBwMAsWLKBx48bc/9D9NOnShMshl9l8ZjOfzPuER9o+wtDooZy4dILHVzxO1eCqxEbGcl/sfbYGZIC+t/Wl7219HQvOaAXbIZQCH+oAdL3SJIll+Xe01hOVUn9zQTxClFppz74dOdC740y/rGfHzrjK2bdvH0uWLPlj3QWzieiu0VTpVIXI0EgOHDjAhZwL1HyzJuZVlqkpKvhXoHlEc9uspi1rtOTYk8eoXam2dCn1ATdNEkqpd4AntdbfFnxea/2eq4ISvseV9filPft25EDvyfXgZUl+2dnZbNiwgV4//cjVF55j5syZfPLJJ1QbXY2IXhGc9z9PijmFu7+6m/ti76PbsG6EB4Tzcs+XaVK9CbE1LFNTBPgFWM6sYyDQP9Anp6ewa/JkoyNwqZuOk1BKTQNaAfdqra8qpfoBk7XWXd0RYGnIOAnP5Ck9dq6PyRMP9I6YPRsef9ySJCpUsP89nz59mlnfzGLplqXsObuH3FtyqRIJraPvYO4dcwkICOChdQ+hULYpKWJrxNI8ojkVA+13XfX1ahdf5tA4Ca31/yml7gN+VEplA5nARCfHKHyYp/TYKcjXxmYkJcETT1i+Yz8/eOcdaNM+m11n9rHz9E52nd7FxZyL3FfpPss6CyOBxhBYO5BmlZrR5ce9dBsxmgYNGgDw/f3fG1kc4UH8bvYCpVQf4GEsySEC+LvWer2rAxO+w1UL2jtTUhJMn265dcbr3Mmszfxv7SGy6y/BrM1oDfNOPUfoq6G0+qgVD3zzAG9veptfT/7K7W1vZ9q0aSx+dDGHzX8n++VcUp7dy6eJ8GCbhyxXA6XpnTNliuU9+W0P+fd9rIdPeVaS6qYfgJe01huUUrHAv4GntNY/uCPA0pDqJs9Vkuodo6qASlod5gnVZlprlFIkHUvi0+2fWha9SdtDZm4mAOrdg6gLz+DXbCl5NfMIuRxC96bdeWzEYwyKH2T/g51RVSTVTV7L0eqm3gXuJyul+gOLgC7OC7FoSqk44F+AP/CJ1nqGq/cpXONm1TtGHoBLWh3mzmqza3nX2HVmV6HprJPTklk0fBHd6nbj+KXjfLvvW+pWqEuzrGb4n/PnLwP/xemna7N9S00aNPg78fHxdO/e3WnrLojyqdTTcmitT1mroFxKKeUPfADcCRwHflVKLdFa/+bqfQv3M7LdoqS9nVzR/TXXlGubmiI5LZm4RnF0q9uNbae20XWOpW9IxcCKtIhoQULjBCpXqMxPP/3Ex69+Q+Z6M9uythEQEMAdd9zB6D5tCOofBMwqfSDO6KHj4718yqsyzd2ktc5ydiBF6AAc0FofAlBKfQkMBiRJ+CAjRxp37mxp6F20CO6+235ycqT7q9aaoxePYtZmGtzSgIysDHp+3vOGqSmqh1SnW91utKrRisX3LiYmIobsM9msWL6CUQNHUbNGTZ578yO+/+47lEogKCiBxMR+3HlnFce+BGe0IUg7hE/y5An+agPHCjw+DnQ0KBbhYkaOP8jvGZSTA+vXQ2xs8YmipLF9tOUjtp/abqsuupxzmQdbP8icwXOoGlyVRtUaEdcwzrb6WdPwpgQHBJOTk8PPP/7M2sS1PJ34NIcOHQKgXr16/OlPf6JSpQfx83sYs9kfkwm2bIE773TOdyHE9Tw5SRQ1VPOGVjGl1DhgHEDdunVdHZNwIaO6pZa1qisrN4vf0n+zrX62O303ERUjmDdsHgDvbHqHtMw0YmvE8kCrB4iJjKFTVCfAsrjNouGLbJ91+vRpjh46SpMmTUhPT6dfv362dReeffbZQusu9O1bgenTPWN+J+H7PDlJHAcKrkYSBZy8/kXWmWhng6V3k3tCE76kYFWXvz8cPWq5ushPFCaziQPnD5Cclsypy6f4W0fLjDQDFwxkzeE1AAQHBNM8ojktIlrYPnfTXzZRpUKVIqem0FoXWnfh119/ZdCgQXz77bfUrl2btWvX0qFDhyLXXSh41TXqwBTqdJ7i7K9ECBuPXZlOKRUA/A70AU4AvwL3aa332HuPdIEVZZWUBJ9/oZn79UnyLtxKhSDF3z6fzffnP+K39N9s8xJV8K/A5UmXCfQPZOWBlWTmZBITGUOjao1umNL6erm5ubaeRnFxcaxatQqlFB07diQhIYGBAwfSsmXLkgedPxNpWf+HfXyxHJfw0e+suC6wHpskAJRSA4B3sHSBnaO1frW415fHJOHN00sYHfvhjMOsOrjKVlX065FkssiAN4/jf7U2A6fM5mr9RZbprK3TU0RHRN84NUUxB47U1FSWLVtGYmIiv/zyC8ePHyckJISvvvqKa9euERcXR2RkZNkKkH+FUtb/YRnXUHo++p0VlyTQWvvMT9u2bXV58vPPWoeEaO3vb7n9+WejIyo5d8V+Lfea3nFqh/73zn/r5757Tvef119vP7Vda631vJ3zNFPQladX1l0+7aKHfPKIDuz6nvYLPVe6mOCGp77//nsdExOjsbSj6UaNGuknn3xSnz171vFCTZ5s2ef1P5Mnl+5ziohb3ISPfmfAFm3nuGr4gd2ZP+UtSbz2muUgC5bb114zOqKSc3bsJrNJHzh3QC9OWaxT0lO01lonHUvS/lP9NVPQTEEHvhyoW37YUq89vFZrrfWFrAs6NSNVm81m2+f8/LMllmITxHUH4/OgFyxYoEeOHKlXr16ttdZ6y5Ytunfv3vrNN9/U+/btc6xw1+/bkQThrARTnpSD76y4JOHR1U2lVd6qmzxhmoiyKmvsWmtyzbkE+QdxKfsST658kuS0ZH5L/802NcXLPV/mxTte5HzWed5Kesu2+lnjao0J9HfC6GOlyH7hBd579VUSgQ2ACQgHZg4ezOhvvnF8HyWMA5DqJnfy0e/Ma9skSqu8JQkwvl7fESWJffPxzew6s8s2Ijk5LZm7o+/mo4SPMJlNNHy3IQ2rNbS1G+RPaR0WFObUWLOzs/npp5/IyMhg+L33YjaZiIqKIjIykvidOzn/6M+MHNmBbt2Kb7wuxNFGUKUcWy7TRw94LuWj35kkCeHR1m/M5b/r9lE9ejfZVZIJCQzh/3r8HwAN/tWA1AuphAaGEhMZQ0xkDHGN4rin+T0uj+vMmTMsX76cxMREvktM5EpODo2xdLkDuARcfGgydeZMJcBfl/5qztEDjqNJxkd76riUj35nkiSER9Bac+TiEY5cOMId9e8AIP7jMSw/9h/wzwUsU1P0ua0Pq0atAixXEhGhEdSvWh8/ddOZ7R2Ob8eOHbRu3RqlFIMGPczSpZ8QGRnFsGEJxMfH07t3byqGhtoO7tOnw6TnFQqNvz+88gpMmlTCHfroWanwPsUlCdf+14lyb9WBVTyy9BG6fNqFKjOq0OBfDYibH4fJbAIg4ExH1KanYNE8/P7fTiYHZNoSBEDHqI7cdsttLksQmZmZfPvtt4wbN46oqCjatGnDzp07SUqC7757Gj+/7Vy6dJQHHviQhISEPwa3WddRmPS8pV1Ao8gzKUYdmFL8DmX9BcfI9+R2ciUhHJKVm8We9D2WsQYF2g12PbqLiNAIpv00jbc3vV1orEFsjVg61u6Iv5+/IY3vJpMJf39/kpKS6NWrF9nZ2VSqVIm4uDji4+MZPHgwH35YlRdftEzVccMVwvVVDkox/TVtv23FXhWFXEmUnnxnLiHVTcJheeY8Dpw/YEkEZ5J5oNUDNKzWkLnb5/LQkoeAP6amiI2MZVrvaURVjiLXlEuAX0CRU1Pkc3Xje15eHps2bbJNgTF8+HDuvPMlvvsuk5SUFxk3LoFu3boRFBRUKKYSJ6+iDlwFE4O9A5sc8EpPvjOXcGjRIVG+aK05cfkEFfwrEBEawe603dy/+H5S0lNsU1P4KT9a12xNw2oNubPhnSwavoiYyBga3tLwhqkpStLl1GUT+02ZwriTJ1m0aBHnz58nICCAbt26AU2sCSCUoKC3+Mc/LIng+phKPCttUesoTJ1686oRWX+hZPKnH8mXf8LhSM8uUWJyJVHOZeVm8dmOzwqtfnbh2gVm9JnBhG4TOH3lNGO+GWPrXhoTGUN0eDQhgSFGh16I1pq9e/eybNkyDh48yIcffghKMfqBBwCI/+IL+mVkULVqVaZPx35Vkj2l7dVi78qpvBzYXNULSK4kXEKqm8q5a3nXSElPsSWB3Wm76RTViZfueIkcUw5hr4VRwa8i1U0xtKkdS5/YGHo36E10RLTRod/Utm3b+Pzzz0lMTLStu9C6dWuSkpIIDgn544BS4OBSpnaQkhycrj/jLciH/s9KxFUHc0kSLiHVTeWEWZs5eP4gu9N2k2vOZXiL4QA0fb8pRy8eBSDIP4jo8GgC/QJtj//X/Sh/GlCD4zmKtCB4dg1ERxhWjGKdPn2a5cuXEx8fT40aNfjll1+YPXv2H+su7NtHnXfegRDrlU7BM3rr/c6TJ7NmzRTnt4PYa4copj3G7vtF0aSKzv3szdfhjT/lZe4ms9msz109Z3s87cdput3sdjpkWohtnqIm7zWxbf9s+2f6y+Qv9Z60PTonL+eGz/PkOaDMZrPeunWrnjp1qm7fvr1twrzPP/9ca631lStXdGZmZtFvLjgZW1kmZnNkzp6C+yvpHD/ePnlcOZjjyFchczd5t71n97L+yPpC7QbX8q5xceJF/JQfz373LDvO7CAmIsa2FGbziOaEBoWW6PM9bQ6ozMxMzp07R926dTlx4gRRUVEopejQoQMDBw4kPj6eVq1aFdtjCrjxbN6Rv/XSvr8sVwW+VJXiS2UpB6S6yQvkmHLYd3ZfoUTw76H/pmpwVebvms+09dNsU1MMaTqE2Bqx5JnzCPIP4p/9/unQvo1cXzpfwXUX1q5dS1xcHN988w21a9dm8eLFdO3alYgIB+rA3F1NUdIEIT13hIeTKwk3M2szRy4cITktmU5RnYgMjeTL3V9y/+L7yTPnARDgF0DT6k35373/o0n1Jpy8fJLsvGzqVa3n8qkp3EVrbbsSGDFiBAsXLgSgUaNGxMfHM3ToUO644w7HduLMs1l3tBf40tm3tK94FendZBCzNuOn/Ei9kMqrP71Kcloye9L3cCXnCgAL71nI8BbD2ZO2h/nJ8y2jkSNjaRrelCD/oJt8uvfJyMhg1apVJCYmsn79evbu3UtISAhffPEFZ8+eJSEhgSZNmji2E3s9jLzhzNyXkoRwPhcmXkkSLpZnzmP7qe2FqoqSzyQzsdtEnuj0BKkXUmk3u52tvSA/GbSs0bLE7QbebOPGjbzwwgts2LABk8lEeHg4AwYM4PXXX6dmzZqu27FRB92y/jPL2bcojgv/niVJOEmeOY/95/bbkkCjao0Y3Xo0WblZhE0Pw6zNhASE0CKyBbGRsYyIGUG/hv3I/45v2tDqA/LXXUhMTGTYsGHccccd/Prrr/zlL38hISGBhIQEOnTogL9/KdZdKCujkkT+fuWg75m89fdiUJIwvNuqM3+c1QXWbDbroxeO2tZC1lrrPp/30UGvBNm6mPpN9dNjvx1r275i/wr9+9nfdZ4pzykxeJOcnBw9Z84cPWzYMB0WFqYBHRwcrN9//31jAzOq62V+V1Zv79Lqq7zp9+KmbsVIF9ib+3rP1/xw+AdbldHF7IvERMaQ/NdkACaunohZm21VRdER0QQHBDszfM913ZmX1prt27dz6tQp4uPjMZvN1K5dG39/f9vVQu/evf+YVrs88Oa2kIKccZbt6Wfq3tr2I9VNjnMkSdy/+H6W7ltqazeIjYylVc1WdKnTxclReiGlyLxyhTVr1pCYmMiyZcs4efIkdevWJTU1FaUUJ06c4NZbby0XVWrF8vZk4YwDkScehL399wKSJJzBkSRxNfcqIQEhcpAr4OjRo9SpUwfl58cj48Yxe/bsQusu9O/fn8jISKPDLBtXnu3m/zN74sHyZnw1SRTk6fHZY1DvJt/odO8EFQMrlvsEYTKZ2LhxI5MmTSK2Rg3q1avHTj/Ln8jjs2ezBjj7j3/w1VdfMXr0aO9NEGB/Ij5n8Lb5hZyxWp6suOd6Bn2XciXhAq5eRMcVtm3bxp133mlbd6F79+4kJCQwcuRIatSs6b1nxva4oyyeXjdflPJwJeGNvxcXkysJN8qfB+nFFy23SUlGR1SY1pZ1F2bOnEnPnj2ZPn06AE2bNmXQoEF89dVXnD17lh9++IGnnnqKGh9+aHmjL5whuvts1xu/o/JAfi+lYujcTUqpPwFTgGigg9Z6S4Ftk4CxgAn4u9Z6lSFBltK6dZaJ8kwmy+26dZ5zNTFp0iS+/vprDh48CEDLli2pXr06AKGhocydO/fGN+X/Q+UfYD35DPFmSrKkqDv27cmcUVXmbdVtoliGVjcppaIBM/D/gGfyk4RSqjmwAOgA3AqsBpporU3FfZ4nVDd5yoyqZ86cYfny5aSkpPDGG28AcM8993Dt2jXi4+OJj4+nbt26pftQb08SBbm7LL703Qmf47GzwGqtU6DIkciDgS+11tnAYaXUASwJw8Mqb25k5Iyqv//+OwsXLiQxMZFffvkFgDp16jB58mRCQ0P5+uuvHWuc96UzRF8qixAu5KltErWBYwUeH7c+dwOl1Dil1Bal1Jb09HS3BHcznTtb1kx2dYLIzMxkyZIl5Jf7hx9+YPLkyfj5+fHKK6+wY8cOjhw5QmioZX4oh3tveUN1SUm5oyzObgNx9/fvS79vUWYur25SSq0GiprF7QWt9bfW16yjcHXTB0CS1nqe9fGnwHKt9aLi9uUJ1U2uduTIEdu6Cz/88APZ2dl88sknjB07losXL5Kdne3dXVN9lTf2GpIqsnLD0OomrXXfMrztOFCnwOMo4KRzIvIuJpOJ8+fPExERwenTp6lfvz5gWXdh/PjxxMfH0717dwCqVKninqC8pRFWCOEwT61uWgKMUEpVUEo1ABoDvxgck9tkZGTw5ZdfMmrUKCIjI3n44YcBqFmzJnPnzmXfvn3s37+ft956iz59+hAU5Oa1J1w5EM1XlbUNxIhuu+4eFCcnHB7N6N5NQ4H3gAjgArBDa32XddsLwENAHvCE1nrFzT7PF6qb/vrXv/Lxxx9jMpmoXr06AwYMYNiwYQwZMsTo0P4g1RDG8NXqJvl7MpzHDqbTWi/WWkdprStorWvkJwjrtle11g211k1LkiC8TU5ODqtXr+aJJ56gZcuWZGVlAdC6dWuee+45fv75Z86cOcMXX3zhGQlCpl0oH+T3Ka5nbw5xb/xx1noSrrRt2zY9bNgwXalSJdu6CwMGDNBHjx41OrSSM3o+fqPWiTCaO8pd8Hc7ebLr9ummdRJEySDrSRhDa82OHTtITEykd+/edO3alV9//ZWhQ4faBrT17dvXmHUXHGl8Nrp6wOj9+7Lrv1t3fNfy+zScxw6m80Umk4nly5eTmJhIYmIiJ0+eRClFYGAgXbt2pV27dhw7dsz4GWenTi17kpCBaL7l+rUW8v825fcs8NzeTV7lyJEj/Pjjj4BlwNq4ceNYsGABnTt3Zu7cuZw6dYqJEyfathueIBxlRL21tIm4zpQpf1T4wB/JIT9xuPq7lmTk0aS6qQxMJhObNm2yXS3s3r2bqKgojh49ilKKvXv3ctttt7m/a+rN+MLqXOC71ROeMP7EiOomYThZmc7J/va3v/H+++8XWnchISGBJk2auHzfTuPN//zujN2dB25P+J1cX15PiEm4nCQJJ9uyZQuHDh2iX79+VK1a1eX7cwlv/uf31QO3J/5OPOHqRricx46T8Fbt2rVj+PDh3psgwLvrgX3poOXpbS2eEocwjCSJ8kr++e1z54H7+kbj/Pvy+xEeQqqbhChOea9uEuWCVDcJ4Q28uQpQ+CxJEkIUx50HbqliEh5IkoQQxZEDtyjnJEkIIYSwS5KEcD05GxfCa0mSEK4nK9kJ4bUkSQghhLBLkoRwDU8fSSyEKBEZTCdcTwaJCeHRZDCdEEKIMpEkIVxPRhIL4bUkSQjXk3YIIbyWJAlRvkkCE6JYkiRE+SZjOIQoliQJIYQQdhmaJJRS/1RK7VVK7VJKLVZKVS2wbZJS6oBSap9S6i4DwxS+RsZweC75HXgcQ8dJKKX6AT9orfOUUq8DaK0nKKWaAwuADsCtwGqgidbaVNznyTgJUWoyhsOzyO/DEB47TkJr/Z3WOs/6cBMQZb0/GPhSa52ttT4MHMCSMIQQQriRJ7VJPASssN6vDRwrsO249bkbKKXGKaW2KKW2pKenuzhE4XNkDIfxpPrPo7m8ukkptRqoWcSmF7TW31pf8wLQDhimtdZKqQ+AJK31POv2T4HlWutFxe1LqpuE8HJS3WSI4qqbAly9c6113+K2K6VGAwlAH/1HxjoO1CnwsijgpGsiFMJHTZkiZ+PCYUb3booDJgCDtNZXC2xaAoxQSlVQSjUAGgO/GBGjEHZ5+gHYG8eASPWfxzG6d9MBoAJwzvrUJq31o9ZtL2Bpp8gDntBaryj6U/4g1U3CrTy9asTT4xMew5N7NzXSWtfRWre2/jxaYNurWuuGWuumJUkQQgikEVg4nSf1bhLC83n6QXjKFMvVQ/4VRP59T4lPeB1ZdEiIsvL06hxPj094DI+tbhJCuJA0AgsnkCQhRFl5+kFYqpiEE0iSEKKs5CAsygFJEkIIIeySJCGEEMIuSRJCCCHskiQhhBDCLkkSQggh7JIkIYQQwi5JEkIIIeySJCGEr5PxHMIBkiSE8HXeuK6E8BiSJIQQQtglSUIIX+TpU5oLryFThQvh62TKcHETMlW4EEKIMpEkIYSv8/QpzYVHkyQhhK+TdgjhAEkSQggh7JIkIYQQwi5JEkIIIeySJCGEEMIuSRJCCCHs8qnBdEqpdOBICV4aDpx1cTju5mtl8rXygO+VydfKA75XppKWp57WOqKoDT6VJEpKKbXF3uhCb+VrZfK18oDvlcnXygO+VyZnlEeqm4QQQtglSUIIIYRd5TVJzDY6ABfwtTL5WnnA98rka+UB3yuTw+Upl20SQgghSqa8XkkIIYQoAUkSQggh7CpXSUIp9YpSapdSaodS6jul1K0Ftk1SSh1QSu1TSt1lZJwlpZT6p1Jqr7VMi5VSVQts87ryACil/qSU2qOUMiul2l23zVvLFGeN+YBSaqLR8ZSFUmqOUipNKbW7wHPVlFLfK6X2W29vMTLG0lBK1VFKrVVKpVj/3v5hfd4ry6SUClZK/aKU2mktz1Tr846XR2tdbn6AygXu/x34yHq/ObATqAA0AA4C/kbHW4Ly9AMCrPdfB1735vJYY48GmgLrgHYFnvfKMgH+1lhvA4KsZWhudFxlKEcPoA2wu8BzbwATrfcn5v/9ecMPUAtoY71fCfjd+jfmlWUCFBBmvR8IbAY6OaM85epKQmt9qcDDUCC/1X4w8KXWOltrfRg4AHRwd3ylpbX+TmudZ324CYiy3vfK8gBorVO01vuK2OStZeoAHNBaH9Ja5wBfYimLV9Fa/wScv+7pwcDn1vufA0PcGZMjtNantNbbrPcvAylAbby0TNriivVhoPVH44TylKskAaCUelUpdQwYCbxkfbo2cKzAy45bn/MmDwErrPd9oTzX89YyeWvcJVFDa30KLAddINLgeMpEKVUfuB3L2bfXlkkp5a+U2gGkAd9rrZ1SHp9LEkqp1Uqp3UX8DAbQWr+gta4DzAcez39bER/lEX2Db1Ye62teAPKwlAk8uDxQsjIV9bYinvOYMhXDW+MuF5RSYcAi4Inrahq8jtbapLVujaVGoYNSKsYZnxvgjA/xJFrrviV86X+AZcBkLGd3dQpsiwJOOjm0MrlZeZRSo4EEoI+2VjziweWBUv2OCvLoMhXDW+MuiTNKqVpa61NKqVpYzmC9hlIqEEuCmK+1/p/1aa8uE4DW+oJSah0QhxPK43NXEsVRSjUu8HAQsNd6fwkwQilVQSnVAGgM/OLu+EpLKRUHTAAGaa2vFtjkleW5CW8t069AY6VUA6VUEDACS1l8wRJgtPX+aOBbA2MpFaWUAj4FUrTWbxXY5JVlUkpF5PduVEqFAH2xHN8cL4/RrfJu7gGwCNgN7AKWArULbHsBSy+UfUB/o2MtYXkOYKnv3mH9+ciby2ONeyiWs+9s4AywygfKNABL75mDwAtGx1PGMiwATgG51t/PWKA6sAbYb72tZnScpShPNyzVfrsK/P8M8NYyAS2B7dby7AZesj7vcHlkWg4hhBB2lavqJiGEEKUjSUIIIYRdkiSEEELYJUlCCCGEXZIkhBBC2CVJQgghhF2SJIRwEetU1Hda709TSr1rdExClJbPTcshhAeZDLyslIrEMoHcIIPjEaLUZDCdEC6klPoRCAN6aq0vK6VuwzJyvIrW+h5joxPi5qS6SQgXUUrFYlncJltb1ixAW9aVGGtsZEKUnCQJIVzAOuPmfCyLvmR603KrQhQkSUIIJ1NKVQT+BzyttU4BXgGmGBqUEGUkbRJCuJFSqjrwKnAn8InWerrBIQlRLEkSQggh7JLqJiGEEHZJkhBCCGGXJAkhhBB2SZIQQghhlyQJIYQQdkmSEEIIYZckCSGEEHZJkhBCCGGXJAkhhBB2/X8s70+vclg9BwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x_true = data[data[\"y\"]][[\"x_1\", \"x_2\"]]\n", "x_false = data[~data[\"y\"]][[\"x_1\", \"x_2\"]]\n", "\n", "plt.plot(x_true[\"x_1\"], x_true[\"x_2\"], '.b')\n", "plt.plot(x_false[\"x_1\"], x_false[\"x_2\"], '+r')\n", "\n", "r = np.linspace((-1,-1), (1, 1)) * 3 * x_scatter / np.linalg.norm(slope_data)\n", "x_r = r * (np.array([[0,1],[-1,0]]) @ slope_data) - intercept_data / np.linalg.norm(slope_data)**2 * slope_data\n", "plt.plot(x_r[:,0], x_r[:,1], '--k');\n", "\n", "inferred_slope_data = predictor(trained_graph.parameters_y.location.slope)\n", "inferred_intercept_data = predictor(trained_graph.parameters_y.location.intercept)\n", "\n", "s = np.linspace((-1,-1), (1, 1)) * 3 * x_scatter / np.linalg.norm(inferred_slope_data)\n", "x_s = s * (np.array([[0,1],[-1,0]]) @ inferred_slope_data) - inferred_intercept_data / np.linalg.norm(inferred_slope_data)**2 * inferred_slope_data\n", "plt.plot(x_s[:,0], x_s[:,1], '--g');\n", "\n", "\n", "\n", "plt.xlabel(\"$x_1$\");\n", "plt.ylabel(\"$x_2$\");\n", "plt.legend([\"$y$=true\", \"$y$=false\", \"$p_{true}=1/2$\", \"inferred $p_{true}=1/2$\"]);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7.7" } }, "nbformat": 4, "nbformat_minor": 4 }