Logistic Regression

LogisticRegression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression\n",
    "- http://hosho.ees.hokudai.ac.jp/~kubo/stat/2013/e/kubostat2013e.pdf\n",
    "\n",
    "ロジスティック回帰は、誤差構造が二項分布、リンク関数がlogit関数であるGLMである。    \n",
    "ロジスティック回帰は、データ（の誤差）が二項分布に従うと仮定した場合の一般化線形モデル。0/1データに対する二値分類問題に用いられる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "まず、二項分布は、成功確率qである事象のN回の試行において、y回成功する確率、を表す。これは式にすると以下。<br><br>\n",
       "$$p(y|N,q)= \\left(\\begin{array}{cc} N \\\\ y \\end{array} \\right) q^y (1-q)^{N-y} = \\frac{n!}{y!(N-y)!} q^y (1-q)^{N-y}$$\n",
       "<br><br>\n",
       "次に、ロジット関数とは何か。ロジットとは、0から1の値を取るpに対し、以下で表される値を指す。これをpを変数をするロジット関数と呼ぶ。<br><br>\n",
       "$$logit(p)=\\log(\\frac{p}{1-p})=\\log(p)-\\log(1-p)$$ <br><br>\n",
       "また、これはオッズ(成功確率p/失敗確率1-p)の対数であることから、対数オッズとも呼ばれる。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "まず、二項分布は、成功確率qである事象のN回の試行において、y回成功する確率、を表す。これは式にすると以下。<br><br>\n",
    "$$p(y|N,q)= \\left(\\begin{array}{cc} N \\\\ y \\end{array} \\right) q^y (1-q)^{N-y} = \\frac{n!}{y!(N-y)!} q^y (1-q)^{N-y}$$\n",
    "<br><br>\n",
    "次に、ロジット関数とは何か。ロジットとは、0から1の値を取るpに対し、以下で表される値を指す。これをpを変数をするロジット関数と呼ぶ。<br><br>\n",
    "$$logit(p)=\\log(\\frac{p}{1-p})=\\log(p)-\\log(1-p)$$ <br><br>\n",
    "また、これはオッズ(成功確率p/失敗確率1-p)の対数であることから、対数オッズとも呼ばれる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "ロジスティック回帰では、データ(の誤差)が二項分布に従うと仮定すると書いた。これは、あるx点において、そのx点ごとに特有のqをパラメータ\n",
       "とする二項分布に従ってyが生成される、と言い換えることができる。ここで、各x点におけるqを変数xによって以下のように表せると仮定する。(リンク関数)\n",
       "<br><br>\n",
       "$$\\log(\\frac{q_i}{1-q_i})=\\beta_0 + \\beta_1 x_{i1} + \\beta_2 x_{i2} + ... + \\beta_m x_{im}$$<br>\n",
       "(mは説明変数の数)<br><br>\n",
       "ここで、右辺(線形予測子)を$$z_i$$とおくと、式を以下のように変形できる。<br><br>\n",
       "$$\\log(\\frac{q_i}{1-q_i})=z_i$$<br>\n",
       "$$\\frac{q_i}{1-q_i}=exp(z_i)$$<br>\n",
       "$$q_i = \\frac{exp(z_i)}{1+exp(z_i)}$$<br>\n",
       "$$q_i = \\frac{1}{1+exp(-z_i)}$$<br><br>\n",
       "これ(一番下の関数)をlogistic関数(=sigmoid関数)という。logitとlogisticは互いに逆関数の関係である。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "ロジスティック回帰では、データ(の誤差)が二項分布に従うと仮定すると書いた。これは、あるx点において、そのx点ごとに特有のqをパラメータ\n",
    "とする二項分布に従ってyが生成される、と言い換えることができる。ここで、各x点におけるqを変数xによって以下のように表せると仮定する。(リンク関数)\n",
    "<br><br>\n",
    "$$\\log(\\frac{q_i}{1-q_i})=\\beta_0 + \\beta_1 x_{i1} + \\beta_2 x_{i2} + ... + \\beta_m x_{im}$$<br>\n",
    "(mは説明変数の数)<br><br>\n",
    "ここで、右辺(線形予測子)を$$z_i$$とおくと、式を以下のように変形できる。<br><br>\n",
    "$$\\log(\\frac{q_i}{1-q_i})=z_i$$<br>\n",
    "$$\\frac{q_i}{1-q_i}=exp(z_i)$$<br>\n",
    "$$q_i = \\frac{exp(z_i)}{1+exp(z_i)}$$<br>\n",
    "$$q_i = \\frac{1}{1+exp(-z_i)}$$<br><br>\n",
    "これ(一番下の関数)をlogistic関数(=sigmoid関数)という。logitとlogisticは互いに逆関数の関係である。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで、実際のデータについて考えてみる。  \n",
    "https://github.com/tushuhei/statisticalDataModeling/blob/master/data4a.csv  \n",
    "xは体サイズ、fは試肥処理の有無(C:なし、T:試肥処理)、yは8個の種子中何個麦芽したか。(Nはこの8。)  \n",
    "  \n",
    "つまりここでは、体サイズxと試肥処理の有無ごとに、麦芽する確率qを算出したい。このqが算出できれば、将来新しい種子について体サイズ/試肥有無がわかれば、麦芽する(1)/しない(0)を予測分類することができるはずだから。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "      <th>f</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>9.76</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8</td>\n",
       "      <td>6</td>\n",
       "      <td>10.48</td>\n",
       "      <td>C</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   N  y      x  f\n",
       "0  8  1   9.76  C\n",
       "1  8  6  10.48  C"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv('./data4a.csv')\n",
    "df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAD7CAYAAAChScXIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFgBJREFUeJzt3XtspFd5x/Hf49gjv7vOOnFjQGGDJxfCrkgWdqEualUx\naRMIUAEttAVKW9TyB2rdXdGqgl7UuFcVpEqlRfxBG4WL8CaC0kLFpVvEDih/gF2yl4AdoKWzJKXs\nDC0gNlqx2+zTP2a89trzzrzXmTne70eysOc9Puc5r80vs2dn5zF3FwAgDGPDLgAAkByhDQABIbQB\nICCENgAEhNAGgIAQ2gAQkPGiJjIzXjsIACm5u6UZX+gzbXffkR/33Xff0Gtgf+yP/e28jyw4HgGA\ngBDaABAQQjuBWq027BJKxf7Cxv6uLpb1XGXbRGZe1FwAcDUwM/kw/yISAFAuQhsAAkJoA0BACG0A\nCAihDQABIbQBICCENgAEhNAGgIAQ2gAQEEIbAAJCaANAQAhtAAhIotA2s7ea2ZfN7LSZfcjMKmUX\nBgDYru+7/JnZjZIelrTP3S+Y2UOSPuHuH9gyjnf5Q2KtVkuNRkPValWzs7OlfH+3MVsfa7VaOnHi\nhCTp4MGDlx/LMnfcmKmpKZ07dy7xfJL6fl+/ubvNF7d+q9XS8ePHdfbsWd1999264YYb+u5tbW1N\ny8vLmp+f1/79+1Pdl156/cyS3Mci1hukLO/yl6Qdzo2Szki6Xu2ekv8s6e4u4xxIYmnpQY+iGZ+e\nPuRRNONLSw8W/v3dxmx9bGHhiFcq0y7d5tIun5iY8oWFw5nmjhsTRXe6FHkU3Zxovkpl2icmpjyK\nbul8353bvq/f3Jvnm5i41iuV6dhal5Ye9PHxa1261aVdLk34+Pi1Pfe2sHDEpcil212KfGHhcOL7\n0kuvn1mS+5hW3nqL0MnNdC3KEg2SDkv6gaSzkj4YM2ZA20TIms2mR9GMS6dccpdOeRTNeLPZLOz7\n48ZMTl636bHjneDZGCNd79Jk51q6ufuNkWZcOp5w7HRn/PY1+s29urq66Xqzs6futTabTZ+c3Hx9\n+z3ZWu/q6mqX+xb5ww8/XMrP9cqfWfx9TCvv72FRsoR2327sZnadpFdJmpP0fUkfMbM3uPvS1rGL\ni4uXP6/VanScwDaNRkOVSlXnzx/oPHJAExNzajQaif54muT7u40ZG9sr6YeS1h/bLWnvpq8PSKpK\n+k7nWvK5k4xp/99nd8KxT5cUXVHb+vdJ6jn38vLypusrkm7uOs/s7KwajYbMnrnlntwUO16SlpeX\nt42R9urYsWOF/1y3/8zi72NaeX8Ps6rX66rX6/km6Zfqkl4r6e82ff3Lkt7dZdxg/tOEoPFMm2fa\naX6uPNPOcDwiaV7So5ImJZmk90n6zS7jBrZRhG39LHHPnoO5zj57fX+3MVsfW1g43DnTbp/nbj7T\nTjt33JgousOlyCcnq4nmWz/TnpysevsM944eZ9rd594838TElFcq07G1bj/THvfx8Wt77m1h4XAn\nuJ/t3c60i/y5prmPaeWttwhZQjtRj0gzu0/S6yRdlHRC0pvd/eKWMZ5kLkDi1SO95pN49QivHunx\nPUUFLaENAOnQ2BcAdjhCGwACQmgDQEAIbQAICKENAAEhtAEgIIQ2AASE0AaAgBDaABAQQhsAAkJo\nA0BACG0ACAihDQABIbQBICCENnJptVpaWVlRq9Ua2HrHjh3TsWPHtq1ZZC1p5uo1Nu7aoO9bv7XL\nrGeYe92R0nZNiPsQnWuuOoPuZr209KBPTFzb6bBym1cq0107kOetJc1cvcbGXRtmF/AkXeqLrGcU\nOp6PMpXVjT3RRIT2VWXQPfba/Qyv69rz8Mq+iPlqSbOvXmPjrhVZa1F729qHsah6RqUP4yjLEtoc\njyCT9W7WcR3Dy1jvmmuerq3dxcfG9l7uQF5ELWn21Wts3LUia02rW01jY3s797X4egb9O3K1ILSR\nSbVa1YULDUmnO4+c1sWLZy73JCxjvaeeOivpP69Y89KlJzQ/P19YLWn21Wts3LUia02rW02XLj3R\nua/F1zPo35GrRtqn5nEf4njkqjPobtbtM+2pzpn2rV3PtIuoJc1cvcbGXRtmF/AkXerLONMeZsfz\nUaayurEnQWPfq9Ogu1l3655eRi1p5uo1Nu7aMLuAJ+lSX/Z6aKMbOwAEhG7sALDDEdoAEBBCGwAC\nQmgDQEAIbQAICKENAAEhtAEgIIQ2AASE0AaAgBDaABAQQhsAAkJoA0BACG0ACEii0DazaTP7sJmt\nmdlXzOzHyi4MALDdeMJx75L0SXf/eTMbl7SrxJoAADH6vp+2me2RdMLdb+0zjvfT3mHSvnn92tqa\nlpeXNT8/r/379+daT1LixgLDaMSQp7FBnnq77X1zUwgp/r5h9GR5P+0kbcSeJ+mLkh6Q9Iik90qK\nuowrsSkPBm29TdT09KFEbaIWFo64FLl0u0uRLywczrxepTLtExNTXdfeWtfCwpFUdeYVd1+S3q+0\n97XX9y4sHPFKZdql21za5WNjkVcq0wO7F8hPGdqNJQntF0i6KOmFna//WtIfdxk3qH2iZM1m06No\nxqVTLrlLpzyKZrzZbHYdv7q62gnsjfFS5Kurq5nXk653qXnF2tvHHd+2bq8684q7L6urq4nuV9r7\n2vt7t+59sPcCxcgS2knOtJ+Q9Li7/1vn649Ielu3gYuLi5c/r9VqqtVqKZ7zY1Q0Gg1VKlWdP3+g\n88gBTUzMqdFodP0j9/LysqSbJG2Ml/ZqeXk50TFJt/WkqqSGpB+9vLakLeN2b1u3V515xd2X5eXl\nRPcr7X3tvfZuSXu1sffB3gtkU6/XVa/X802SJNklfU7S7Z3P75P0ji5jBvOfJpSOZ9rJ6+SZNvJQ\nGccjvnGuvSLppKSPSpruMmZQ+8QArJ+f7tlzMOGZ9uFOaDw715n2nj0HL59pd1t7a10LC4dT1ZlX\n3H1Jer/S3tde37uwcLhzpn1r50x70iuV6YHdC+SXJbTpxo5YvHqkf528egR5ZHn1CKENAEOSJbT5\nZ+wAEBBCGwACQmgDQEAIbQAICKENAAEhtAEgIIQ2AASE0AaAgBDaABAQQhsAAkJoA0BACG0ACAih\nDQABIbQBICBJ2o2hYMN+/+ey1y96/vX5pqamdO7cuULrzvLe2HH15H0/8SR1Aak6JvT6EJ1rEsnT\njbuI9cruXl70/tbni6I7XYo8im4urO4sndXj6rnnnpfl6kafpC7sPCqr3ViiiQjtvvL0CCxmvXL7\nCBa9v+69I2dcOp677iz9HuPreV+uHplJ6qLX486UJbQ50x6g9Y7a3TpmD2a9+I7d5ayXb/5u80lz\nknbnrjuu1vXO6t32EF/PN3RlZ/SNbvRF1VXW7wjCQ2gPULVa1YULDUmnO4+c1sWLZy73RCx/vScl\nPV7a+kXvr9t80hlJT+auO67W+fn52D3E13OLpCe2PP6E5ufnC6urrN8RBCjtU/O4D3E8kkiebtxF\nrFd29/Ki97dxhnyHS5FPTlYLP9NO01k9rp6XvGT9TDtbN/okdWHnEd3Yw8CrR7LNx6tHsNPQjR0A\nAkI3dgDY4QhtAAgIoQ0AASG0ASAghDYABITQBoCAENoAEBBCGwACQmgDQEAIbQAICKENAAEhtAEg\nIIQ2AAQkcWib2ZiZPWJmHy+zIABAvDTPtI9IWi2rEFw9Wq2WVlZW1Gq1co0paq1RnBuIkyi0zWyv\npJdL+vtyy8FOd/ToQ5qb26d77nmL5ub26ejRhzKNKWqtrMqcG+glURMEM/uwpD+XNC3pd9z9lV3G\n0AQBPbVaLc3N7dP588fVblx7WlF0l86ceexyd5YkY4paq8x9AElkaYIwnmDSV0g66+4nzawmKXaB\nxcXFy5/XajXVarU0tWCHW+80fv789k7j62GXZExRa5W5D6Cber2uer2ea46+z7TN7C8kvVHS/0mK\nJF0r6aPu/itbxvFMGz3xTBu4UpZn2mk7rr9Y0sdjrhXfqhg7TpJO40V1Iy+zqzkd01EEld2N3cxe\nLM60kVOSTuNFdSMvs6s5HdORF93YASAgdGMHgB2O0AaAgBDaABAQQhsAAkJoA0BACG0ACAihDQAB\nIbQBICCENgAEhNAGgIAQ2gAQEEIbAAJCaANAQAhtAAgIoQ0AASG0R0Sr1dLKyoparVam61nnXr+2\ntraWef40a+bZRwhrA6VL2+om7kO0G8tsvXXV9PShrq2r+l3POvf6tSi6xaXIo+jOQlpnxa2ZZx8h\nrA2kpQztxgjtIWs2mx5FMy6dcsldOuVRNOPNZjPR9axzb1w77lK2+dOsubq6mnkfIawNZJEltDke\nGbJGo6FKpap2V29JOqCJiTk1Go1E17POvXFtt6Rs86dZc3l5OfM+QlgbGBRCe8iq1aouXGhIOt15\n5LQuXjyjarWa6HrWuTeuPSkp2/xp1pyfn8+8jxDWBgYm7VPzuA9xPJLZ+nnrnj0He55px13POvf6\ntcnJaudM+45Cz7S3rplnHyGsDaSlDMcjdGMfEa1WS41GQ9VqVbOzs6mvZ517/drU1JTOnTuXaf40\na+bZRwhrA2lk6cZOaAPAkGQJbc60ASAghDYABITQBoCAENoAEBBCGwACQmgDQEAIbQAICKENAAEh\ntAEgIIQ2AASE0AaAgBDaABCQvqFtZnvN7LNm9hUze9TMDg+iMADAdn3f5c/MniHpGe5+0symJH1J\n0qvc/bEt43iXPwBIoZR3+XP3b7v7yc7n5yStSXpmthLDEWrn7lGte1TrAkKT6kzbzKqSni/pi2UU\nMyqOHn1Ic3P7dM89b9Hc3D4dPfrQsEtKZFTrHtW6gBAlboLQORqpS/pTd/9Yl+s74nik1Wppbm6f\nzp8/rnYj2NOKort05sxjI93tZFTrHtW6gFGQ5XhkPOHE45I+IumD3QJ73eLi4uXPa7WaarVamlpG\nwnpH7/Pnt3fuHuWQGdW6R7UuYBjq9brq9XquORI90zazD0j6jrv/do8xPNMeolGte1TrAkZBKX8R\naWY/IemXJP2UmZ0ws0fM7N6sRY662dlZ3X//exRFd2nPnkOKort0//3vGfmAGdW6R7UuIFQ09o0R\naufuUa17VOsCholu7AAQELqxA8AOR2gDQEAIbQAICKENAAEhtAEgIIQ2AASE0AaAgBDaABAQQhsA\nAkJoA0BACG0ACAihDQABIbQBICCENgAEJFG7satFEe/5nHeOot53eljvX91qtXTixAlJ0sGDByWp\nax1Z6uM9uQFJ7l7IR3uqcC0tPehRNOPT04c8imZ8aenBgc9RRA1FzpNl3Upl2qXbXNrlY2ORVyrT\n2+rIUt+w9gSUqZOb6bI27TfEThRwaDebTY+iGZdOueQunfIomvFmszmwOYqooch50tq+7nGXom11\nrK6upq5vWHsCypYltDnT1kbH8HbjWWlzx/BBzVFEDUXOk1aj0dDY2E2b1t0tafPX7TqWl5dT1zes\nPQGjiNCWVK1WdeFCQ9LpziOndfHiGVWr1YHNUUQNRc6TVrVa1aVLj29a90lJm79u1zE/P5+6vmHt\nCRhJaZ+ax30o4OMR940z0z17DuY+0846RxE1FDlPlnXbZ9q3ds60J71Smd5WR5b6hrUnoEzKcDxC\nY99NePVIfrx6BEiObuwAEBC6sQPADkdoA0BACG0ACAihDQABIbQBICCENgAEhNAGgIAQ2gAQEEIb\nAAJCaANAQAhtAAgIoQ0AASG0ASAgiULbzO41s8fM7Gtm9rayiwIAdNc3tM1sTNK7Jb1U0nMlvd7M\n9pVdWB6tVksrKytaW1vTysqKWq1Wz3FZr+cdn1SvebOsWVado7YmsCP165Ig6UWSPrXp67dLeluX\ncSX3eEhmvcNJFN3pUuRRdHPXTif9unun7f5dVrfwXvOG0tWcTupAdyqjG7uk10h676av3yjpb7qM\nG9A243Xr2i3NuHT8iu7d/bp7p+3+XVa38F7zZllzGF3N6aQOxMsS2uNFPmtfXFy8/HmtVlOtVity\n+r7Wu3afP7/RtVuak7T7cvfu2dnZruPSXE+ybq/xefazuQt52jXLqjPrHmgZhqtNvV5XvV7PN0m/\nVFf7eOTTm74e2eMRnmnzTBsIiUo6HrlG0r+r/ZS1IumkpP1dxg1so71snGnf4VLkk5PVnmfWcd29\n03b/LqtbeK95Q+lqTid1oLssoZ2osa+Z3SvpXWq/2uR+d//LLmM8yVyDsN61e2pqSufOnYvt3t2v\nu3fa7t9ldQvvNW8oXc3ppA5sRzd2AAgI3dgBYIcjtAEgIIQ2AASE0AaAgBDaABAQQhsAAkJoA0BA\nCG0ACAihDQABIbQBICCENgAEhNAGgIAQ2gAQEEIbAAJCaCeQuz3QiGN/YWN/VxdCO4Gd/kvD/sLG\n/q4uhDYABITQBoCAFNpurJCJAOAqMrQekQCA8nE8AgABIbQBICC5QtvMbjezE2b2SOd/v29mh4sq\nbhSY2VvN7MtmdtrMPmRmlWHXVCQzO2Jmj3Y+gv/Zmdn9ZnbWzE5veux6MztmZl81s38xs+lh1phH\nzP5e2/kdfcrMDg2zvjxi9vZOM1szs5Nm9g9mtmeYNeYRs78/MbNTnfz8tJk9o988uULb3b/m7gfd\n/ZCkF0h6UtI/5plzlJjZjZJ+S9Ihdz8gaVzS64ZbVXHM7LmSfl3SCyU9X9LPmNktw60qtwckvXTL\nY2+X9Bl3f46kz0r6vYFXVZxu+3tU0s9K+tzgyylUt70dk/Rcd3++pK9r5/3s3unuz3P3g5I+Iem+\nfpMUeTxyt6T/cPfHC5xzFFwjabeZjUvaJelbQ66nSPslfdHdf+juT0n6vKSfG3JNubj7w5K+u+Xh\nV0l6f+fz90t69UCLKlC3/bn7V93965JSvQph1MTs7TPufqnz5Rck7R14YQWJ2d+5TV/ulnRJfRQZ\n2r8o6WiB8w2du39L0l9J+qak/5L0PXf/zHCrKtSXJf1k5/hgl6SXS7ppyDWV4WnuflaS3P3bkp42\n5HqQza9J+tSwiyiamf2ZmX1T0hsk/VG/8YWEtplNSHqlpA8XMd+oMLPr1H6WNifpRklTZvaG4VZV\nHHd/TNI7JP2rpE9KOiHpqaEWNRi8zjUwZvYHki66+9Kwaymau/+huz9L0ofUPo7tqahn2i+T9CV3\nbxU036i4W9I33P1/O8cHH5X040OuqVDu/oC7v9Dda5K+J+lrQy6pDGfN7OmS1PmLnuaQ60EKZvYm\ntf8UuGOeMMVYkvSafoOKCu3Xa4cdjXR8U9KLzGzSzEzST0taG3JNhTKz2c7/Pkvtv8zaCc9kTFee\n735c0ps6n/+qpI8NuqCCbd3f1mshu2JvZnavpN+V9Ep3/+HQqirO1v3dtunaq5UgX3L/i8jOWegZ\nSbe4+w9yTTaCzOw+tV8xclHt44M3u/vF4VZVHDP7vKQZtff3VnevD7eifMxsSVJN0o9IOqv238b/\nk9pHdzep/bv6C+7+vWHVmEfM/r4r6W8l3aD2n5ZOuvvLhlVjVjF7+31JFUn/0xn2BXf/jaEUmFPM\n/l4h6TlqH0uekfQWd//vnvPwz9gBIBz8i0gACAihDQABIbQBICCENgAEhNAGgIAQ2gAQEEIbAAJC\naANAQP4fcoYcAipTwW8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105bf8b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(1,1)\n",
    "chart = ax.scatter(df.x,df.y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "# binarize f \n",
    "df['f_'] = df['f'].map(lambda x: 0 if x=='C' else 1) \n",
    "# add constant to the data\n",
    "x_f_c = sm.add_constant(df[['x','f_']])\n",
    "# define target\n",
    "df['y_failure'] = df['N']-df['y']\n",
    "# modeling with GLM (logit link function is automatically selected for Binomial)\n",
    "model = sm.GLM(df[['y','y_failure']], x_f_c, family=sm.families.Binomial())\n",
    "result = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>  <td>['y', 'y_failure']</td> <th>  No. Observations:  </th>  <td>   100</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>        <th>  Df Residuals:      </th>  <td>    97</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Binomial</td>      <th>  Df Model:          </th>  <td>     2</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>         <td>logit</td>       <th>  Scale:             </th>    <td>1.0</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>        <th>  Log-Likelihood:    </th> <td> -133.11</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Thu, 21 Feb 2019</td>  <th>  Deviance:          </th> <td>  123.03</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>20:26:40</td>      <th>  Pearson chi2:      </th>  <td>  109.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>8</td>         <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th> <th>[95.0% Conf. Int.]</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>  -19.5361</td> <td>    1.414</td> <td>  -13.818</td> <td> 0.000</td> <td>  -22.307   -16.765</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x</th>     <td>    1.9524</td> <td>    0.139</td> <td>   14.059</td> <td> 0.000</td> <td>    1.680     2.225</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>f_</th>    <td>    2.0215</td> <td>    0.231</td> <td>    8.740</td> <td> 0.000</td> <td>    1.568     2.475</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:     ['y', 'y_failure']   No. Observations:                  100\n",
       "Model:                            GLM   Df Residuals:                       97\n",
       "Model Family:                Binomial   Df Model:                            2\n",
       "Link Function:                  logit   Scale:                             1.0\n",
       "Method:                          IRLS   Log-Likelihood:                -133.11\n",
       "Date:                Thu, 21 Feb 2019   Deviance:                       123.03\n",
       "Time:                        20:26:40   Pearson chi2:                     109.\n",
       "No. Iterations:                     8                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [95.0% Conf. Int.]\n",
       "------------------------------------------------------------------------------\n",
       "const        -19.5361      1.414    -13.818      0.000       -22.307   -16.765\n",
       "x              1.9524      0.139     14.059      0.000         1.680     2.225\n",
       "f_             2.0215      0.231      8.740      0.000         1.568     2.475\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAD7CAYAAAChScXIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOX1+PHPk5AhF0KiYVEBmxFwZakgRou2horW7atf\na11qtWhtrT+N0Eq/Vb/qV2y1dWu1tba2SsUtgCIIVlGWEtGiBhUSJbghATdIUBECgWzn98czA0mY\nSWa/Mzfn/XrNyyRz53nOHfBwc+a5zzEiglJKqcyQ5XYASimlIqdJWymlMogmbaWUyiCatJVSKoNo\n0lZKqQyiSVsppTJIj0QNZIzRtYNKKRUlETHRHJ/QK20R8eTj5ptvdj0GPT89Pz0/7z1ioeURpZTK\nIJq0lVIqg2jSjkBJSYnbISSVnl9m0/PrXkysdZW9BjJGEjWWUkp1B8YYxM0PIpVSSiWXJm2llMog\nmrSVUiqDaNJWSqkMoklbKaUyiCZtpZTKIJq0lVIqg2jSVkqpDKJJWymlMogmbaWUyiCatJVSKoNo\n0lZKqQwSUdI2xvzSGPOOMabKGPOEMcaX7MCUUkrtrctd/owxA4FXgMNEpNEYMwt4TkQe7XCc7vKn\nIlZXV0dNTQ1+v5/+/fsn5fWhjun4s7q6OlauXAnA6NGjd/8slrHDHZOXl0d9fX3E4wFdvq6rsUON\nF27+uro6li5dyqZNm5gwYQL9+vXr8tzWrFlDRUUFxcXFHH744VG9L53p7M8skvcxEfOlUiy7/EXS\nDmcgsB7YF9tT8llgQojjRKlIzCwrk0LHkTEFBVLoODKzrCzhrw91TFnZTHGcQikoGCOOUyiTSq+W\nPjk50gtkGEiBzyeTS0tjGjvcMSMdRxyQXo5fHKdQyspmdjpegc8nvpze0tvxiwMy0nH2miN4/JDA\n2B2PaXuePXN6S4HPFzbWmWVl0qdHDxkKkosRQ5b06dGj03MrLZ0s4AgcIuBIaemkveYNd66d6ezP\nzHGGCDjiOCNjGjvS+VItkDeja1EW0UEwCdgGbAIeC3NMik5TZbLa2lopdBypBBGQSpBCx5Ha2tqE\nvT7cMbm5+whUiv3xUskF2TfwvIAsBXHafB/N2F0dk4sjsFQcp7DLY3vSU3IJPUfw+KUghSFira6u\nFscpDJxnreSSGzbW2tpa2Tc3N6rzr66uDiTs4PtYKeDIK6+80mZe+/OO59rln2turlQF5n4b5KDc\nXBnUM18G8LQcwD4ymBfFz0cylGdlVM8C+eLll0WqqkRWrhR54w2R114T+c9/RJYtE1m6VGTxYpEX\nXxR5/nmR+fNF5s4VmT1bZNYs+fqBB+TnOTmyIca/h4kSS9Lushu7MWYf4CygCPgamG2MuVBEyjoe\nO3Xq1N1fl5SUaMcJtZeamhr8Ph+jGhoAGAUU5eRQU1MT0a+nkbw+1DGDs7JYy36B7wB60w/DAKTN\nT+DANkdEOnYkxwwkh4/oTU5OUZfH9qUHTpg4APw+H70bGvCHOKaiogKfz09DwyhgBYPoySh2hoy1\npqaGQca0O//BHcYc0aMHn65YQf9DD4WdO1n71FMcT18cNuKwjlx24pDPV7/7Hb8UB+FfODyFQwMF\nzVlw2WWQlwcNDfaxc+eer9t8X7h9O7U7d5IdmHsEsGrnTppopIWf0kI9zVxGC9m0kI007cS5+GLo\n3Ruys6FHD/vfUI8QzzV9/TVnBf68w/05JkN5eTnl5eVxjdFl0gYmAB+JyJcAxpg5wDig06StVCh+\nv5+axkaqsP+jVAHrm5p2114T8fpQx3zS2kqLbAp8NwrYzmaE7e1+Ah+3+T7SsSM55jOagO00Na3v\n8tgvaMaQFXaOmsZGtgM1IWItLi6msTH4jJ+N7KQaOAL4CBjT0MBhS5bAggUM/+QTbty1i61AfuCY\npUALkA20AvO3bSPv0kttcszN5URj2IfPaOBmdtKPBhppYDNjc3KoadrMNtbSQBFb2c5n7OCC8eNh\nwABwHPvIzd3zdZvvv6yvZ/gRR7B4587d5zPecdghPdm5cw5wDvDc7rN1eo5nfUUFTowJtrWujouK\nilja1BTT38NYdbyYveWWW6IfpKtLcaAYeBvIBQwwHbgqxHEp+oVCZbpgLXF0fn5cNe3OXh/qmGB9\nND9/tDhOoVxderXkBWraQwM17UmBmna0Y4c7ZkSg7uzkFnVZ0x6dn7+7pu3kFokDMiJMTbtvbq4U\n9+wp3wH5pc8n1+bkyOozzxS55BL5bNQ3pdJkyybTQxpBPgd5NytLyrOyZP2xx4pccYXIDTeI/PGP\n8toVV8i5WVlyLsgYkMEYGZidLUf16dNJTXtSoERycMiadvD9jbWm3fZ9ffzxmdKz5yjx+c4TuFR6\n9vzvhNe0Y/17mAjEUB6JqEekMeZm4AKgCVgJ/FREmjocI5GMpRTo6pHOxoPA6hGfj+YPPuAbjY0U\nbN4Ma9fufkhNDc29etE8aBA7+/Yld8gQnIMOggMOgAMO4KvcXGp27WLw6NGQnZ32q0d27ID334eK\niq2sWLGNzZv7sm5dLh98AAUFLQwYsIPBg1s5/fTP+cEP+nbr1SPa2Fcpt4nARx/BG2/A6tV7HjU1\nMGgQDBsGQ4fueQwbBkOGQK9ebkceFRHYtAnefXfvx6ZN9rQOO6z945BDoE8ftyNPHk3aSmWCr7+G\nFSvg9dfhtdfsw+eD4mIYORKOOAKGD7cZq2dPt6ONSVMTvP22PbUVK2DNGpucc3L2TsyHHQZ+v/2M\nsLvRpK1Uumlpgepqm72CSbqmBkaPhmOPtY9jjoHBg92ONC7bt8Mrr8DSpfCf/8DKlTYRH3OM/bdo\n+HCbnPv1czvS9KJJWym3idjLygUL4IUXbKLef//2CXrUKHvJmcGam22SXrLEJupVq2DMGBg/Hr79\nbZuo8/PdjjL9adJWyg3bt9vstWABPP+8/dmpp8Ipp8Dxx3vm8rKlBZYtg1mzYM4c+MY34OSTbaIe\nN86uClTRiSVpR7JOWynVUX09PPccPPUULFoERx0Fp51mk/YRR4CJbjuJdNXaCsuX20Q9e7b9peH8\n822VZ8gQt6PrnjRpKxUpEXj5ZXjwQZg/315e/uAH8Pe/Q9++bkeXMCJQUWET9VNPwT77wHnnwUsv\n2c9Glbs0aSvVldpaeOQReOghu8ThZz+De+7xTNkDbKJeudIm6ieftItWzj/fluWHD3c7OtWWJm2l\nQmlthcWL7VX1okVw9tnwz3/aq2uPlD4Adu2CGTPsv0HbtsEFF8Azz9jPSj10mp6iH0Qq1VZDAzz8\nsM1ieXlw+eVw4YVQUOB2ZAn1xRfwwANw//0wYgRcc439UDFLe1mllH4QqVSstm6Fv/0N7r0Xjj7a\nJu7jjvPc5eann8LvfgdlZfaXhxdftPfzqMyh/66q7m3rVrj1Vnt7eFUVLFxoP2Q8/nhPJezNm+FX\nv7IJulcvu5T8n//UhJ2JNGmr7qm1FaZPh0MPtfdXL18OTzzhuSy2bRvccos9ze3b4Z134K677NI9\nlZm0PKK6n7fegtJSe7fI/Pm2HOIxjY22Xn377bZWXVFhf5lQmU+Ttuo+vvwSbrgB5s6F226DSy/1\n5CdvL74IkybZJL14sed+eej2vPc3VqVUXV0dK1asoK6uLmXzLVy4kIULF+41Z9hYWlrs0r3DD7fr\nrKurbRusThJ2NOfV2bHhnkvG+7Zxo70J5qqr4A9/sDdnhkrYoeZO5p9jqv+OeF60XRPCPdDONd1O\nqrtZzywr26t7enDOsLG8/rrI2LEi48bZBrARzhPpeXV2bLju5PF2Le+otVXkwQdF+vcXuf56kR07\nwh8bau5Ex9PVfGoPktWNPaKBNGl3K/F2VY9lvn1yc9t1T2/bgbxjLINzc2XHxIkiBxwg8sgjNrMl\n+Lw6O7a2tjZkd/L23dKj71re0bp1IuPHixx9tEhlZdfn1nHu3Nx9EhpPV/MlamyviCVpa3lExWR3\nF/HA9x07hidjvv2yszmI9t3CB2dlUVFRsVcsrzQ2sr2uzq5t+/GPI16+F815dXZsTU0NPp+/XbQ5\nOUW7u6V3/Hks79vjj9vPUE85BV591d7F2NW5dZw7O3sAWVnte9DHGk8k8yVq7O5Mk7aKSdsu4pD8\nbtZ+v59NLS2sC8wVnPOT1laKi4vbxdIM3JKdjUyfHvWdjNGcV2fH+v3+Nl3R7bNNTes7dEvf8/No\n3revv4aLLrKfpS5aBL/+dWRdX0LF1NJSS2vrx8QTTzTzJWrsbi3aS/NwD7Q80u2kupv1zLKyvbqn\nB+d85v775eWsLKnIzpYRPXvGFUs059XZseG6k8fTtfytt0QOOkjk5z8X2b49+nMLNXe8XdSjnU/t\nQbK6sUdC9x7pnlLdzTpU93Reegl++EO2X3wx1WefjX/o0Lhjiea8Ojs23HOxvG8zZtilfPffb1eJ\nxCqSLvWJ5HbH83SmnWtU9yIC991n6wSPPWbvIvGglha4/nrbhGDuXPjmN92OSCWKbhiluo+dO+GK\nK+wm0K++6tk2Kl99ZbdLbW62dzV6aAtvFSP9IFJlnro6OPFEu5nG8uWeTdirV9vVIYcfbu9y1ISt\nQJO2yjTvvWe7mo8fb9useLSb7PPPQ0kJ3HST3S22h/5OrAL0r4LKHK+/DmedZWvYl13mdjRJM22a\n3SJl/nz41rfcjkalG03aKjMsWGBvkpk+HU4/3e1okuZPf7JNc5Yt0ya6KjRN2ir9PfkkXH215y89\n77wT/vEPu4KxqMjtaFS60qSt0tujj8J119lb/7q6TzuD/fa3tgfDSy/BoEFuR6PSmSZtlb7+8Q/4\nzW9gyRK7hMKDROyHjXPnQnm5dpRRXdOkrdLT3/5m266Ul8OwYW5HkxQidt+QRYvsaerNgioSmrRV\n+pk2DX7/e5vJPLoGWwQmT7bLzP/9bygsdDsilSk0aav08thjcPPNsHSppxN2aaltVblkSdQbEapu\nTpO2Sh/PPGPrBUuWwMEHux1N0lx/PaxYYfs35ue7HY3KNJq0VXp45RW4/HK7HvuII9yOJmn++Ee7\ncvHllzVhq9hEdBu7MabAGPOUMWaNMWa1MeaYZAemupF33oFzzrFr3o46yu1okuaJJ+wt6S++CH37\nuh2NylSRXmn/CXheRM41xvQAeiUxJtWdfPwxnHaavQQ96SS3o0mahQvhmmvsh44HHuh2NCqTdbmf\ntjEmH1gpIkO7OE730/aYaDevX7NmDRUVFRQXF3N4JOuqv/wSjj8efvpTuOaadvMBETcWcKMRQzSN\nDd54A0491a7FPv74+JoChDr3tk0hIPz7ptJPLPtpR9JG7JvA68DDwFvAPwAnxHFJbMqjUi3YRmtM\nQUFErcQml5aKA3IIiAMyqbS08wl27BA57jiRKVP2mq9PTo4U+Hwh5+4Y1+TS0qjijFe49yXYVqug\nYEy7tlrr1tmG8HPnSqfHRaLja0tLJ0tOTh+BXgLDJDu7t/h8BTGNrdxBDO3GIknaRwFNwNjA9/cC\nt4Q4LlXnqZKstrZWCh1HKu3qNKkEKXQcqa2tDXl8dXW1OIHjgsc7INXV1aEnaG0V+fGPRc47T6Sl\npd18tSD7dhgrOHfHuJYG5ok0zmS9L9XV1eI4hQKVYp+qFMcplLVr62T4cJF7793z+lDHRRLv3q9d\nKpArsG/gZ7Vtvo5ubOWeWJJ2JDXtT4CPReSNwPezgWtDHTh16tTdX5eUlFBSUhLFNb9KFzU1Nfh9\nPkY1NAAwCijKyaGmpibkr9wVFRUcGDguePzgwM9Dlkn+8hdYtcreWZKV1W6+FcBBHcYKzg20i6s3\n7DVvZ3HGK9z7UlFRgc/np6FhTyQ9ehzExIk9OP5429cx+PqOx+XkFEUU796v7Q0MAPoFItn7nYt0\nbJU65eXllJeXxzdIJJkdeAk4JPD1zcAdIY5JzT9NKumSeqVdXi6y334iH30Ucj6vXGn36HGPjBu3\nS3btav96vdJWbZGM8ojsqWuvAFYBc4CCEMek6jxVCgRrt6Pz8yOqFU8K1LQP7qymvWGDyP77iyxc\n2Ol8eYGadqi5O8Y1KVDTjjTOeIV7X4L15vz80ZKTc7n0779N6ur2fn3b42KtaQdfW1o6SXJy8gI1\n7aGSnd1LfL6CmMZW7oglaWs3dhVWQlePNDTAd74D555r73rsYj7IzNUjzz9fx5Qph7F0aRYjR0b3\n+ljm1tUjmS2W1SOatFXyicCll9rEPXMmmOhWOGWKzz+H4mJbsj/rLLejUZkglqStt7Gr5PvrX+3u\nSK++6tmE3dRkf4m4/HJN2Cq59EpbJdeyZTabLV8OQzu9PyujXXMNfPABzJsHWRFtDqGUXmmrdPPJ\nJ3DBBbZlmIcT9uzZ9m7HN9/UhK2ST6+0VXLs3Gk/ePz+922PR496/3047ji7OeHYsW5HozKNfhCp\n0oMIXHYZbNtmO6l7tI69YwcceyxcdRX8/OduR6MykZZHVHp4+GF4/XX78GjCBvjFL2DkSPvho1Kp\noklbJVZ1tV2H/dJLkJfndjRJM2uW7Yj21lue/ndJpSFN2ipxGhrg/PNtF/Xhw92OJmnWrrU9Hl98\nEfr0cTsa1d1oTVslzpVXwldfQVmZZy8/m5rsB48/+pHtpq5UPLSmrdwzb55dQrFqlWcTNsBvfgP9\n+u3ZuU+pVNOkreL32Wd2+cTTT0NBgdvRJM3y5fDgg57/d0mlOb0VQMWntRUuuQSuuMLWDTxq2za4\n+GJ44AHYf3+3o1Hdmda0VXzuu8+2GX/lFejh3V/cLrvMXl0/9JDbkSgv0Zq2Sq0PPoBbbrF1Aw8n\n7DlzoLzclkWUcpt3/09LY27v/5yQ+VtabFnkppvgkEM6nS9R8efl5VFfX5/Q962rzuqOM4Qrr+zL\n3Ll7lveFiyfqbvQxxKVUVB0TOnugnWsiEm2X80TPl7Du5XfeKXLCCSItLZ3OF+/5Bccb6TjigPRy\n/AnrytJVZ/X8/DGSlbVQzj77nd2vCT7nOEMEHHGckeI4hXLSSacKOAKHCDhSWjop5rji6diuMgvJ\najcW0UCatLsUbe/FRM+XsJ6Kq1eL9OvXrs9jqPniPb9Q4+XiCCyNu/9h5P0et0tu7oDdPSrtc0sF\n2h4zN5Cw9/RnBCd8N/ou4oq1j6TKPLEkbV09kkK7u3kHvu/YaTzZ83XWvTxiTU0wcSLceiscdFCn\n88V7fqHGG0gO0Ht3p/FYhYs12Fl9z7vUC59vEDU1Nbs7ott3su0xX2P7z7fvR19RURFTXO3nHxX3\nuSpv0aSdQn6/n5rGRqoC31cB65uadvdETPZ824GPA/PGPP/tt0NhYchdkhJ9fqHG+4wmYDtNTevj\net/CxVpcXMyuXeuB+t3PBOfy+/00NtZg38ka9ryTBcAntH9nP6G4uDimuOwce8aK91yVx0R7aR7u\ngZZHIhJtl/NEzxdX9/KVK21ZZMOGiOdLVE17RKCm7eQWJbym3THWiRPfkKys16RPn6P2mitYb87N\n9Qdq2iPEcQrl5JODNe2DE1bT1q7q3od2Y88MGbl6pLERjj7a9tWaODGq+RIVfypWj6xbZ0/z2We/\npEePtZ12g9fVIype2gRBJc+NN0JVld1jxKP3cIvAhAnwve/Z3WWVSja9uUYlR0WF3XSjstKzCRvs\nKW7bZn+ZUCpdadJWnWtosOWQP//Z05tufPwx3HCDbWzg4Zs7lQdoeUR17le/gg0bbK9HjxKB006z\n+13deKPb0ajuRMsjKrGWL7ebQb39ttuRJNUjj8DGjXDttW5HolTXNGmr0HbutFvb/fnPdtd/j/rs\nM/uh48KFkJPjdjRKdU1vrlGh3XYbHHYY/OAHbkeSNCLw//6f7d9w5JFuR6NUZPRKW+2tstLu9u/x\n1SIzZ9omvR4u1ysP0qSt2mtutmWR22+HgQPdjiZpNm2CX/wC/vUv6NnT7WiUipyWR1R799xj+zz+\n5CduR5JUpaV2O/Cjj3Y7EqWio1faao8PPoA77rA303i4LDJ7tl0Q8+ijbkeiVPR0nbayWlvhu9+F\ns86CX/7S7WiS5osvYORIm7jHjXM7GtXdxbJOW8sjynroIXv346RJbkeSVFOmwLnnasJWmUvLI8ou\nVg7ew52d7XY0SbNokW3Q+847bkeiVOwivtI2xmQZY94yxsxPZkDKBaWldsHyiBFuR5I027fb9dh/\n+xvk5bkdjVKxi+ZKezJQDeQnKRblhjlzYM0amDEjZVNGsld0ovaTDo7z8MPDGTeuF6eeGvNQSYtR\nqahE0ikB2wBvEVACzA9zTJJ6O6ik+eorkYEDRV5+OWVTRtJpPFHdyIPj9O79I4GN8sADc+INP+Ex\nqu6NZHWuMcY8BdyGbYY3RUTODHGMRDKWSiOXX273If3rX1MyXV1dHUVFh9HQsBTbuLYKxxnP+vXv\n7r5SjeSY2ObagOOMjnqcWM9DqUgkZZc/Y8zpwCYRWWWMKQHCTjB16tTdX5eUlFBSUhJNLCqVysvh\nhRdS+qlcsNN4Q8PencaDyS6SY2Kb6xsxjRPreSgVSnl5OeXl5fEN0tWlOPA7YAPwEfA5tk31oyGO\nS8EvEyohduwQOfhgkXnzUjptbW2tOE6hQKXY7ZoqxXEKpba2NqpjIp2rZ8+xAk1xjRPreSgVCWIo\nj0Tbcf0EtKad+a67TuS881yZOpJO44noRt7SInL44ZskJ+e6pHQ1147pKhFiSdpR3RFpjDkBrWln\ntlWrbOfaykrX2oelYvXIgw/a+4Xmzavj44+Ts8JDV4+oeGk3dtW55mY45hi4+mq7W5JHff45jBoF\n//63vWVdqXSlt7Grzt1zD+y7r23U61HBxgZXXKEJW3mT3sbeXXz4YbfYwe/JJ+1mhbNmuR2JUsmh\n5ZHuQAQmTLAtx6dMcTuapKmrs1fX8+bZKpBS6U7LIyq0adNg61aYPNntSJJq8mS46CJN2MrbtDzi\ndZ9+CtdfD0uW2LsfPWr+fFixwi6KUcrLtDziZSK2qcHo0XDLLW5HkzRbttgNCp94Ak44we1olIpc\nUm5jVxls1iz46CPbpsXDpkyx/zZpwlbdgSZtr6qrs+3G588Hn8/taJJm4UJYvFgbG6juQ8sjXnXh\nhTBwINx9t9uRJM22bXa1yN//bm/yVCrTaHlEWc8+2y0+lbv+etuLWBO26k40aXvNli1w5ZXw+OPQ\nq5fb0STNsmUwd66WRVT3o+URr7noInur+n33uR1J0uzYAUceCXfdZT+AVCpTaXmku3vqKVsWWbnS\n7UiS6uabYcwYTdiqe9Irba/4/HO7Hnv+fCgudjuapKmogDPPhLffBt0NVWU6vY29uxKBn/7U9nz0\ncMLetQt+8hO4915N2Kr70qSdJurq6lixYgV1dXXRP//gg7BxI9x0U9SvDT63Zs2aTuePVrg5uzrP\nztx6KwwdCuefn/q5lUob0ba6CfdA243FbGZZmRQ6jowpKJBCx5GZZWWRP//hhyL9+omsXh312MGW\nWb0dvzggIx0n5PzRCo5bUDCmXSuucD+PxGuviey3n8jnn6d+bqWShWT3iOx0IE3aMamtrZVCx5FK\nW+SQSpBCx9ndJLbT55ubRY47TuSee6Iee09z2qWSS/j5YzmfUE1vq6urY26GW19v+xDPnp36uZVK\npliStpZHXFZTU4Pf52NU4PtRQFFODjU1NV0/f/fd9hb1SZOiHrumpgafzw/0ZhDh54/lfOy4e0bM\nySmioqIi5M8jmed//geOPRbOOSf1cyuVbjRpu8zv91PT2EhV4PsqYH1TE36/v9Pnh9bX26Q9fTpk\nhf5j7Gxsv99PY2MNsJ1PCT9/LOdjx90zYlPTeoqLi0P+vKt5XngBnnsusmXniZ5bqbQU7aV5uAda\nHolZsO48Oj+/05p28PknH3lEZORIkenT4xo7WOd1covEARmR4Jp2fv7okHXljj8PZ/NmkUGDRJYs\nSf3cSqUCMZRHdJ12mqirq7PlDL+f/iHWs7V7/u674f33Yc6ciPo9djZ28Lm8vDzq6+vDzp+o8+nq\nPINE7CqRwYPhj39M7dxKpUos67Q1aWeaV16Bc8+FqipPL1Z+7DG4/XZ4803IzXU7GqWSQ29j97pt\n22DiRLsXqYcT9nvvwTXX2A5pmrCVak+vtDPJ5ZdDS4tt1OtRDQ12pciVV8LPf+52NEoll15pe9lz\nz8GiRZ7fI/uaa+Dww+2/T0qpvWnSzgSffAKXXQZPPgn5+W5HkzRPPmn/XXrrrYg+X1WqW9LySLpr\nbobx4+GUU+CGG9yOJmnWroVvfcuuyx4zxu1olEoNXT3iRdddB6tWwfPPh72JJtPt2gXjxsEll8DV\nV7sdjVKpozVtr/nXv+CJJ2y9wKMJG+DXv4aiIigtdTsSpdKfJu10tWGDrWPPmePp5X1z59q+DVrH\nVioyWh5JR42NcMIJcPbZ9jLUo2pqbM+GZ5+FY45xOxqlUk9r2l4xZYq9w2T+fM+WRRob4TvfsTd3\nTpnidjRKuUNr2l4wZw7Mnu35Ovb//q+t+lxzjduRKJVZNGmnk9Wr7W2ACxZA375uR5M0Tz4JTz8N\nb7yhdWylotXlpZwxZrAx5t/GmNXGmLeNMaF33Ffx2bLF1rD/8AcYO9btaJKmshKuusr+QuHhf5eU\nSpoua9rGmP2B/UVklTEmD3gTOEtE3u1wnNa0Y9XcDGecAYccAn/+s9vRJM3mzfaDx9tugx/+0O1o\nlHJfLDXtLq+0RWSjiKwKfF0PrAEGxRZi5khZ524RmDzZ1gmi3Tg6hHTtOP7pp3V873tbOeOMHZqw\nlYpDVJ90GWP8wJHA68kIJl3MmjGDw4qKuOKkkzisqIhZM2Ykb7L77oOXXoKZM6FHfB8xpDTuKJSV\nzaKoaAFVVSt58MFvMGPGLLdDUipjRbzkL1AaKQd+KyLzQjzvifJIXV0dhxUVsbShgVHYroLjHYd3\n169PfLeT556Dn/0Mli+HOPsVpjTuKOMaOPB+mpuvBRygCscZz/r172r3GNXtJW3JnzGmBzAbeCxU\nwg6aOnXq7q9LSkooKSmJJpa0sLuDeUMD0L47eUKTTFUVXHopzJsXd8KGFMYdpWnTvqal5QpswraR\nBTuha9KqnTE9AAANEUlEQVRW3U15eTnl5eVxjRHRlbYx5lFgs4iEXVWrV9pR2LjR3gJ4xx1wwQUJ\nGTIdr7SXL4czz2ylvv5Edu36EwQi0yttpaykfBBpjDkO+BHwXWPMSmPMW8aYU2INMt3179+fv06b\nxnjHYUx+PuMdh79Om5a4BPPFF3DyyXaX/wQlbEhB3FF67z34/vfhsceyePjhK3Cc8eTnj8FxxjNt\n2l81YSsVI72NPYykdO7esgVOPBFOOgl+//uk3FmSDh3HP/sMjjsO/u//bAUoXeJSKt3o3iPpbOtW\nm6zHjbNL+zx6K+CXX9o9RS66yG4FrpQKT5N2uqqvt51nvvlN+MtfPJuwN22y9wiNH2/L9R49TaUS\nJik1bRWnHTvgv/7Ldqu97z7PZrLVq20X9dNP14StVDLplXYy7dwJZ54J++8P06d7dte+xYvhwgvt\ntikXX+x2NEplDi2PpJPGRrt8Ii8PHn887rsd09VDD8GNN9qd+77zHbejUSqz6H7a6aKpCc4/H3w+\neOwxTybs1la7J/bTT8OyZXavK6VU8nkvm7itudkunWhqsvuP5uS4HVHCNTTAj39s7xF69VXo18/t\niJTqPrxZZHVLc7NdmLxli+0+4/O5HVHCbdpkV4f4fLaWrQlbqdTSpJ0oX31lP3SsrbUtxnNz3Y4o\n4aqr4VvfsqsXH38cevZ0OyKluh9N2olQVQVHHw2HHgr/+hf06uV2RAm3eDGUlMDUqfahS/qUcocm\n7XjNmGFvTb/lFrjnHk/WsKdNgx/9CJ56ytaylVLu0Q8iY9XUBNdea7dWXbzY3u3oMa2tcMMNtjy/\nbJn9RUIp5S5N2rHYtMku6XMcWLECCgvdjijhGhpg4kT4/HNdIaJUOtHySLRef93Wr7/9bVu/9mDC\nrq2F737XVnp0hYhS6UWTdjQefNDuI3LfffDb30J2ttsRJVx1td1D5OSTdYWIUulIyyNthN3zedcu\nKC2F//wHXn650+JuvPtGJ2rf6VjGWbLE7iFy112xf+BYV1fHypUrARg9ejRAyDhiiU/35FYKEJGE\nPOxQmausbKY4TqEUFIwRxymUsrKZ9okNG0SKi0XOOUdk69ZOx5hZViaFjiNjCgqk0HFkZllZVDHE\n+/ouz6UT06aJDBggUl4e05S7583J6SPQS2CYZGf3Fp+vYK84Yokvltcole4CeTO6XBvtC8IOlMFJ\nu7a2VhynUKBSQAQqxXEK5au5c0UOOEDkjjtEWlu7HKPQcaTSDiCVIIWOI7W1tRHHEM/ruzqXcOO0\ntIhcf73I0KEi774b1VR7zZubu4/AvoG5a9t8vSeO6urqqOKL5ZyUyhSxJG2taWN/fff5/NjGs5DL\nwdzc4qP3ZZfBo4/Cr3/d5d0ku7uhB75v2w090hjieX24c2nb/byjjRvh3HPtcr7XXotvSV9NTQ3Z\n2fsBBwXmrmnz9Z44KioqIo4vlnNSyus0aQN+v5/GxhqgitN4jnc4hGEtX7J16VKYMCHiMWoaG6kK\nfF8FrG9qwu/3p+T1bccJnktwpKam9e3G2bXL1q1HjIBhwxKzQsTv99PSsglYF5jb3+brPXEUFxd3\nGV8s56RUtxHtpXm4BxlcHhEReeZP98m8rBz5MKunnOHrE1PNNFiTHp2fH1dNO9bXBwXrv/n5o9vV\nf1tbRZ59VmTYMJEzzhB5//2Yhu903pycvEBNe6hkZ/cSn69grzjCxRfLOSmVyYihPKJNEBobbaPd\nu+5i+89+xpozzqDo0ENjXp2QrqtH3n0XfvELWL/e3m1/yikxD93lvLp6RKnIaOeaaIjAwoU2kw0Z\nYtdeDxnidlQJt2UL/OY3thfDDTfAVVd5cnsUpTKSdq6JhAg89xzceit8/TX8/vdw1lme27aupQX+\n+U+46Sa7Y+zq1TBggNtRKaXi1X2SdmsrPPOMTdbNzbax4TnneO6uxpYW2zDnttugTx9YsAACVQql\nlAd4vzzS0mL3FL3tNntP9k032VvRPdYZvaEBHnkE7r7bXlFfe629wvbYLxBKeYqWR9pqbISZM+F3\nv7ObOt15p/30zWNZ7Msv4YEHbEl+7Fh4+GE4/njPnaZSKsCbSXvTJjjySBg+HO6/325Z56Es1tpq\n9wmZNg1eeMGW5BctsuuulVLe5t3yyNq1MHSo21Ek1Lp1MH26ffTrBz/5Cfzwh57cHVapbkGX/HlM\nayu8955tQlBWBpWVdhe+Sy+1v0gopTKbJu0MJgIffwwVFbYZTkUFvPUW9O9vey58//v2g0Xd31op\n79CknUE2b7bJOZigV6ywZffiYvs4+mj7wWLfvm5HqpRKFk3aaaq+3l41B5PzihV21cdRR+1J0Ecf\nDYMHe+rzUqVUFzRpp4HGRnj77fZljnXrYOTIPQm6uBgOPthzS8WVUlHSpJ1iwQ8K25Y53nnHbmHS\nNkGPGAE+n9vRKqXSjSbtJAp+UNi2Bv3mm3bpXbC8UVxsbxnPy3M7WqVUJkha0jbGnALci22aME1E\n7ghxjCeStoj9kPDTT22SXrlyz5W0SPsr6LFj428eoJTqvpKStI0xWcD7wInAZ8AK4AIRebfDcWmT\ntIP7Lufl5VFfXx92/+WO+zN/+KG9iTIvDwYNgv79Gxk06AtKSnoxYUIBBx7Y+QeFydrvubNxM2Vf\nat0LW6m9xZK0I+lIcyywoM331wHXhjgu6q4NyRDs/jLEccQBGek4IbvAhOp83twssn27fT7a7t+J\n6qTeUWdxxNKhPFlxdkY7qSsVGsnoxg6cA/yjzfcXAX8OcVyKTjO8YEfzpSCFgY7moTqbd9X5PNru\n34nqpB5q3HBxxNKhPFlxxnoOSnV3sSTthG4YNXXq1N1fl5SUUFJSksjhuxTsaN67oQE/hOxs3r9/\n/z2dzxsawj7v8/lpaNi7+3eoX+27Gi+e8wkXBxBVjMmMM9Zz0DKJ6m7Ky8spLy+Pb5Cusjq2PPJC\nm+/TtjyiV9p6pa1UJiFJ5ZFs4EOgCPABq4DDQxyXshPtTLBm68/NFQdkRBc17XCdz6Pt/p2oTuod\ndRZHLB3KkxVnZ7STulKhxZK0o1ny9yf2LPm7PcQxEslYqRDr6pFon4/3+Ejp6hGlvElvrlFKqQwS\nS9LW3S+UUiqDaNJWSqkMoklbKaUyiCZtpZTKIJq0lVIqg2jSVkqpDKJJWymlMogmbaWUyiCatJVS\nKoNo0lZKqQyiSVsppTKIJm2llMogmrSVUiqDaNJWSqkMokk7AnG3B0pzen6ZTc+ve9GkHQGv/6XR\n88tsen7diyZtpZTKIJq0lVIqgyS03VhCBlJKqW7EtR6RSimlkk/LI0oplUE0aSulVAaJK2kbYw4x\nxqw0xrwV+O/XxphJiQouHRhjfmmMeccYU2WMecIY43M7pkQyxkw2xrwdeGT8n50xZpoxZpMxpqrN\nz/Y1xiw0xrxnjHnRGFPgZozxCHN+Pwj8HW0xxoxxM754hDm3O40xa4wxq4wxTxtj8t2MMR5hzu83\nxpjKQP58wRizf1fjxJW0ReR9ERktImOAo4DtwNx4xkwnxpiBwNXAGBEZBfQALnA3qsQxxgwHLgPG\nAkcCZxhjhrgbVdweBr7X4WfXAYtF5FDg38D1KY8qcUKd39vA2cBLqQ8noUKd20JguIgcCXyA9/7s\n7hSRb4rIaOA54OauBklkeWQCsFZEPk7gmOkgG+htjOkB9AI+czmeRDoceF1EdolIC7AM+L7LMcVF\nRF4Bvurw47OARwJfPwL8d0qDSqBQ5yci74nIB0BUqxDSTZhzWywirYFvXwMGpzywBAlzfvVtvu0N\ntNKFRCbt84EZCRzPdSLyGfAHYAPwKbBFRBa7G1VCvQN8O1A+6AWcBhzockzJMEBENgGIyEZggMvx\nqNj8BFjgdhCJZoy51RizAbgQ+L+ujk9I0jbG5ABnAk8lYrx0YYzZB3uVVgQMBPKMMRe6G1XiiMi7\nwB3AIuB5YCXQ4mpQqaHrXDOMMeYGoElEytyOJdFE5EYR+QbwBLYc26lEXWmfCrwpInUJGi9dTAA+\nEpEvA+WDOcA4l2NKKBF5WETGikgJsAV43+WQkmGTMWY/gMAHPbUux6OiYIy5BPtboGcumMIoA87p\n6qBEJe0f4rHSSMAG4FhjTK4xxgAnAmtcjimhjDH9A//9BvbDLC9cyRja13fnA5cEvp4IzEt1QAnW\n8fw6PpfJ2p2bMeYU4H+AM0Vkl2tRJU7H8xvW5rn/JoL8EvcdkYFa6HpgiIhsi2uwNGSMuRm7YqQJ\nWz74qYg0uRtV4hhjlgGF2PP7pYiUuxtRfIwxZUAJ0BfYhP00/hls6e5A7N/V80Rki1sxxiPM+X0F\n3Af0w/62tEpETnUrxliFObf/BXzAF4HDXhORK10JME5hzu904FBsWXI9cIWIfN7pOHobu1JKZQ69\nI1IppTKIJm2llMogmrSVUiqDaNJWSqkMoklbKaUyiCZtpZTKIJq0lVIqg2jSVkqpDPL/AUwfG8Zj\nyxpMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x103a64278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "import math\n",
    "df_c = df.query('f_==0').sort_values(by='x')\n",
    "df_t = df.query('f_==1').sort_values(by='x')\n",
    "\n",
    "# get params\n",
    "b1,b2,b3 = result.params\n",
    "# plot\n",
    "fig, ax = plt.subplots(1,1)\n",
    "raw_c = ax.scatter(df_c.x, df_c.y) # plot original data C\n",
    "raw_t = ax.scatter(df_t.x, df_t.y, c='r') # plot original data T\n",
    "fitted_c = ax.plot(df_c.x, df_c.x.map(lambda x:8*1.0/(1.0+math.exp(-(b1+b2*x))))) # plot the fitted line for C\n",
    "fitted_t = ax.plot(df_t.x, df_t.x.map(lambda x:8*1.0/(1.0+math.exp(-(b1+b2*x+b3)))),c='r') # plot the fitted line for T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Iris dataで試してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VPW5x/HPo2yxAopVqKx1AVygaC2utRGVTStWq0Vb\nbWmxeivV1lsrXgvkauvaWrFavahFLHAtV1sEFEWrURKX4AJBE4RKQ0Ep4FqhqIDP/eMMkMSZ5CRz\nMmcm5/t+vebFzJwnv/Pkx+Q88/udzdwdERFJpl3iTkBEROKjIiAikmAqAiIiCaYiICKSYCoCIiIJ\npiIgIpJgkRQBM7vHzNaZWWWG5eea2ZLUo8zMBkSxXhERyU5UI4GpwLAGlq8Ejnf3LwG/BO6KaL0i\nIpKFNlE04u5lZta7geXP13r5PNA9ivWKiEh24tgnMBaYH8N6RUSknkhGAmGZ2QnAGOC4XK5XRETS\ny1kRMLOBwBRguLu/10CcLmYkItJE7m7N+bkop4Ms9fjsArNewIPAee7+RmMNuXtBPiZNmhR7Dso/\n/jyUf2E+Cjn/bEQyEjCzmUAxsJeZ/QOYBLQD3N2nABOALsDvzcyALe4+OIp1i4hI80V1dNC5jSy/\nALgginWJiEh0dMZwhIqLi+NOISvKP17KP16Fnn9zWbbzSVEzM8+3nERE8pmZ4XmwY1hERAqMioCI\nSIKpCIiIJJiKgIhIgqkIiIgkmIqAiEiCqQiIiCSYioCISIKpCIiIJJiKgIhIgqkIiIgkmIqAiEiC\nqQiIiCSYioCISIKpCIiIJJiKgIhIgqkIiIgkmIqAiEiCqQiIiCSYioCISIKpCIiIJFgkRcDM7jGz\ndWZW2UDMrWa2wswWm9mgKNYrIiLZiWokMBUYlmmhmY0A9nf3A4ELgTsjWq+EsGHDBhYtWsSGDRsy\nxlRXVzNt2jSqq6uzbitMTFhh25oxYwajRo1ixowZWceFbau8vJxJkyZRXl6eMSZsv86bN4+xY8cy\nb968rGKass4o84/yMxYmr7Cfi6jjWiV3j+QB9AYqMyy7E/hWrdfVQNcMsS7RmTnzfi8q6uKdOx/u\nRUVdfObM+z8TM27cpQ5FDn0dinzcuEua3VaYmChzd3fv0eOLqfwPdCjynj37NDsubFsnnzyiTtzQ\noSM+ExO2Xw899LA6bQ0YMKhZMU1ZZ5T5R/kZC5NX2M9F1HH5LLXdbN62u7k/+JmGGi4Cc4Fjar1+\nAjg8Q2yLdFISrV+/3ouKujgscXCHJV5U1MXXr1+/I6aqqir1R7czBoq8qqqqyW2FiYkyd3f36dOn\np81/+vTpTY4L21ZZWVnauLKysib369y5c9PGzZ07t0kxTVlnlPlH+RkLk1fYz0XUcfkumyLQJndj\njvBKSkp2PC8uLqa4uDi2XApZTU0N7dr1YfPmgal3BtK2bW9qamrYe++9AaioqAB6AjtjoAcVFRUc\ndNBBTWorTEyUuQPMmjUL6FEv/+7MmjWLb3/7202KC9vWggUL0sYtWLCAY489Fgjfr7Nnz07b1uzZ\nszn11FNDxzRlnc3P/yvMnv06//73QWzaBJs2QWXlB8DFdeI+/fQaxo/fwl57wSefQGXlLsB99dqa\nw+jRHenWDbZuDR4rV/YEnq0Xt5jTT9+DvfaCbdtg8+aOfPRRJdB9R8xHHy2jf/89aNMGPv002JRv\n2bInmzf/Hei0I27z5n/Qu3cHdtkliHGHTz/di48/fgtovyOuuZ/ZXCotLaW0tDSaxppbPeo/aNp0\n0DI0HdTiNBJoWlwSRgJbtrg/8ECFwwkONam4tQ6TfcSIf/rXv+5+7LHu++33kcM/HbalYrY4rPYv\nfvEjP+ww9+OOcx82zH3kyI98113/5PB2Km6dt2nza7/iio1+/fXuv/2t+4QJax0ucliVivm7w1l+\nxx3/8EcfdX/iCfennnK//fYlDsc7LE/Fve5whP/xjy95dbX766+7v/DC296+/SCHqlTMa96hw8G+\ndOnb/s9/uq9b575+vfuyZRu8Q4f9HF5NxS31Dh16+8qVG/zDD90//NB940b3mpr13qHDvg6ViR0J\nRFkE+gBLMywbCTycen4U8HwD7bRIJyXV9vnOTp0Oa2C+9hKvPQ/b2NxvQ22FiYkyd3f3nj37pPI/\nwBuaxw8TF7atoUNH1IlLP6cerl8HDBhUp6108/1hYtzdf/Sjnzgc7vAjh8t80KC/+tlnux99tHv3\n7u5t2wb/du78msMDDn90uN779p3iU6e6P/SQ+zPPuL/6qvuYMb9w2COSz0XYvgjTr2E/F1HH5bNs\nioAFP58dM5sJFAN7AeuASUC7VGJTUjG3AcOBTcAYd385Q1seRU6y04YNG6ipqaFPnz4Zh7jV1dVU\nVFQwePDgOlMHzWkrTEyUuUNwRM+sWbM4++yz60zdNCcubFvl5eUsWLCAoUOH7phGqS9sv86bN4/Z\ns2dz+umn15niaShm40aorITFi4PHK6/Aa6/BPvt8zJ57rufggz/HwIFd6NULevaEXr3gC1+Atm2j\nzz/Kz1iYvMJ+LqKOy1dmhrtbs3423za4KgIin+UebPAffRRefjnY6K9eDYccAoMG7XwMHAgdO8ad\nreSaioBIK/TJJ/D00zBnTvBo0wZOOQUGDw42+P37B++JZFME9BESySPvvQfz58NDD8GCBcGG/rTT\n4JFH4OCDwZr1Zy6SmUYCIjH7+GP44x9h5kx48UUoLoZRo4Jv/d26xZ2dFAKNBEQK0ObNcNddcNNN\ncOihcOmlcPLJsNtucWcmSaIiIJJjGzfCHXfAzTfDkUfCX/4CRxwRd1aSVCoCIjny/vtw221w660w\nZAg89lhwNI9InFQERFrY22/D5MnBt/9TToFnngl2+IrkA91URqSFuMPdd0O/frBuHVRUwLRpKgCS\nXzQSEGkBH34IF10UnOBVVgYNnCArEiuNBEQiVlkZ7OjdbTd44QUVAMlvKgIiEXEPDvk88USYMCF4\nrsM9Jd9pOkgkArWnfxYu1Ly/FA6NBESytGRJ3ekfFQApJCoCIs3kDlOmwEknwcSJmv6RwqTpIJFm\n+PRTGDs2uNaPpn+kkKkIiDSRO/zkJ/DGG/D88/r2L4VNRUCkia67LrjO/9NPqwBI4VMREGmCu+8O\n5v7Ly2GPPeLORiR7up+ASEgPPRQcBvr009C3b9zZiOyk+wmItLCFC4MdwfPnqwBI66JDREUasXQp\nfPObwZ2/dN1/aW1UBEQasGoVjBwJt9wS3PVLpLVRERDJ4O23Ydgw+NnP4Jxz4s5GpGVEUgTMbLiZ\nLTOz5WZ2RZrlncxsjpktNrOlZva9KNYr0lI2bgxuAHPGGcG9f0Vaq6yPDjKzXYDlwInAW8AiYLS7\nL6sVcyXQyd2vNLPPA68DXd19a5r2dHSQxGrLFvj616F79+CQUGvWMRciuZPN0UFRjAQGAyvcfZW7\nbwHuB0bVi3GgY+p5R+CddAVAJB9MnBhs+P/nf1QApPWL4hDR7sDqWq/XEBSG2m4D5pjZW8DuwLci\nWK9I5F58Ef7wh+CS0G10ALUkQK4+5sOAV9x9iJntDzxuZgPdfWO64JKSkh3Pi4uLKS4uzkmSkmwf\nfwxjxsDNN0PXrnFnI5JZaWkppaWlkbQVxT6Bo4ASdx+eej0ecHe/oVbMPOA6dy9Pvf4rcIW7v5im\nPe0TkFhMnAiLFwdnBmsaSApJ3GcMLwIOMLPewFpgNFD/gLpVwElAuZl1BfoCKyNYt0gkFi+GO+8M\n/lUBkCTJugi4+zYzGwcsINjRfI+7V5vZhcFinwL8ErjXzCpTP/Zzd38323WLRGHLlmAa6MYbYd99\n485GJLd0ATlJvGuugeeeg4cf1ihAClM200EqApJoS5fCkCHw8svQs2fc2Yg0T9znCYgUpK1bg2mg\na69VAZDkUhGQxPr1r2HPPYNLRIsklaaDJJGqquBrX4NFi6BPn7izEcmOpoNEmmDbNvj+9+Hqq1UA\nRFQEJHF++1soKoILL4w7E5H4aTpIEmX5cjjmGKiogP32izsbkWhoOkgkBPdgJ/DEiSoAItupCEhi\nzJkD778PF18cdyYi+UMXy5VE2LYNrroKrrsOdt017mxE8odGApIIM2dC585w6qlxZyKSX7RjWFq9\nTz6Bfv1g2jQ4/vi4sxGJnnYMizRgyhTo318FQCQdjQSkVdu0CQ44AB55BA47LO5sRFqGRgIiGUye\nHIwAVABE0tNIQFqtd98N9gWUl0PfvnFnI9JyNBIQSePGG+Eb31ABEGmIRgLSKq1dC4ceCkuWQI8e\ncWcj0rJ0ZzGRen70I9htt+CeASKtnYqASC1vvAFHHgnLlsHnPx93NiItT/sERGqZOBEuuUQFQCQM\njQSkVamshKFDYcUK6Ngx7mxEckMjAZGUq66CK69UARAJK5IiYGbDzWyZmS03sysyxBSb2Stm9qqZ\nPRXFekVqKy8PRgIXXRR3JiKFI+vpIDPbBVgOnAi8BSwCRrv7sloxnYFngaHu/qaZfd7d387QnqaD\npMncobgYvvc9GDMm7mxEcivu6aDBwAp3X+XuW4D7gVH1Ys4FHnT3NwEyFQCR5lqwANavh/POizsT\nkcISRRHoDqyu9XpN6r3a+gJdzOwpM1tkZvpTlci4w9VXw6RJ0Ea3SRJpklz9ybQBDgeGAJ8DnjOz\n59z9b+mCS0pKdjwvLi6muLg4BylKoVq4MBgFnHVW3JmI5EZpaSmlpaWRtBXFPoGjgBJ3H556PR5w\nd7+hVswVQAd3/+/U67uB+e7+YJr2tE9AmmTECDjjDLjggrgzEYlH3PsEFgEHmFlvM2sHjAbm1It5\nCDjOzHY1s92AI4HqCNYtCffKK8ERQeefH3cmIoUp6+kgd99mZuOABQRF5R53rzazC4PFPsXdl5nZ\nY0AlsA2Y4u5V2a5b5Prr4bLLoH37uDMRKUw6Y1gK1vLlcOyx8Pe/w+67x52NSHzing4SicWNN8LF\nF6sAiGRDIwEpSGvWwMCBwTWC9tor7mxE4qWRgCTOzTcHZwerAIhkRyMBKTjvvAMHHhgcFaS7holo\nJCAJ87vfwZlnqgCIREEjASkoH34I++0Hzz4bjAZERCMBSZApU2DIEBUAkahoJCAF4+OPg1HAvHlw\n2GFxZyOSPzQSkES4777gsFAVAJHoaCQgBWHrVujfH6ZOha9+Ne5sRPKLRgLS6j3wAHTrpgIgEjUV\nAcl77sGF4q68Mu5MRFofFQHJe/PnB4Vg5Mi4MxFpfVQEJO9ddx2MHw/WrBlPEWmIioDktWeegbVr\ndetIkZaiIiB57Ve/CkYBuoG8SMvQIaKStxYtCq4R9Le/Qbt2cWcjkr90iKi0Sr/6FVx+uQqASEvS\nSEDy0tKlMHQorFwJRUVxZyOS3zQSkFbn2mvhpz9VARBpaRoJSN7ZfgP5lSuhY8e4sxHJfxoJSKty\n/fUwbpwKgEguaCQgeWXVKjj88OAG8l26xJ2NSGGIfSRgZsPNbJmZLTezKxqI+4qZbTGzM6JYr7Q+\nN94IF1ygAiCSK1mPBMxsF2A5cCLwFrAIGO3uy9LEPQ5sBv7g7n/O0J5GAgm1di0ccghUV0PXrnFn\nI1I44h4JDAZWuPsqd98C3A+MShP3Y+ABYH0E65RW6Oab4TvfUQEQyaUoTsbvDqyu9XoNQWHYwcz2\nBU539xPMrM4yEYB33oF77oElS+LORCRZcnVFlluA2vsKGhy2lJSU7HheXFxMcXFxiyQl+WPy5OAS\nET17xp2JSP4rLS2ltLQ0krai2CdwFFDi7sNTr8cD7u431IpZuf0p8HlgE/BDd5+Tpj3tE0iYDz6A\n/feH55+HAw6IOxuRwpPNPoEoRgKLgAPMrDewFhgNnFM7wN332/7czKYCc9MVAEmm3/8ehg9XARCJ\nQ9ZFwN23mdk4YAHBjuZ73L3azC4MFvuU+j+S7Tql9di0CW65BZ58Mu5MRJJJJ4tJrG65BRYuhAcf\njDsTkcKVzXSQioDE5uOPg30BDz0EX/5y3NmIFK64zxMQaZZ774UBA1QAROKkkYDEYuNG6NcP/vxn\nOPLIuLMRKWwaCUjB+c1v4PjjVQBE4qaRgOTcW28F00AvvQR9+sSdjUjh045hKShjxwZXCb3xxrgz\nEWkd4j5ZTCS0ykqYOxdefz3uTEQEtE9Acuzyy+EXv4A99og7ExEBFQHJoUcfhZoauOiiuDMRke1U\nBCQntm6Fn/0s2A/Qtm3c2YjIdioCkhNTp8Jee8Fpp8WdiYjUpqODpMV9+GFwYticOXDEEXFnI9L6\n6GQxyWs33QRDhqgAiOQjjQSkRa1ZA1/6ErzyCvTqFXc2Iq2TThaTvDVmDHTrBtddF3cmIq2XThaT\nvLR4McyfD8uXx52JiGSifQLSItzhP/8TJk2CTp3izkZEMlERkBbxyCPBheIuuCDuTESkISoCErmt\nW4PLQ9x0E7TRhKNIXlMRkMhdcw307g2nnBJ3JiLSGH1Pk0gtXAhTpgSHhFqzjlUQkVzSSEAi8/77\n8J3vwF13BYeFikj+03kCEgl3GD0a9tkHfve7uLMRSZbYLxthZsPNbJmZLTezK9IsP9fMlqQeZWY2\nIIr1Sv6YNg2qqnS3MJFCk/VIwMx2AZYDJwJvAYuA0e6+rFbMUUC1u39gZsOBEnc/KkN7GgkUmBUr\n4Jhj4Kmn4NBD485GJHniHgkMBla4+yp33wLcD4yqHeDuz7v7B6mXzwPdI1iv5IFPPoFzzw1OClMB\nECk8URSB7sDqWq/X0PBGfiwwP4L1Sh6YNAm6doWLL447ExFpjpweImpmJwBjgOMaiispKdnxvLi4\nmOLi4hbNS5rnySfhvvuCawTpcFCR3CktLaW0tDSStqLYJ3AUwRz/8NTr8YC7+w314gYCDwLD3f2N\nBtrTPoEC8M47MGgQ3HMPDB0adzYiyRbrpaTNbFfgdYIdw2uBCuAcd6+uFdML+Ctwnrs/30h7KgJ5\nzh3OOAP22w9+85u4sxGRWC8l7e7bzGwcsIBgH8M97l5tZhcGi30KMAHoAvzezAzY4u6Ds123xOOu\nu6CmBu6/P+5MRCRbOllMmqS6Gr76VSgrg/79485GRCD+Q0QlIdauhW98I7hLmAqASOugIiChvPUW\nnHACnHee7hEg0pqoCEij3nwTiovhu9+Fq66KOxsRiZKKgDRozZqgAHz/+3DllXFnIyJRUxGQjFav\nDgrABRfA+PFxZyMiLUFFQNL6xz+CAnDRRfDzn8edjYi0FN1ZTD5j1SoYMiS4HtBll8WdjYi0JI0E\npI5Vq4KjgMaNUwEQSQIVAdmhpiaYArr0UvjpT+PORkRyQUVAAHj11WAEcNllQREQkWRQEUi4rVvh\n2muDAlBSAj/+cdwZiUguacdwgi1dCmPGQJcu8NJL0KtX3BmJSK5pJJBAW7bAL38ZHAF00UXw2GMq\nACJJpZFAwlRWBt/+995b3/5FRCOBxNiyBa65Bk48MTj+f/58FQAR0UggERYvDr79d+sGr7wCPXrE\nnZGI5AuNBFqpTz+Fhx+G4cNh2DC45BJ45BEVABGpSyOBVuaDD2DqVLj9dujUKdj4z54NHTrEnZmI\n5CMVgVaiqgpuuy247++wYTBtGhx9NFizbjgnIkmhIlDAtm4NpnhuvRVeew1++MPgzN999407MxEp\nFCoCBcQdVqyAJ56Axx+H0lLo1y+42NtZZ0H79nFnKCKFxtw97hzqMDPPt5zitG4dPPlksNF/4olg\nh+/JJ8NJJwWHe3brFneGIhI3M8PdmzX5qyKQR/71L6iuDub3KyuDjf+qVcGVPU86Kdj49+2reX4R\nqSv2ImBmw4FbCA45vcfdb0gTcyswAtgEfM/dF2doq9UXgXffDTb02zf42x/vvQf9+8PBBwePr30N\nvvIVaKNJOxFpQDZFIOvzBMxsF+A2YBhwCHCOmfWvFzMC2N/dDwQuBO7Mdr1R2LBhA4sWLWLDhg0N\nxpWXlzNp0iTKy8sbbGvhwpeoqHiHsjKYNQsmT4YrroDzzw++yR9yCOy++yd067aZsWPf47nnoHv3\n4NLNCxcGI4EXX4T77oNRo6pZvnwaK1ZUZ1znjBkzGDVqFDNmzGgw/3nz5jF27FjmzZuXMWbixIn0\n69ePiRMnNthWmLgw/QXh+j/s71hdXc20adOors7cX2GFbSvs50ckr7l7Vg/gKGB+rdfjgSvqxdwJ\nfKvW62qga4b2PBdmzrzfO3To6h07jvT27Qf7tdc+6uXl7k884T53rvusWe7Tprn37/87h8sdbnGY\n7N27P+pnnul+4onuRxzhfuCB7p06bXb4yOEjN6vxvn03+Jlnuo8b537tte5Tp7o/9pj7AQec5dDd\n4UCHIh8wYFDa3MaNu9ShyKGvQ5GPG3fJZ2J69PhiKiZoq2fPPmnbOvTQw+rEpVtn27afqxPTrl1R\n2rbCxJ188og6MUOHjsjY/0VFXbxz58O9qKiLz5x5f7N/xzD9FVbYtsLkL5Irqe1m87bhzf1B37nR\nPhOYUuv1d4Bb68XMBY6p9foJ4PAM7bVQN+20fv16Lyrq4rDcg2NuPnKzpf7lL3/iJ5zgPnKk+5ln\nug8bts7hbof1qbg3Hf7Dr7662h97zP2FF9yfffYd79Chv0NlKmaJFxV18fXr19dZ59y5c1MblyU7\n4qDI586dWyeuqqoqbVxVVdWOmOnTp6eNmT59epPXOWHChLQxEyZMqNNWmLiysrK0MWVlZRn6f2dc\n/T4L+zuG6a+wwrYVJn+RXMqmCOTlZSNKSkp2PEpLSyNvv6amhnbt+gAHpt5pT8eO53PHHYt58sng\ncgsPPABHHnk7cAOwdypuX+Bxtm79X4YOhcGDoU2bN2jffjdgQCpmIG3b9qampqbOOmfPng30AAbu\niIPuqfd3qqioAHrWi+uRej8wa9astG0F7zdtnX/605/SxgTv06S4BQsWpI0J3t9pZ//vjKvfZ2F/\nxzD9FVbYtsLkL9KSSktL62wns9Lc6rH9QTAd9Git12Gmg5YR43RQ2G9yYb7Zhm1LI4Gm9b9GAiLh\nEfN00K7A34DeQDtgMXBQvZiRwMO+s2g830B7LdZRtW2f0+3U6bAG53SHDt0+x32AZ5rjDtvWgAGD\n6rSVeZ/AJV57LjzdvHTPnn3qtJVpvjzMOtu1K6oTk2mfQJi4MP3lHq7Pwv6OYforrLBthf0/F8mF\nWItAsH6GA68DK4DxqfcuBH5YK+a2VLFYQob9AZ7DIuAefKOrqKho9BtcWVmZT5w48TPfaJvT1ty5\nc/0HP/jBZ0YA9VVVVfm9997b4Dfa6dOn+2mnnfaZb8fNWeeECRO8b9++nxkBNCcuTH+5h+uzsL9j\nmP4KK2xbYf/PRVpaNkVAJ4uJiBS4WM8TEBGRwqUiICKSYCoCIiIJpiIgIpJgKgIiIgmmIiAikmAq\nAiIiCaYiICKSYCoCIiIJpiIgIpJgKgIiIgmmIiAikmAqAiIiCaYiICKSYCoCIiIJpiIgIpJgKgIi\nIgmmIiAikmAqAiIiCaYiICKSYCoCIiIJpiIgIpJgWRUBM9vTzBaY2etm9piZdU4T08PMnjSz18xs\nqZldks06RUQkOtmOBMYDT7h7P+BJ4Mo0MVuBy9z9EOBo4GIz65/levNSaWlp3ClkRfnHS/nHq9Dz\nb65si8AoYFrq+TTg9PoB7v5Pd1+cer4RqAa6Z7nevFToHyLlHy/lH69Cz7+5si0C+7j7Ogg29sA+\nDQWbWR9gEPBClusVEZEItGkswMweB7rWfgtw4Bdpwr2BdnYHHgAuTY0IREQkZuaecbvd+A+bVQPF\n7r7OzLoBT7n7QWni2gDzgPnuPrmRNpufkIhIQrm7NefnGh0JNGIO8D3gBuC7wEMZ4v4AVDVWAKD5\nv4iIiDRdtiOBLsAsoCewCjjb3d83sy8Ad7n7qWZ2LPAMsJRgusiB/3L3R7POXkREspJVERARkcIW\n2xnDZraLmb1sZnMyLL/VzFaY2WIzG5Tr/BrTUP5m9jUzez+1/GUzS7cTPTZmVmNmS8zsFTOryBCT\nt/3fWP4F0P+dzez/zKw6dRLlkWli8rL/G8s9n/vezPqmPjMvp/79IN3Jq3nc943m35z+z3afQDYu\nBaqATvUXmNkIYH93PzD1IbsTOCrH+TUmY/4pz7j7aTnMpyk+Jdih/166hQXQ/w3mn5LP/T8ZeMTd\nz0odNLFb7YV53v8N5p6Sl33v7suBwyD4EgesAf5SOyaf+z5M/ilN6v9YRgJm1gMYCdydIWQUcB+A\nu78AdDazrhlicy5E/hAcSpuvjIb/7/O6/2k8/+0xecfMOgFfdfepAO6+1d3/VS8sL/s/ZO6Qp31f\nz0nAG+6+ut77edn3aWTKH5rY/3FNB/0WuJzM5xV0B2r/cm+SX2cZN5Y/wNGp4eTDZnZwjvIKy4HH\nzWyRmV2QZnm+939j+UP+9v8XgbfNbGpquD7FzIrqxeRr/4fJHfK372v7FvC/ad7P176vL1P+0MT+\nz3kRMLNTgHWpS0kYhfGtYYeQ+b8E9HL3QcBtwOwcphjGse5+OMFo5mIzOy7uhJqosfzzuf/bAIcD\nt6d+h38TXIOrEITJPZ/7HgAzawucBvxf3Lk0RyP5N7n/4xgJHAucZmYrCSrZCWZ2X72YNwkOO92u\nR+q9fNBo/u6+0d3/nXo+H2hrweG0ecHd16b+3UAwpzi4Xkg+93+j+ed5/68BVrv7i6nXDxBsWGvL\n1/5vNPc87/vtRgAvpT4/9eVr39eWMf/m9H/Oi4C7/5e793L3/YDRwJPufn69sDnA+QBmdhTw/vZr\nFMUtTP615xDNbDDBobjv5jjVtMxsNwsu4YGZfQ4YCrxaLyxv+z9M/vnc/6l+XG1mfVNvnUhwgEFt\nedn/YXLP576v5RwyT6XkZd/XkzH/5vR/nEcH1WFmFwLu7lPc/REzG2lmfwM2AWNiTq9RtfMHvmlm\n/wFsATYTzN/li67AXyy4PEcbYIa7Lyig/m80f/K7/wEuAWakhvUrgTEF1P8N5k6e972Z7UawU/WH\ntd4rlL5vNH+a0f86WUxEJMF0e0kRkQRTERARSTAVARGRBFMREBFJMBUBEZEEUxEQEUkwFQERkQRT\nERARSbAD2l99AAAABklEQVT/BwqkJK4XWBSNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a04d828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import load_iris\n",
    "d = load_iris()\n",
    "y = d.target\n",
    "X = d.data\n",
    "\n",
    "# make the target 2 labels (remove label=2)\n",
    "i = np.where(y<2)\n",
    "y = y[i]\n",
    "X = X[i]\n",
    "\n",
    "# keep only 1 explanatory variable (as all variables makes a perfect separation)\n",
    "x = X[:,0]\n",
    "\n",
    "# modeling\n",
    "model = sm.GLM(y, sm.add_constant(x), family=sm.families.Binomial())\n",
    "result = model.fit()\n",
    "b1,b2 = result.params\n",
    "\n",
    "# sort the x data for drawing line\n",
    "x_plot = sorted(x)\n",
    "# draw\n",
    "fig, ax = plt.subplots(1,1)\n",
    "raw = ax.scatter(x,y)\n",
    "fitted = ax.plot(x_plot, list(map(lambda x: 1.0/(1.0+math.exp(-(b1+b2*x))),x_plot)))"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}