newton method for logistic regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton Method for Logistic Regression\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "Logistic regression is a simple model for binary classification. With a tiny modification, it can be applied to high dimension classification (Softmax regression). Logistic regression is usually in the form below.\n",
    "\n",
    "$$h_\\theta(x)=\\frac {1}{1+e^{-\\theta^Tx}}$$\n",
    "\n",
    "$$J(\\theta)= y log(h_\\theta(x)) + (1-y) log(1-h_\\theta(x))$$\n",
    "\n",
    "In gradient descent\n",
    "\n",
    "$$\\theta:=\\theta-\\alpha\\frac{\\partial J(\\theta)}{\\partial \\theta}=\\theta-\\alpha  (y-h_\\theta(x)) \\cdot x$$\n",
    "\n",
    "In Newton method\n",
    "\n",
    "$$\\theta:=\\theta-H^{-1} \\cdot \\nabla_\\theta {J(\\theta)}=\\theta-\\frac {\\frac{\\partial J(\\theta)}{\\partial \\theta}} {\\frac{\\partial^2 J(\\theta)}{\\partial \\theta^2}}=\\theta- \\frac {(y-h_\\theta(x)) \\cdot x}{h_\\theta(x)(1-h_\\theta(x)) \\cdot x^2}$$\n",
    "\n",
    "NM doesn't only skip the learning parameter $\\alpha$ of gradient descent, but also comes with quicker convergence than GD. NM's mea culpa is the inverse of a Hessian matrix (Fisher information matrix in statistics). When dealing with high dimensional matrix, the computation cost of an inverse can be quite expensive.\n",
    "\n",
    "Reference to Newton method\n",
    "\n",
    "http://thelaziestprogrammer.com/sharrington/math-of-machine-learning/solving-logreg-newtons-method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir('d:/')\n",
    "from sklearn.linear_model import LogisticRegression as LR\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.decomposition import PCA\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#formula of sigmoid function\n",
    "def sigmoid(theta0,theta1,x):\n",
    "    return 1/(1+np.e**(-1*(x*theta1+theta0)))\n",
    "\n",
    "#newton method to fit logit regression\n",
    "def newton_method(y,x):\n",
    "    \n",
    "    #stop is something set by us\n",
    "    #the concept is similar to gradient descent\n",
    "    #when our mean square error is smaller than stop\n",
    "    #we reach convergence therefore stop\n",
    "    #plz refer to my gradient descent file for explanation\n",
    "    # https://github.com/je-suis-tm/machine-learning/blob/master/gradient%20descent.ipynb\n",
    "    stop=0.000000001\n",
    "    \n",
    "    #we use pca to reduce our four dimension iris data into one\n",
    "    #reference to pca\n",
    "    # https://github.com/je-suis-tm/machine-learning/blob/master/principal%20component%20analysis.ipynb\n",
    "    pca_data=pd.DataFrame(PCA(n_components=1).fit_transform(x))[0]\n",
    "    x_train,x_test,y_train,y_test=train_test_split(pca_data,y,test_size=0.3)\n",
    "    \n",
    "    #when we read csv,we may get str value instead of float\n",
    "    #just convert to float from str\n",
    "    y_train=copy.deepcopy(y_train.apply(lambda x: float(x)))\n",
    "    y_test=copy.deepcopy(y_test.apply(lambda x: float(x)))\n",
    "\n",
    "    #set our parameter theta as a vector\n",
    "    #theta is the current parameter vector\n",
    "    #theta_prev is the last parameter vector\n",
    "    #we initialize theta_prev with the first example of training set\n",
    "    theta_prev=np.mat([y_train.iloc[0],x_train.iloc[1]]).T\n",
    "    theta=np.mat([0,0]).T\n",
    "    \n",
    "    #get the fitted value of y\n",
    "    fitted=sigmoid(theta_prev.item(0),theta_prev.item(1),x_train.iloc[0])\n",
    "\n",
    "    for i in range(1,len(y_train)):\n",
    "        \n",
    "        #calculate gradient and hessian matrix\n",
    "        #hessian is just partial derivative matrix of gradient\n",
    "        residual=y_train.iloc[i-1]-fitted\n",
    "        gradient=np.mat([residual,residual*x_train.iloc[i-1]]).T\n",
    "        hessian=np.mat([fitted*(1-fitted)*(x_train.iloc[i-1]**2),\\\n",
    "                       fitted*(1-fitted)*x_train.iloc[i-1],\\\n",
    "                       fitted*(1-fitted)*x_train.iloc[i-1],\\\n",
    "                       fitted*(1-fitted)]).reshape(2,2)\n",
    "\n",
    "        #newton method formula\n",
    "        theta=theta_prev+np.linalg.inv(hessian)*gradient\n",
    "        \n",
    "        fitted=sigmoid(theta.item(0),theta.item(1),x_train.iloc[i])\n",
    "        \n",
    "        if (y_train.iloc[i]-fitted)**2<stop:\n",
    "            print('converged')\n",
    "            break\n",
    "            \n",
    "        theta_prev=copy.deepcopy(theta)\n",
    "    \n",
    "    #calculate fitted y of test and train samples\n",
    "    test=sigmoid(theta.item(0),theta.item(1),x_test)    \n",
    "    curve=sigmoid(theta.item(0),theta.item(1),x_train)\n",
    "    \n",
    "    print('train accuracy: {}%'.format(len(\n",
    "        y_train[curve==y_train])/len(y_train)*100))\n",
    "    print('test accuracy: {}%'.format(len(\n",
    "        y_test[test==y_test])/len(y_test)*100))\n",
    "    print('intercept: %d'%(theta.item(0)))\n",
    "    print('coefficient: %d'%(theta.item(1)))\n",
    "    \n",
    "    #plot logit curve and scatters of train and test samples\n",
    "    plt.plot(sorted(x_train),sorted(curve),c='#009ab0',\n",
    "             lw=2,label='logit',alpha=0.5)\n",
    "    plt.scatter(x_test,y_test,c='#128a08',marker='o',\n",
    "                alpha=0.5,s=85,label='test')\n",
    "    plt.scatter(x_train,y_train,c='#ff0659',marker='*',\n",
    "                alpha=0.2,s=85,label='train')\n",
    "    plt.legend(loc='best')\n",
    "    plt.title('newton method')\n",
    "    plt.show()\n",
    "    \n",
    "    return theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#using sklearn library\n",
    "def logit(y,x):\n",
    "    \n",
    "    pca_data=PCA(n_components=1).fit_transform(x)\n",
    "\n",
    "    x_train,x_test,y_train,y_test=train_test_split(pca_data,y,test_size=0.3)\n",
    "\n",
    "    m=LR()\n",
    "    result=m.fit(x_train,y_train)\n",
    "    print('train accuracy: {}%'.format(m.score(x_train,y_train)*100))\n",
    "    print('test accuracy: {}%'.format(m.score(x_test,y_test)*100))\n",
    "    print('intercept: %s'%result.intercept_.item(0))\n",
    "    print('coefficient: %s'%result.coef_.item(0))\n",
    "\n",
    "    curve=sigmoid(result.intercept_.item(0),result.coef_.item(0),x_train)\n",
    "\n",
    "    plt.plot(sorted(x_train),sorted(curve),\n",
    "             c='#009ab0',lw=2,label='logit',alpha=0.5)\n",
    "    plt.scatter(x_test,y_test,c='#128a08',marker='o',\n",
    "                alpha=0.5,s=85,label='test')\n",
    "    plt.scatter(x_train,y_train,c='#ff0659',marker='*',\n",
    "                alpha=0.2,s=85,label='train')\n",
    "    plt.legend(loc='best')\n",
    "    plt.title('logit')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=pd.read_csv('iris.csv')\n",
    "\n",
    "#make it a binary classification problem\n",
    "df=df[df['type']!='Iris-virginica']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert text type into numerical type\n",
    "df['regressand']=np.select([df['type']=='Iris-setosa',df['type']=='Iris-versicolor'],['0','1'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=pd.concat([df['sepal length'],\n",
    "                 df['sepal width'],\n",
    "                 df['petal length'],\n",
    "                 df['petal width']],axis=1)\n",
    "y=df['regressand']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\ipykernel_launcher.py:2: RuntimeWarning: overflow encountered in double_scalars\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged\n",
      "train accuracy: 95.71428571428572%\n",
      "test accuracy: 96.66666666666667%\n",
      "intercept: -36557931289088960\n",
      "coefficient: 51749885124383608\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4XNWd5//399ZepV2yJVkLkrENdsAsMQ4Emi0hgSyk+YUkwDA/0nSH9NM/ZnqedDodujuZdKanJwmZTKYn6SSkw2QPTSDpkAQCIWC2QMCAAa94k63FtvalVHvd8/vjlqSyLFslqeSqkr6v59Ej3brbuaVbHx2de+65YoxBKaXU0mIVugBKKaXyT8NdKaWWIA13pZRagjTclVJqCdJwV0qpJUjDXSmlliANd6UKRES+KyL/mKdtfU5EfpiPbamlQcNdlTQRMSKyptDlmI2IfFREni10OdTyoeGulFJLkIa7WjQi0iEinxSR10VkRET+TUT8WfPfJyLbRGRYRH4vIhszr/+JiPwya7l9InJ/1nSniJwvIk9nXnpNRMIi8pHM/I9l1hkUkYdEZFXWukZE/lxE9orIkIh8XUTkJOX/nIj8VER+KCJjIvKGiKwTkbtEpDdTjndlLV8pIt8RkSMi0i0i/ygiLhFZD3wTuCRTzuGs3VSLyK8z2/+DiJyZtb23i8hLmffuJRF5e9a8dhF5KrPeb4G6Of+C1NJmjNEv/VqUL6ADeBFYBdQAu4A/z8y7EOgF3ga4gNsyy/uA1cAwTuWjETgEdGfWWw0MAVZm2gBrsvZ5NdCf2b4P+D/A01nzDfAroApoBfqAa09S/s8BMeDdgBv4PnAQ+DvAA3wMOJi1/L8D3wJCwMrMsX88M++jwLPTtv9dYBDYnNn+j4D7MvNqMsf5HzPzbs5M12bmPw98JXOMlwNjwA8L/TvXr+L50pq7Wmz/bIzpMcYMAr8Ezs+8/jHgW8aYPxhj0saY7wFx4GJjzAGcsDofuAJ4FOgWkbMz088YY+yT7O8/APcaY14xxsSBu3BqzG1Zy3zBGDNsjDkMPJlVppk8Y4x51BiTAn4KrMisnwTuA9pEpEpE6oHrgP9ijBk3xvQC/wu4aZb352fGmBcz2/9RVlneC+w1xvzAGJMyxvwE2A28X0RagYuAzxhj4saYp3HeW6UmuQtdALXkHc36OYJTiwc4A7hNRP5T1nxv1vyngCuBNZmfh3GC/ZLM9MmsAl6ZmDDGhEVkAGjC+c9gpjKVnWJ7x7J+jgL9xph01jSZ9Vfh1OaPZLXyWEDnKbZ9qrKswvmPJdshnONYBQwZY8anzWuZZV9qGdFwV4XSCfx3Y8x/P8n8p4D3A+3AP+GE+3/ACfevnWK7PTh/OAAQkRBQC3Tnocyn0onzn0ddphY+3VyHXz3uODJagd8AR3Da6kNZAd86j32oJUybZVShfBv4cxF5mzhCIvJeESnPzH8KuAoIGGO6gGeAa3GC+tWs7RzDaYef8GPgTzIXXH04fxj+YIzpWMyDMcYcAR4D/qeIVIiIJSJnisgVWeVsFhFvjpt8GFgnIreIiDtzsXgD8CtjzCFgK/APIuIVkctw/hAqNUnDXRWEMWYrTrv713AuFO7Dueg4Mf9NIIwT6hhjRoEDwHNZzSLgXPT8XqbHzYeNMb8DPgM8iFPDPZPZ273z5f/FaVraiXNMD+BcEAZ4AtgBHBWR/tk2ZIwZAN4H/BUwAHwKeJ8xZmLdW3AuRg8C/xXnYq9Sk8QY/U9OKaWWGq25K6XUEqThrpRSS5CGu1JKLUEa7koptQQVrJ97XV2daWtrK9TulVKqJL388sv9xpgVsy1XsHBva2tj69athdq9UkqVJBGZfufyjLRZRimlliANd6WUWoI03JVSagkqqoHDkskkXV1dxGKxQhdlUfn9fpqbm/F4PIUuilJqiSqqcO/q6qK8vJy2tjZO8nCckmeMYWBggK6uLtrb2wtdHKXUElVUzTKxWIza2tolG+wAIkJtbe2S/+9EKVVYRRXuwJIO9gnL4RiVUoVVdOGulFJq4Yqqzb0YlJWVEQ6H57Xun/3Zn/GJT3yCDRs28E//9E/87d/+bZ5Lp06b3hFYWXnKRQajA7x+7DWGY0NU+avZWH8eNWPuWdeby/5m3Eegdub1TlXmUy13quneEQbLU+x7dgv9wRjBpkZaKlrpHD1MpPsIlf4qmsqb6PAOEuk+QrCp8cQyTmwTTizfxL4y34df3c626t4TjzdTjuz3oqWilX27t7I1tYfaQzEqL1hPmbectG3jsiwwQtqkCQwmSdT5SNs2gaEEiZoArv7xmcs6w35OWGZnJ2xomfz9PNf5NK8fex0MbKw/j0tb/+jkv6PTqGDjuW/atMlMv0N1165drF+/viDlmbCQcJ/LdorhWNVJhGPw4l7YvBbK/CfMTtkpfr77AZ7qeAIE3JablJ3CGzFcP3oul33wdtwVp3os6+z7O9k+MHBF29XccPaNuCOpqfXg5GXO3v705abvO2s6Zad49sF7ecJ+jU3PpOhttvjN6i460j20uVfxjsMtDA0f4VhqkIH1FbxtrJWDay0Sfpkqo+V2tvnUTpzCv2WqfBP72tBCavshfpN+icRPn+HFqz3Ear2Tx3vVisuofn2IX1VtJxEULLHYP7CPI8f2ce6hcoZWurn0FT/3n3eY8VoXVf4akukkCFQmvKw74GZ7Wwyv5WbdAQ8HVsXZPNzM+PkNJAOuybISjvHsg/fyUMUbJIIy83veG4YfPU3q5rfz4NBj/OD1/8vRsR4ssTACxtg0hFZx68aPcuOGjxz/O5rhXJoPEXnZGLNptuVmrbmLyL04T4TpNcacM8N8Af438B6cB/x+1BjzyvTlSo0xhk996lM88sgjiAh///d/z0c+8hFs2+bOO+/kqaeeor29Hdu2uf3227nxxhu58sor+fKXv8wDDzxANBrl/PPP5y1veQs/+tGPCn04ai4GxmAg7Hyf4QP5890P8MTBx2ipOAOX5Zp8PRSOsnfvVsIvlvG+d96+oP2dbB9pO80TBx8D4EOhd0ytBycvc/b2py83fd9Z07/Z+0uO7N3KuaE6WgcjxJPDUBVmZWM9HAmTONhN3ahNZWUtb3SPIslx1qxZx3C5a6qMG25ytnmoD4TjyzexrwPHeHH74xzue5W3hivYMFrJkTNCk8f70PPf5+z9HlZfdiHjFQF29G3n8GgHNWEXVWE35VFYOeJl/Vgdz1QcIZw4TFO5U7OO9w1SM95CTcRmMNpPzXgrwbEarMFxzOA4q85eP1nWQE+EI3u3svqCFsYrAjO/58MXQM8QW578N74nDzIaG6WhrBERp4XbNjbD8SG+/9p3cFmu439HeQr3XOXSLPNdnEehnewxXtcBazNfbwO+kfm+IJ97efdCNzHzdt96dk7L/exnP2Pbtm289tpr9Pf3c9FFF3H55Zfz3HPP0dHRwRtvvEFvby/r16/n9tuP/yB/4Qtf4Gtf+xrbtm1bjENQ+ZZKw2DWf1mdA1BX7nwPZD3ytKaMweQwT3U8QUvFGbiNRWBo6lnY1f02pqGBN7c9x9C6q6n2V0+uh3sqoGfb31BsiJe3PU7rqlasrGCXtKF81HCW3cSrr/yWawP1pC0/qW0HAeNcQdt2AJJJGI87K5X5YfvhqXlpA+EIPL0TNjTBoX5n+uX90D8CO7shFif8wJP07/49dZ5aanYl6BewjyS4JtHESKWhYriK1MgYZa4qxt3CBfuC9AXCNL8yRvLMctaaVezavYX9I2dR2ZuAdArStrPfpmoYizrNFeMJxrftYcfwbtb1VjFQ5qHy9RjJMSEZsIibBL5uFwdDNud0J2A8TfxwLxceXUFNzA/YtBy2OFrpZt3hFQwnDCC43W76KhKcPdBKr0nTvMdPm6udcCrF6v1BelcGYecIodg4a8ub+MOh3xIaNTQ0NOPut7B9TouGlbAJjtuskSZe7nqKi2JeaPCz7/kthFZZrPSuQsTiWK1N2guShvp4gPDYOE+8/Dib6s+gvDpI2aE+/NPOpePOiUUwa7gbY54WkbZTLPIB4PvGad95QUSqRKQx88DgkvXss89y880343K5qK+v54orruCll17i2Wef5UMf+hCWZdHQ0MBVV11V6KKqhUrbsP8o9Aw5Yeh2QXUIhsbh1YNO88GqaqhczevHXgMBl+XCSthUdccp708RD1oYC2IVHjxHDYe3PEV16KzJ9Y77IM+yv54DW1kRM5hGi+yHxVppM7k/XzLNw9497CpvBGPAAJaA24KxMejP1NLryp0gTdlgG2fZcAySaegbcMoiQP8QdPTASAQsYXxokFSiHSvlwq4Vxn0p3OMN+JIuXGHo8dqEm5IYtwtxW7hTkLBsDoa9pDqDuNIGMx7i//x2K1WVK8DvARvo6IY39jv79bvBhujwMKO0s73ORdLvAmOwDhnEwEgwxYGqC4gEDJ12NWUdBnvwYlJpGPBYhL0pDq22SXqFUNzCm0jjS1vEPDaJUJxDZT4kmaY85safdpP02hx1u4gGXbiiaQYPB/G5ffS63YyVGyqCK3EnbOgEd9KQ9AiB0RT+sM24FeQr5SOky7z0VLWRGlmFlfYwWmZzqCJNCrBsqIlYeMZT7I77+GLqMA0rWnl/wsNbp51LBQ/3HDQBnVnTXZnXTgh3EbkDuAOgtbX1lBvNtYa9WE52LUKfObsE+Txw0Rp48wh09DpBC1Dud0JyfTOsawS3i+HYkNOODKS9FkfWB4l3xqk8miRS7nxYE0EXkdgYbJpaby77O9pq0S1e6r3Hd2bL3p910M3uFj8QpHwkihegqgyqAhBJOBf8BOcPRksIhqMwnPlv4YwVTtD3jUJD9dQflmPDsLoeLCF2eJRUMkko5gQcqQRj/ghecWPSEPEmibji+MRHVcJD0g0Rf4rYCigvDxKtcBOORVk5bqh1WeDN3I1dZYPfCyE/hHwwGqWXUaqPJUhW+nC7DUm/hZUyxEMuIu44uKL4bRBvlHQA+quHIGVTE/VSNWYx7IlTFvcQ9dhI2qajOsZQMElV3MdIWQz/uI035cJv/Ag2kaBQlvIzGkwwHnThDgRIhyOIG8okBR6DN2oTrXKRqnYRtr14D8epPJLE50ljUlEsK0pFIsVQXYLuRhvL5TwZHReM1UCAJHUxG58nSa3bhT/og3DyuHNpseUj3GfqtD1jAhpj7gHuAeeCah72vWguv/xyvvWtb3HbbbcxODjI008/zd133008Hud73/set912G319fWzZsoVbbrnlhPU9Hg/JZFKHGCgVbhdsaIahsNNsYoxTu6wMOq9nVPmrnYtsGcYlDLT58Y+lsdIGY4BUGl9NzXHrzWV/rqY2ErtemHG1if15jhrEFcSEgnxgPMUaccHaZmcb1cClmcrRc7udT+iKWtjW4exnXauzXPQQnNkELgtqqmA8BWtbwGWx5/AuXrE6aI8HKetKMZYKsz/QTxA/R+sN48lxyg9FCXnLaR4NMlhls792HGltodGXoHtdiM7Rw3xkvJ0LB/yQEkDA44cLz3LKJsCOTnaXCSN/2IfdWIntEkYbPFjGsO2qOs78/SEOjhwgbRnOrDgTseAnq17krDct3G4vl3YESaQTuN1eumtipNJJtq7pAoRz+2sJxEOE02HW95Th9/hJmiTl4XL83hC760dod60mvclD45O7YFw4q9yF2AbxQvfa0NQbvxLOuO8A50eqMW4XTx/dQzg+ypbNLs6Y4fc0WjPCBV31fDBwHhdWV2V+t55TnxN5lo9w7wJasqabgZ48bLegbrjhBp5//nnOO+88RIQvfelLNDQ08MEPfpDf/e53nHPOOaxbt463ve1tVFae2M3pjjvuYOPGjVx44YV6QbVURBMwGgURJ3zTaaeGG01Mtr1vrD+PB3beR9pOT17odMdtfOM2Bki5DP6oYY276bj15rK/jevX84DhuH1kk2gSf8RgWdUMRxIEByNOOQfCTlBPlBmmtm/b0JdpqukPO/vqG4O+sFODjsQz02PgdtGcqOJgX5rQQJLQSBq/+BiOe/BjYcQmYQXwjBkqPD4qRy2MgVqPi0C8Ao/YSDSJJ25YY9dD3xC4M1GTSkP3oFOTF6csLaMh7BGbtDtB2ufC9gjGEkJ9CapMFUcjQtotNJZXgYG6eIhQNEFF0s2KcR+JlIXb6ybts4iko1RE3EjAS1XCizst+FIBVkR9+MbdJCyoDASROFQEXTRUVNE7HGdlqhIDuMdTuG1BjPN7Tfmc/548w3HKRmxaaqtBXKyMBXGHxwmN2YyXH/8flm1sPAlYmapkjbsZxuIznkuLLR/h/hBwp4jch3MhdaSU29snui+KCHfffTd33333cfMty+LLX/4yZWVlDAwMsHnzZs4991wAtmzZMrncF7/4Rb74xS+etnKrPBgMO23Tb2mG1jo43A87upzadaAGgJpALVe0XX1cTxb/aAqxDb2tHnYGeniv71IqEmXHrTeX/dXEfCfsY0LaTjPac4hL6s5ja8MKsNMEwmlA4JwWZ3sTZTZMbd9lQfeQE6gbW6FrAJI2rCx3mohe3OdM11fCqmpCR4Zo7m1ncHAfY6srGG0OsJcoq/8wSvWgh47aJD3rPQx7R9ncUcGQO066qYZUXRD3cIqBI4e4uuUKKg65oL4arsh0+316l3OtoSYEaxuBJkJHywjW1/P7tmN4mhqxEOr3x6k9EMFl3Ay0euirsVnjC1F/OMmm5JmMyF7CVpSecjf7V8WJu9O4cbO600XToB9PQx0Vbj8vVXVTGayhLeJhfCiKXRuguymFazSKb1UjXuMmsv8QV9RfyNjqcn4df44N0VXUH07iG0uR8nlJ22lS+w7RVt5G6Mrz4eJ1VN3XR/SXv8XfPUD4rOBxvWVG4sO0xau4sP5CKs4966Tn0mKbtZ+7iPwEuBKoA44B/xXwABhjvpnpCvk14FqcrpB/YoyZ9RFLxdrPPRdXXnklw8PDJBIJPvWpT/HRj350ztsolWNdVvpGnQuS1Vl91IcyAbyiYvKlyT7oh54AoHJUSFg20TLhijMy/aFHYiesN5f9pWqDx+3DLW5SxmkOuqb87bx3/Q186fAxEmNR7trQjs9lTe1voswwtf2JfTkH4HwPR6Es4KzTNzo1nVkvtaeLZ3b8hocrdjgvpWF8fyfurhHMqko8Z7fSNdjB2GAf9b6VtFS2EgkZbJdwWeOlvHft9bhH485F3YljHAo71xpqyqAq5JTp2AipSJRfR57l2SPPYaUNwbiLQH+cvlYPmzZchTGGZzufIhA2VAzZbOcww4cPk/LajDR4sUZiJOJRau0gNVYNg4EYthtiIYvgiE3aBXVRH8l0ku76FOukmabQKnC7uDqwiWsuuQmqQ5PveSBs8NoWIxVOPr43/VauWf9+3GeumjwHHnn8//Kr1x/g5YquyXA3GBpCjdzecjPXb/gg7tqs/+pnOJfmI9d+7noTU4Esp2NdqibuHh2JDVPpr5r5zsxF2kfKtvnHV9/EEuEzF6xb1PGKppdh4g7Vk00v5H041Xs6Uzl29m3njUwPpvbqdsq85di2wbIEjGAbG0ssEHPC6ycr61x+rxN3qL6RuUP13JnuUM0zDfcit5yOVeXfaCLJV97YT5nHzSc3ril0cdRplGu468BhSpWgaMrpAR9w6UdYzUzPDKVKUCQT7sHT0F9alSYNd6VKUDQ9Ee46sKuaWUmH+2B0gC0dT/Dvux9kS8cTDEYHFrS94eFh/uVf/mVe6371q18lEoksaP9K5Wqi5h5wl/RHWC2ikvyzf7LhUB/Yed/xQ43O0US4/8Vf/MWc1/3qV7/KrbfeSjAYnPO6Ss1VNNOdUZtl1MmUZLjnNBzqhpvmvN1Pf/rT7N+/n/PPP59rrrmGlStXcv/99xOPx7nhhhv4h3/4B8bHx/nwhz9MV1cX6XSaz3zmMxw7doyenh6uuuoq6urqePLJJ/N2rErNZLLm7tJwVzMruXAfjA5MDrk6/dZsl+WipeIMnjr0BO9ov2bOfU2/8IUvsH37drZt28Zjjz3GAw88wIsvvogxhuuvv56nn36avr4+Vq1axa9//WsARkZGqKys5Ctf+QpPPvkkdXV1eTtWpU5mqs1dw13NrOQa7LKHXJ3JxOuvH3ttQft57LHHeOyxx7jgggu48MIL2b17N3v37uXcc8/l8ccf52/+5m945plnZhxXRqnFNtXmruGuZlZyNffsIVdPxi1uRmLDC9qPMYa77rqLj3/84yfMe/nll3n44Ye56667eNe73sVnP/vZBe1LqbnSrpBqNiVXc58+5OpMUiZFpb9qztsuLy9nbMwZOe/d734399577+RAYt3d3fT29tLT00MwGOTWW2/lk5/8JK+88soJ6yq12KLa5q5mUXI195mGXM2WttOTy81VbW0tl156Keeccw7XXXcdt9xyC5dccgngPPD6hz/8Ifv27eOv//qvsSwLj8fDN77xDcAZ4ve6666jsbFRL6iqRac1dzWbkhxb5qc77ztpb5nO0UNc3f6uefWWOZ10bBk1X7Yx/LdX9mCAz154FtYiDhqmik+uY8uUXM0d4IazbwSYcTjUq9vfNTlfqaUolk5jAL/L0mBXJ1WS4e623Hxow028o/2aRR9yValiozcwqVyUZLhPqAnUcmXb1YUuhlKnlXaDVLkoud4ySi13kZTTBKk1d3UqGu5KlZho2mmW0W6Q6lQ03JUqMdoNUuViaYR770heNjPfIX/f8573MDy8sDtilcpVVNvcVQ5KP9zDMdh+2Pm+QCcL93RmkKaTefjhh6mqmvsdsUrNh9bcVS5KurcMAANjMBB2vpf5F7Sp7CF/PR4PZWVlNDY2sm3bNnbu3Mkf//Ef09nZSSwW4y//8i+54447AGhra2Pr1q2Ew2Guu+46LrvsMn7/+9/T1NTEL37xCwKBQD6OVClgquYe1DZ3dQqlF+6pNAyGp6Y7B6Cu3Pke8E69XlMGc6zZZA/5u2XLFt773veyfft22tvbAbj33nupqakhGo1y0UUX8cEPfpDa2uP71e/du5ef/OQnfPvb3+bDH/4wDz74ILfeeuu8D1ep6SJpbZZRsyu9cE/bsP8o9Aw5NXW3C6pDMDQOrx50mmdWVUPl6jmH+3SbN2+eDHaAf/7nf+bnP/85AJ2dnezdu/eEcG9vb+f8888H4K1vfSsdHR0LKoNS00W1WUbloPTC3eeBi9bAm0ego9cJdoByP/SPwfpmWNe44GAHCIVCkz9v2bKFxx9/nOeff55gMMiVV15JLHZiO7/P55v82eVyEY1GF1wOpbLpTUwqF6UX7uAE94ZmGAo7zTTGODX6yqDz+jydatjekZERqqurCQaD7N69mxdeeGHe+1FqvowxekFV5aQ0wx0gmoDRKIg4YZ9Og22c17Pb3ucge8jfQCBAfX395Lxrr72Wb37zm2zcuJGzzjqLiy++OF9HolTOkrYhbQxuETxW6Xd2U4undMN9MAwpG97SDK11cLgfdnQ5tflAzbw3++Mf/3jG130+H4888siM8yba1evq6ti+ffvk65/85CfnXQ6lZqK1dpWrnP70i8i1IrJHRPaJyKdnmN8qIk+KyKsi8rqIvCf/RZ3G64bLzoLV9U7NfXW9M+0p3b9XSs1GH4ytcjVrEoqIC/g6cA3QBbwkIg8ZY3ZmLfb3wP3GmG+IyAbgYaBtEco7ZUXFia9Vly3qLpUqNL2YqnKVS819M7DPGHPAGJMA7gM+MG0ZA0ykbSXQM98CFerJUKfTcjhGtTi0G6TKVS7h3gR0Zk13ZV7L9jngVhHpwqm1/6f5FMbv9zMwMLCkw88Yw8DAAH7/wu6mVctTRB+MrXKUSwP1TM/xmp6+NwPfNcb8TxG5BPiBiJxjjLGP25DIHcAdAK2trSdstLm5ma6uLvr6+nIqfKny+/00N8+/y6ZavrTNXeUql3DvAlqypps5sdnlT4FrAYwxz4uIH6gDerMXMsbcA9wDzgOyp+/I4/Ecd0eoUup42uaucpVLs8xLwFoRaRcRL3AT8NC0ZQ4D7wAQkfWAH1ja1W+lCkC7QqpczRruxpgUcCfwKLALp1fMDhH5vIhcn1nsr4CPichrwE+Aj5ql3HCuVIFEtc1d5SinTuHGmIdxLpRmv/bZrJ93Apfmt2hKqem05q5ypfcvK1VC9IKqypWGu1IlRC+oqlxpuCtVItK2IZ62EcDv0o+uOjU9Q5QqEdGsJzCJzHT7iVJTNNyVKhF6MVXNhYa7UiVCu0GqudBwV6pEaM1dzYWGu1IlQrtBqrnQcFeqRGg3SDUXGu5KlQgdy13NhYa7UiVCx3JXc6HhrlSJ0DZ3NRca7kqVCG1zV3Oh4a5UidCukGouNNyVKhFTNzHpx1bNTs8SpUqAMYZo2nkkcdCd02MY1DKn4a5UCYinbWxj8LksXJYOGqZmp+GuVAnQbpBqrjTclSoB2g1SzZWGu1IlQLtBqrnScFeqBGg3SDVXGu5KlQAdy13NlYa7UiVAa+5qrjTclSoBU89P1Y+syo2eKUqVgKmau97ApHKj4a5UCdCx3NVcabgrVQIiOq6MmiM9U5QqATqujJorDXelSkAklQL0gqrKXU5niohcKyJ7RGSfiHz6JMt8WER2isgOEflxfoup1PKVtG2StsElgtfScFe5mfV/PBFxAV8HrgG6gJdE5CFjzM6sZdYCdwGXGmOGRGTlYhVYqeUmmjX0gIiOCKlyk0s1YDOwzxhzwBiTAO4DPjBtmY8BXzfGDAEYY3rzW0ylli+9gUnNRy7h3gR0Zk13ZV7Ltg5YJyLPicgLInLtTBsSkTtEZKuIbO3r65tfiZVaZrQbpJqPXMJ9pv8DzbRpN7AWuBK4GfhXEak6YSVj7jHGbDLGbFqxYsVcy6rUsqRjuav5yCXcu4CWrOlmoGeGZX5hjEkaYw4Ce3DCXim1QDqWu5qPXML9JWCtiLSLiBe4CXho2jL/DlwFICJ1OM00B/JZUKWWq6mx3LWnjMrdrGeLMSYF3Ak8CuwC7jfG7BCRz4vI9ZnFHgUGRGQn8CTw18aYgcUqtFLLSTSlNzCpucvpbDHGPAw8PO21z2b9bIBPZL6UUnk0eQOTDj2g5kDPFqWKXERzPs7qAAARn0lEQVTb3NU8aLgrVeQmmmX0+alqLjTclSpyE80yWnNXc6HhrlSRm7qgquGucqfhrlQRs40hlk4jgF9vYlJzoOGuVBGLpdMYnGC3dNAwNQca7koVMb2BSc2XnjFKFTF9MLaaLw13pYpYVGvuap70jFGqiE3W3PViqpojDXelipg+GFvNl4a7UkVMH4yt5kvPGKWKmN7ApOZLw12pIqZPYVLzpeGuVBHTpzCp+dJwV6qITd3EpOGu5kbDXakiNnUTk4a7mhsNd6WKlDEm6yYmDXc1NxruShWphG2TNgaPJXgs/aiqudEzRqkipQ/GVguh4a5UkdIHY6uF0LNGqSI1NfSAtrerudNwV6pIaTdItRAa7koVKX0wtloIDXelitTEBVUdekDNh4a7UkVKb2BSC6HhrlSRmhhXRtvc1XxouCtVpLTmrhZCw12pIhXVcFcLkFO4i8i1IrJHRPaJyKdPsdyNImJEZFP+iqjU8qRjuauFmDXcRcQFfB24DtgA3CwiG2ZYrhz4z8Af8l1IpZYjHctdLUQuNffNwD5jzAFjTAK4D/jADMv9N+BLQCyP5VNqWUrbhnjaxhLBp8MPqHnI5axpAjqzprsyr00SkQuAFmPMr061IRG5Q0S2isjWvr6+ORdWqeUie1wZESlwaVQpyiXcZzqzzORMEQv4X8BfzbYhY8w9xphNxphNK1asyL2USi0zE+PKaDdINV+5hHsX0JI13Qz0ZE2XA+cAW0SkA7gYeEgvqio1f9oNUi1ULuH+ErBWRNpFxAvcBDw0MdMYM2KMqTPGtBlj2oAXgOuNMVsXpcRKLQPaDVIt1KzhboxJAXcCjwK7gPuNMTtE5PMicv1iF1Cp5Ui7QaqFyukRL8aYh4GHp7322ZMse+XCi6XU8qbdINVCaR8rpYqQjuWuFkrDXakipG3uaqE03JUqQtrmrhZKw12pIqRdIdVCabgrVYR0LHe1UBruShUhrbmrhdJwV6rIGGMmL6hqm7uaLw13pYpMLG1jAJ/LwmXpoGFqfjTclSoy2g1S5YOGu1JFRrtBqnzQcFeqyOjFVJUPGu5KFRntBqnyQcNdqSKjNXeVDxruShUZvaCq8kHDXakioxdUVT5ouCtVZHQsd5UPGu5KFRkdy13lg4a7UkVG29xVPmi4K1VktM1d5YOGu1JFRrtCqnzQcFeqiCRtm5QxuETw6KBhagE03JUqItm1dhENdzV/Gu5KFRG9mKryRcNdqSKi3SBVvmi4K1VEJmvu2lNGLZCGu1JFJKIjQqo80XBXqohoN0iVLxruShWRqLa5qzzRcFeqiGjNXeVLTuEuIteKyB4R2Scin55h/idEZKeIvC4ivxORM/JfVKWWvsmnMLm03qUWZtYzSERcwNeB64ANwM0ismHaYq8Cm4wxG4EHgC/lu6BKLQdTNXd3gUuiSl0u1YPNwD5jzAFjTAK4D/hA9gLGmCeNMZHM5AtAc36LqdTyoDcxqXzJJdybgM6s6a7Mayfzp8AjM80QkTtEZKuIbO3r68u9lEotE1M3MWmzjFqYXM6gmQa4MDMuKHIrsAm4e6b5xph7jDGbjDGbVqxYkXsplVoGbGOIpW0E8OtNTGqBcmnY6wJasqabgZ7pC4nIO4G/A64wxsTzUzyllo+JJhm/y4Wlg4apBcql5v4SsFZE2kXEC9wEPJS9gIhcAHwLuN4Y05v/Yiq19Gk3SJVPs4a7MSYF3Ak8CuwC7jfG7BCRz4vI9ZnF7gbKgJ+KyDYReegkm1NKnURUhx5QeZRTfytjzMPAw9Ne+2zWz+/Mc7mUWna05q7ySS/JK1UktBukyicNd6WKhD4YW+WThrtSRWKizV1r7iofNNyVKhJ6A5PKJz2LlCoS0ZQNaM1d5YeGu1JFIpJKAdrmrvJDw12pIqFdIVU+abgrVSSiaadZRm9iUvmg4a5UETDGaM1d5ZWGu1JFIGHb2MbgtSzcln4s1cLpWaRUEYjog7FVnmm4K1UEdOgBlW8a7koVgamhB/QjqfJDzySlioA+GFvlm4a7UkVgaix3/Uiq/NAzSakioN0gVb5puCtVBHRcGZVvGu5KFQEdy13lm4a7UkVAx3JX+abhrlQR0JuYVL5puCtVBPQmJpVvGu5KFQFtc1f5puGuVIGlbJuEbWOJ4NM7VFWe6JmkVIFFs4YeEJECl0YtFRruShWY3sCkFoOGu1IFNvEEJg13lU8a7koV2OSDsTXcVR4tnSHoekdgZeWJP891/enr9o4431dWwi+3wvs3Qe8Ig+Updu18gYHOAxy0+hnoOsi25D7GKtNcMnYm69dvpjawAs++PnbVDGL3jRBsbORtTW9nqCJF1+hhBiIDrIqVsaJtLRvrz6PmYAQ2tAAwGB3g9WOv0dexlx5/mNpgLe3JFTS0n8XRg3voC8Wo8lc76wVq8/AGqkLRoQfUYsgp3EXkWuB/Ay7gX40xX5g23wd8H3grMAB8xBjTkd+inkI4BtsPw+a1zvTEz2X+ua2/oQV2dk6tG47BS/sBA7Xl8K+Pk6oL8ezu3/KIZxsVLx2hfO8IB0KDNA77qG4co2vVMIHXojz7/IuEUj4uOVjDsTNGaYiG2FMX5Tnf99hxdhzKggTiwjmHQvRtKOdMbws3v9HG2f/5Dn45soUnOx6n99ghVm4f4422cUQszu0I8mbDOJuHmohvaiXmNzyw8z6uaLuaG86+Ebe1dP5WLyfaDVIthlnTQERcwNeBa4Au4CURecgYszNrsT8Fhowxa0TkJuCLwEcWo8AzGhiDgbDzHaZ+zjXcJ9Y/cOz4dQfG4FAfCPD6Iege5o1/+zf2Bg/gr4hR3jFG8xGLMm8FVXEvKZchjc0ZQyEqo27Kom5aBwIELD/pVBJPzBCqtdg+MEDUb9jIWZyRctE3OIg/FaZ/75t8+6efY+faBGOJMAyEaU3XEqOK/vE+GAzT4vbjGouS6h+hZcM5pO00Txx8DIAPbbhpcd5ftaj0gqpaDGKMOfUCIpcAnzPGvDszfReAMeZ/ZC3zaGaZ50XEDRwFVphTbHzTpk1m69at8yt1Kg2DYQC6xiIc2HYA0ja4LEAgnQaXC9Y1Zl4DKgJTP6dtGI0638NR6Bl21ukdBb8Hkmko98KuIzA87qzTPUAinWZ4uJ+xlW4SoyN4o4akSeFPuUm5UiQFvGlABG/KhcsIguAywpgvScqCkWCS7toYXeUR2qnH9rqI+A3+0SSN3np6x46QOr+Fw6MdtA4GifgB22YkNozXuAhFLOIrA8StFCvXnoUpC5AWmx4Z5E/e+nEqfBXze09VweweDtMTiXH9GQ1cWFdV6OKoIiciLxtjNs22XC7/xzcBnVnTXcDbTraMMSYlIiNALdCfW3HnKG3D/qPQM0Sn2DzROwgpG2zj1LJFwG1BIua8Xu6HxmqYqBml0nBkCEYiTtNLLAkhL6RT0DsGx0ay/gCknT8KIS9J0uBdCQakrApTBknLxgIsI4gRZuqlbAAbgxGwLYPYhtpImlGXhW0sYmIY9tkc9bhw+dtx73NTnl7HkMtiFBtPwsZKNhJ3C+FQmohl4U9adB3y4K7wEStzcaTM4oH9e2iral+Ut1wtvgqPNqup/MnlbDpZXs11GUTkDuAOgNbW1hx2fRI+D1y0Bt48QtO+Hi5f0wT9YejLXPysq4DaMhiPQ30VNNdO1donNNZB1wD0DEE8AcMRKC8DcUMgAMeGoSoEPSPQPQAeF2OxccaIERqFmMSJuZK40hZGbLxpC8uWTHgLFoLbFsRASgyWCKP+JCOBFJFAmhFfksNNKVaZGlYPuBh3J4mUwZA3QvtgGWOBFJ0tQk3EwheP4QunEcvF0fIo9SNBImUQctcQqGtkbIUHKz7CmaEolzboxdVSFPK4WF0RKnQx1BKSS7h3AS1Z081Az0mW6co0y1QCg9M3ZIy5B7gHnGaZ+RR4ktsFG5ppHQrTKkBFGYzGwBhY3ejUvEMhePvZJ99Gaz08t9v507SjE5IGAn54S4vTxi4CZzTCd54AYxhKjdPNAIFx2FPVT8NQmrKoRVncTdiXQrAQ2zAUTNA0EiCQdGGAifr8gRXjDCaSvFIbBpPC5/dC+TjrOrzE7BgDQT81o4PUBGo4Wj9EcyqE7YeGaJzx1BiW203ZYIqgv5zD1SnafRbpM+tYSRLvaC9/1FDFlU0rFvS2KqWWhlzC/SVgrYi0A93ATcAt05Z5CLgNeB64EXjiVO3teRNNOE0nImDb0Je5oNofdpplbOMsE/Ceev1EylnHtsGyoHvIaYcXYDzhtO/bUGYZKlMQikFT2o8nZhOICb6Ui0DCPfn/SyDhpirqxWcL2AIIWIaGUT8e20Vzf4qYN4U3WYV/SKgYEYK2m6qVdXRERjhTajkyPIKUgxhhZSRIeDyBcVmkECrSfobCERoqqjgWt4l7nLd6Y/15i/6WK6VKw6zhnmlDvxN4FKcr5L3GmB0i8nlgqzHmIeA7wA9EZB9Ojf30dNsYDDtt6m9pdppduoecgN3Y6lwU3dEFQ2EI1Jx6/ZUVQBNUBmE46tzaFY467fRlfqgqg8ZKPE21xI7twfNCDz78HK0cIVEWd5pfklAb9wHCqCdO2rLxpl0kLBtxufGmYX9thFjAIMk0NRKku8rNqNgMeWPYDRW8flmSpqGzsZ4Ls8a7kl+3dhLwl1M75MXGQ6dnkPAqH9WJcfwNjfjFi3skwT5vD1e3v0v7uyulJs3aW2axLKi3zIS+UaeGXl029TM4gb2iwgn2iZ9PtX7KntrOUNi5oBqOOWG/94gT8I3VMBIhFXTzyv0/YF/PTt6oGaRrrIsdFUfZ0FlOKOFizJ+iccjPkeoYXstDZczLsaoEFQRJWGl6/eM0xMqJVbs5VpOiatxDZWUd1Y3NXL7qj3j3xTfx9M/v5ZWuF3m2oZvYkT7Slk3luAtLXByqi7JGmlhdfgYpF1hpw6aN79R+7kotE7n2lintcC+giTtIu0c7OTxymIPD+9kzsAsLi7aq1by99VJWBhuIJMc5NHKQcGKcFcGVXNZ6BWOJUbpHO+mP9FMXrKOpouWEO02zt5+9XEtFK52jhxmJDVPpr9I7VJVaZjTclVJqCco13HXgMKWUWoI03JVSagnScFdKqSVIw10ppZYgDXellFqCNNyVUmoJ0nBXSqklqGD93EWkDzhUkJ1DHYs1HHFx0+NeXvS4l6YzjDGzjhBYsHAvJBHZmstNAEuNHvfyose9vGmzjFJKLUEa7koptQQt13C/p9AFKBA97uVFj3sZW5Zt7koptdQt15q7UkotaRruSim1BC3LcBeRu0Vkt4i8LiI/F5GqQpfpdBGRD4nIDhGxRWTJdxcTkWtFZI+I7BORTxe6PKeDiNwrIr0isr3QZTmdRKRFRJ4UkV2Zc/wvC12mQlqW4Q78FjjHGLMReBO4q8DlOZ22A/8P8HShC7LYRMQFfB24DtgA3CwiGwpbqtPiu8C1hS5EAaSAvzLGrAcuBv6/ZfL7ntGyDHdjzGPGmFRm8gWguZDlOZ2MMbuMMXsKXY7TZDOwzxhzwBiTAO4DPlDgMi06Y8zTOA+qX1aMMUeMMa9kfh4DdgFNhS1V4SzLcJ/mduCRQhdCLYomoDNruotl/GFfTkSkDbgA+ENhS1I47kIXYLGIyONAwwyz/s4Y84vMMn+H86/cj05n2RZbLse+TMgMr2nf3yVORMqAB4H/YowZLXR5CmXJhrsx5p2nmi8itwHvA95hllhn/9mOfRnpAlqyppuBngKVRZ0GIuLBCfYfGWN+VujyFNKybJYRkWuBvwGuN8ZECl0etWheAtaKSLuIeIGbgIcKXCa1SEREgO8Au4wxXyl0eQptWYY78DWgHPitiGwTkW8WukCni4jcICJdwCXAr0Xk0UKXabFkLprfCTyKc3HtfmPMjsKWavGJyE+A54GzRKRLRP600GU6TS4F/iNwdeZzvU1E3lPoQhWKDj+glFJL0HKtuSul1JKm4a6UUkuQhrtSSi1BGu5KKbUEabgrpdQSpOGulFJLkIa7UkotQf8/hknRfO+7/P0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "matrix([[-3.65579313e+16],\n",
       "        [ 5.17498851e+16]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#our implementation of newton method is not as good as sklearn\n",
    "#sometimes we encounter singular matrix error\n",
    "#the accuracy depends on how we split train and test samples\n",
    "newton_method(y,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train accuracy: 100.0%\n",
      "test accuracy: 100.0%\n",
      "intercept: 0.38365796999300084\n",
      "coefficient: 2.473980898665508\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEICAYAAAB/Dx7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeQXNd59/nv6Zx7pifnGYAACRAEEwiSophFiknMkhX4lvxqVfSuV1Xe8tqyZVve9dZbLslUqbQuvbYs76tXdkmiimIQJVEUIQokmANABCIRaWYwOXZP53DvPfvHHQQSIOLM9HTP86lCzXTfe/s+t8MPd06fc67SWiOEEKJyOMpdgBBCiLMjwS2EEBVGglsIISqMBLcQQlQYCW4hhKgwEtxCCFFhJLhFRVNK9SmlPnWej/EDpdQ356omIeabq9wFCFFuWuv/9cjvSqmbgJ9ordvLV5EQpyZn3EIIUWEkuEVVUEp5lVLfU0oNz/77nlLKe9zyryulRmaXfVUppZVSF8wu+7FS6r8ppYLA80CrUio9+6+1XMckxMeR4BbV4m+Ba4DLgEuB9cDfASil7gD+HPgUcAFw48keQGudAe4EhrXWodl/wwtQuxBnRYJbVIsvAf+P1npcaz0B/APwX2aXfQ74n1rrXVrr7OwyISqWBLeoFq1A/3G3+2fvO7Js4Lhlx/8uRMWR4BbVYhjoOu525+x9ACPA8b1EOk7xODJdplj0JLhFtXgc+DulVINSqh74e+Ans8ueAP6rUmqVUiowu+zjjAF1Sqno/JYrxLmT4BbV4r8Bm4EdwPvAe7P3obV+Hvhn4CXgAPDm7DaFjz6I1nov9n8Ch5RSCelVIhYjJRdSEEuNUmoVsBPwaq2NctcjxNmSM26xJCilHlBKeZRStcC3gV9LaItKJcEtloo/ASaAg4AJ/G/lLUeIcydNJUIIUWHkjFsIISrMvMwOWF9fr7u7u+fjoYUQoipt2bJlUmvdcCbrzktwd3d3s3nz5vl4aCGEqEpKqf7Tr2WTphIhhKgwEtxCCFFhJLiFEKLCLNily0qlEoODg+Tz+YXaZVn4fD7a29txu93lLkUIUaUWLLgHBwcJh8N0d3ejlFqo3S4orTVTU1MMDg7S09NT7nKEEFVqwZpK8vk8dXV1VRvaAEop6urqqv6vCiFEeS1oG3c1h/YRS+EYhRDlJV9OCiFEhVmwNu7FIBQKkU6nz2nbr371q/z5n/85q1ev5h//8R/5m7/5mzmuTiyI8RloPItrJBy//lluO52bYsfYdhL5ODW+WtY2XUrMX3fm9RxZdpr9xg/3sd06dOJ+PrrdSW5Phw0OvPYyk4E8gbYWOiKdjPZ+wHBqkGFfmrpAHSvHwlzwyZtOrP3IY8LJ65vd35H6skMjBNpauNSxjNrO7o+ticYo07kptm/8DduMA2RjLtY2reWTvss/vB0wHTbYs/stet0TTGWnqAvU0VNqYNXqaz5c75m8jic7llMd36kea57NyyRT69at0x8dOblnzx5WrVo15/s6G+cT3GfzOIvhWMVJpPPwzn5YvwJCvrNbH854W8MyeGbvk2zq2wgKXA4XhmWAhhu7b+GBix7G5XCdup4jy1Z3wO6Bk65jWAa/eu/nDP5+I30XOrGC7qP7uaXxBu5NXILrmgvt7T6yLyOZ5rWnfsRGaztXvlpist3F88sHGE2Psn5/GBPNeyvSBIsuPru1nX13NrDu8tuP1X6kxk27sQ/s4g/Xl85jvPUBz3m20v/6Kwx3OegYVAy0a1r7Ldpvu4V7V96Pa3PvsWOb3eaXkW38cvNPuObFIoebSrx9cQZLW9ww1s4V932B+y58EDbt4c3B1/lNeAeuPSO82jJM1mMQLLr45Ggr2ctb+cSq2+x6s8bpX8d0HjbtAhTcuPrYc/bR+87n/XQaSqktWut1Z7Luac+4lVI/Au4BxrXWa863uMVAa83Xv/51nn/+eZRS/N3f/R1/9Ed/hGVZfO1rX2PTpk309PRgWRZf+cpXePjhh7npppv4zne+w5NPPkkul+Oyyy7j4osv5qc//Wm5D0ecqakUTKXtn2fyQTt+fTjjbZ/Z+yQbezfQEenC6XAevd+0TDb2bgDgs6s/f+p6jiw7NPax6zyz90m2vv8i64w6XDpIKuQ5up/N779I88wYn1jRam/3kX397p0nGNm/mUuC9XTGsxRLCahJU2u4CYzmCHsipGijNgWhqRLegTQba46r/UiN/ROgOLG+qRTv7HyRvcZ2VhgN1CQ9RIoG4aSLkFFk8/sv4kkUuSe+6ti2s9v8Wj1PcDzHikQDEbfFpI6gtYVjOsev3vxPvDMGzZtHGZzcg6+jgE7kWNbUxFQY6ifAGc9hTqWOPdfBW0//Ok6l4PCkfcXR4+o54b6Pe3/MQXCfjTNpKvkx8H3gP+dqp//3lr1z9VAfftwrLzqj9Z5++mm2bdvG9u3bmZyc5KqrruKGG27g9ddfp6+vj/fff5/x8XFWrVrFV77ylQ9t+61vfYvvf//7bNu2bT4OQcwlw4Tp4/4yGpiC+rD90+85dn8sBC7nh9c3TNjRDy6H/dM+dYYdh8HttNc/fttZ07kpNvVtPCG0AVzawYVmG1u2v8inQtdQO5A/Vo/bCYmMvV8U9E/a+9vWBy01dg2ZPIT9kMoR9xXZsu1FVmebyNcqakaKePLW0X11JhvYObWTNZt3EslcAPtG7e3fPkDy7R2Unn+ZC9xRoruzODIl6kZT3J+JMVGcYtlUiBlvgaBDUZd2o9yats1xGmtbGNz5Ask9TUSWdcHgNJimvcOtfTAwCYYFCpLDQxw4/D5XpSOYAYOmQwWyYQcN+3NkYi6uztcwtXMTqQvqCb9agvEEmVe3M7J1C2stTWsmhrsEbSMO1r/rxldU+PI+6qezTL3/e5hRNHvCpEbGMMMhmrcp+tsMwlkHDreL0M4Jmt1t7OvfQNoXIRSptV87tP28bu21Azc7e/W6yRQUS5DKw3PvQUed/bpMz4byjn6YyUIyC5HAsdftVO+neXTa4NZav6KU6p7XKhbYa6+9xhe+8AWcTidNTU3ceOONvPvuu7z22mt89rOfxeFw0NzczM0331zuUsX5MC04OArDcfvD53JCbRDiGfuDm85Day1El9nLjl/f64a+ccgb9v0KcDjA5waloWB8eNtZO8a2g+KE0AZwmJq6EYPQgMnh0iZq2644Vs97h+xgMUwIeO3aGsKQK9r1TMzYtdVFIF9kJH+IjhkTTxjyESeB6RJd76SJTJZI17uIt3jQWZh+9V0iH8TtgFJArkhu+06uOmhAIE0hoCiV8iwf93LBgKaoGnBZTjJeH5m4wo2bdETj8Rq4Xp6iaUST9/2OyNo1EAtDTcA+I9112K4vV4L6MHFjnLZxi6Z0EWdR458xyEUceAoaw+Wk5HMQd1vMjG8l7GmA3QPk3trKxcMluh1BTI8brR1o7eDCPQ78ORcmCl84TN7roODwkAo7KehmAhM+nCUH9RlNKgxTNUDewLfZhaaexwP7iTV3YimFCVgOheVUWIeGsIamsVwOzMYolt+D1iWsnQextu5De1xYLTVorw99cBBr1yF0roiO+LHaYmi3i/qwn69a3pO/n+bRnH05qZR6FHgUoLOz85TrnumZ8Xz5uHZ9uahElfG64aoLYN+IHcK1Qfv+sM8+w1rVDitbjn3IPrr+uuUwFIfBSXt5e539wUxkT9x2ViIfP9YG/BGmx8HIqgCGcpPNp+w6jtRTKME9V9rh2jdhB+14EuojkMrBqg4IeGB5MyxrZOQPh3FtdVII2vvPRV042iDR5kEBmZiLqYCbcdNDd8qESzqPnkVOZXsZjjppjbuwnIqZEuS8ReqmFUl3kV3dWfxZuHA6SrLWYqDFYDxYYkXOYusNPjKXLKdxuhlyBQh67WMoFO3nq7UOlGbo0OuMOH3EVRj3FKR7HFimg3idG9NwUHQ4mPSWCLW30hbpoTAxw/77LmJw1EcornDgIufXKA2hrIOs3yIZ1LjyJp6SAo8Tl8NNwciDy0nBo3EbipmIRSagMcwSbdkwLoeL3rY6ul0hSObsWiMB+znPFqFp9vk1LfC5oLUGPE573YjfPqNWCiwLAkFo7QCHspe7HRQsDTUf836aR3MW3FrrHwI/BPvLybl63Plwww038G//9m98+ctfZnp6mldeeYXHHnuMQqHAf/zHf/DlL3+ZiYkJXn75Zb74xS+esL3b7aZUKsmw9krgcsLqdoin7bNZre0PaTRg33+q9bWGjhhMJWd/rzv1tkCNr9b+gvBjaKdiuAO8mdiJ9aztsldK5qAuZP9ZXijZZ/kdMfvMdna/zou7yex9nYCp0RqUpclHnAxdEqTt/QwOU2OaBq7ORpiI2LU7HdBRh7Otju3LDxPc5iE0VsLpgLH6EiZOBkNZxmMmZrRIWz5CKmySjGicBU2+xsX2a1y0rlhO/2AjiX1DzORzpIC030GmrZZ0d4z0gVH669uZzBbwOf34oib5kBNf2v6JAqU1OSNPXW0tOuSDfJ5cIMqodtJUNInmLBwOE6dp4bM0yYBBPmYQGyvSkHGSCzloSXrJ5hJkYg4mwwatYw5cDoP6gkWpmKUt1IpTKa4KdLGmqQnH3iGcSuHobMJhaZxRcKy7AMeWgzj2DuEwHfYfVjVRHNGI/bvlRKFRLi+Oi9pR16xEKYXjrX0oBQ6nw34dT/GemA9LqjvgEQ888ABvvvkml156KUop/umf/onm5mYeeugh/vCHP7BmzRpWrlzJ1VdfTTR6YlefRx99lLVr13LFFVfIl5OVIFe0w1Cp2SYREyxt33982+TJ1rcsmJj9UmsybbePnmLbtU2X8uTun2Na5kmbS0zLxF3QXOBqg2T+xHrA3nfRsPdnWXYTzWTabgef3e/ayCp2ZzVuSuBy4TDts1Nf0sCbsTAx8RXhAt0ME9P2MQR9kM3TXogSi5sEJ0qEEgZ+5SWZddGQ9lAqFZgJWeQcYeLeRkY8Pg66XSR9Dg6Gmtkx6SQUibB5eByyWbs2sOufSUPCD/kijYQwMiVqLIglilhZqE2UMKIKj7ZwmwbFUppPXHQV9dqFdzJPYUWMX4z3suqtKVpSXiZjFkorYtMO9i8vsW9ZibZDBiuS9Uw3KNpTPsZTY0w0Oxg0LDqGnfQ5TfIei3wxx+WqgXQpwz3LWolkDJjMAho6SvZ/YpaGkgWJHExmwOk82pxEyQC3C3we+z7DhKYsmNp+jHThzN9P82BJBfeRLnxKKR577DEee+yxDy13OBx85zvfIRQKMTU1xfr167nkkksAePnll4+u9+1vf5tvf/vbC1a3OE/TaftLs4vbobPe7imwa9A+q/bHTr2+02E3lyhgbSeUzFNuG/PXcWP3LR/bq2Qg2c9dNdcTSYXsP6s/Wo/G3ndjBGizz+QSObu9ezJ9dL+xnJdL6i9hg2c7geXd1E5Y1PfliYwU0IbJrvpJLr3qdiLpNtgxCU1Ruynj3QME8bG65kLGPDsZuKiRw+0B9pei1BxSDNV4eKUnTa0VYPVQiETEYndHhmYjSvOkZrlqoVF7qSmYRENBopd0EnY5CW/tI1hShDw+QoEw3roWXpvZw+6h91idi5CKeIhaJmP1XtypEoWJKTq7Luby9WvtLzm3HobODi6xrsZ64wVGInk2ftKiLunk5tc9lJTmUGCG5kCMxlA3rK9j9/7tuBNBBjwTFJxeXKaT6bDJe21TXF7qobQ/yaUNq4lccqH9H+5I3H5+13TYz/GuQbtJLJGF5hq4YbYL77Pvwv5pWNkK9872ztu0B5KZD79GZ/p+mgen7cetlHocuAmoB8aA/0tr/T9Otc1i7cd9Jm666SYSiQTFYpGvf/3r/PEf//FZP0alHOuSMZG0P7i1oWP3xWfDuSFy6vWP/A7H1j/VthzXj7t/IwAu5cLQdvPJjV238EDd7bg8npPXA/b+DOtYDUeWHbm/IQITSQyH5pmx54/uJ5xx4Jsukq51sv7iW+0+zFNZSOcg5CcZ9TMwMMlgPEn/6DRvxPezT48B4C1oiokU+VKWmYhBKJcjYOTAmSUS9dLQUs/dnst4oOMWvI0xuwdMffjYMcTTdjuvadn/SRgWhkPz3J5n2HLgNZxFC7xusl4TX0Zzdes13L7+c7is2SkiZms0JuK8uPPXPDHya/rNURSKhoQTZwmCdfV85uIHuWPlZ6BQ4jfJTbx6eBPJwUGm05MUXSY5H/haGuiKdnO3Ws/ty+/EtbL941/HsRn7P+fjj2X/iN3NsasBVrR8+PhqgsdeozN9P52hs+nHvaQG4CyUpXSs4uMdGTk5k08Q9dWcfOTkPO4nWSzRm8rSm8rSl8qSKJZO2NbrKJIrjOAlQ3sozMV1HSTzQ4ykBpnMTlIfqKct0nFetZ/L8zCdm+L1w6/y/mwvnUua1nJdxw0nbHfksYeSA3NWb7lIcJfZUjpWsXhorRnK5PlgJs3eRIqJfPFDy31OB+1B/+w/H61BHwHXkmotXdTmdOSkEGLxMiyLQ8ksH8yk+WAmTbp0rEeL1+mgKxSgJxygOxyg2e+V2SurhAS3EBVoJJvnvckEO6aTFMxjIyajHjcX1YS4MBqiKxTA6ZCgrkYS3EJUiIJpsnM6xZbJBMPZYxfraPZ7uagmzEU1IZrkrHpJWLTBfUZTYp6FRCLBz372M/70T//0rLf93ve+x6OPPkogEDjn/QtxLrTWDGfzbJlMsHM6RdGyz679Tidr6yJcWV9Do99b5irFQlt0wf1xU2I+ufvnH54S8ywlEgn+5V/+5ZyD+5FHHpHgFgtGa80HM2leG51mMJM7en9XyM+V9TWsqg3jPjL4RSw5iy64z3hKzLP013/91xw8eJDLLruM2267jcbGRp544gkKhQIPPPAA//AP/0Amk+Fzn/scg4ODmKbJN7/5TcbGxhgeHubmm2+mvr6el156ac6OVYiP0lqzbybNi0MTR3uFBFxOLquLckV9lHqfnF2LRRbcp5oS0+lw0hHpYlP/Rm7tue2sm02+9a1vsXPnTrZt28aGDRt48skneeedd9Bac++99/LKK68wMTFBa2srzz33HAAzMzNEo1G++93v8tJLL1FfXz9nxyrERw2mc/x+aJz+tH2GHfW4+URTLZfX1eBxytm1OGZRBfeppsSEY/fvGNvOTd23nPN+NmzYwIYNG7j88ssBeyj8/v37uf766/mLv/gL/uqv/op77rmH66+//pz3IcSZms4X+cPwBLvi9pwoAZeTG1vqWFdfK71CxEktquA+1ZSYR7iUi5l84rz2o7XmG9/4Bn/yJ39ywrItW7bw29/+lm984xvcfvvt/P3f//157UuIj5M3TF4amWTzRAJTa1xKcW1TjOuaYvgWYGpQUbkWVXCfbkpMAEMbRH01Z/3Y4XCYVMo+o/n0pz/NN7/5Tb70pS8RCoUYGhrC7XZjGAaxWIxHHnmEUCjEj3/84w9tK00lYq7EC0V+emCQyXwRBVxeF+Xm1noiHpkqWJzeogruM5kS88h6Z6uuro7rrruONWvWcOedd/LFL36Ra6+9FrAv/vuTn/yEAwcO8Jd/+Zc4HA7cbjf/+q//CtjTuN555520tLTIl5PivA2kczx+cJCsYdLo9/JQdwtNgYW9ZqGobIturpJf7P75KafEvKXn9nPqVbKQZK4S8XF2TSd5pm8EQ2uWR4J8blkrXqc0i4gKn6vkgYseBjjplJi39Nx+dLkQlURrzetj07w4NAHAuvoa7uxoki8fxTlZdMHtcrj47OrPc2vPbQsyJaYQ8820NM8NjPLe5AwAt7c3cm1jrQxNF+ds0QX3ETF/3Xl1+RNiMcgbJk8cGuJQKotLKR7saWF17blPti8ELOLgFqLSJQolfnpggIl8kaDLyReWt9Me8pe7LFEFJLiFmAdDmRyPHxwiXTJo8Hn40gUd1Hilq5+YGxLcQsyxPfEUT/cNU7I0PeEAf7SsTQbUiDm1+CdAGJ+Zk4c5Mjvg2brrrrtIJM5vpKZYGrTWvDE2zROHhihZmsvrojxyQYeEtphzizu403nYedj+eZ4+LrhN0zzldr/97W+pqTn7kZpiabG05rnDY2wYHEcDt7Y2cG9Xs3T3E/NicTeVTKVgKm3/DJ3fyLLjp3V1u92EQiFaWlrYtm0bu3fv5v7772dgYIB8Ps+f/dmf8eijjwLQ3d3N5s2bSafT3HnnnXzyk5/kjTfeoK2tjWeffRa/X75sWuoKpskvDg1zIJnBpRT3d7ewJiY9R8T8WVzBbZgwnT52e2AK6sP2T7/n2P2xEJzln5/HT+v68ssvc/fdd7Nz5056enoA+NGPfkQsFiOXy3HVVVfx0EMPUVf34X7j+/fv5/HHH+ff//3f+dznPsdTTz3FI488cs6HKyrfTLHEzw4MMpYrEJjtOdIhPUfEPFtcwW1acHAUhuP2GbbLCbVBiGdga6/dZNJaC9FlZx3cH7V+/fqjoQ3wz//8zzzzzDMADAwMsH///hOCu6enh8suuwyAK6+8kr6+vvOqQVS2kWyenx0YJFUyqPd5+OLydmI+z+k3FOI8La7g9rrhqgtg3wj0jduhDRD2wWQKVrXDypbzDm2AYDB49PeXX36ZF198kTfffJNAIMBNN91EPn9iu7rXe+zqI06nk1wud8I6Ymn4IJHiqd4RipZF92zPEb98CSkWyOIKbrBDeXU7xNN204nW9pl4NGDff46On9b1o2ZmZqitrSUQCLB3717eeuutc96PqH5vj0/zuwH7S8hLYxHu7WqRLyHFglp8wQ2QK0IyB0rZQW6aYGn7fv+5/Sl6/LSufr+fpqamo8vuuOMOfvCDH7B27VouvPBCrrnmmrk6ElFFLK15YXCct8fjANzcWs8NzXUy54hYcItuWlcAhqbhvV64uB066+HwJOwahCt7oDU25/XONZnWtfoUTYuneof5YCaNUynu62pmbV203GWJKlLR07oC4HHBJy+E2pB9e1mT3d5tWOWtSyxJOcPkP/cPMJLN43c6+fzyNrrCgXKXJZawxRncDSfpA3skxIVYQKaleeLQECPZPDGvmy9d0EGd9BwRZbagwa21rvr2wPloehLlobXmtwNj9KayhNwuvryyk6hcE1IsAgs25N3n8zE1NVXVwaa1ZmpqCp9Prh9YDd4ej7NlMoFLKb6wvE1CWywaC3bG3d7ezuDgIBMTEwu1y7Lw+Xy0t597t0WxOBxKZnhhcByAB7pbaAvKaEixeCxYcLvd7g+NVBRiscoaBs/0jaCBG1vquFjmHRGLzOKeHVCIBaa15tm+UVIlg66Qnxtb6stdkhAnkOAW4jibJxN8MJPG53TwYE8rjir/Ml1UJgluIWaN5wq8MGC3a3+ms1m+jBSLlgS3EIBh2SMjDW1fuUbatcViJsEtBPD7oQnGcgXqvB7u7Gg6/QZClJEEt1jy9iXSvD0ex6kUD/W04HHKx0IsbvIOFUtaqljil/0jANzSWk+r9NcWFUCCWyxZWmt+2T9K1jBZFg7wiabFP/OkECDBLZawt8bjHExmCLicPNDdUvXz6IjqIcEtlqSRbJ4Xh+zpF+7raiYsXf9EBZHgFkvOkYsimFqzvqGWC2vC5S5JiLMiwS2WnBcGx5nMF2n0e7mtvaHc5Qhx1iS4xZKyO548OlXrQ90tuB3yERCVR961YsmYKZb4df8YALe3N9IUkHnTRWWS4BZLgqU1T/cOkzNNVkZDXNVQU+6ShDhnEtxiSXh1dIr+dI6Q28V9Xc3S9U9UNAluUfUOp7O8PDyJAh7sbiHoXpzXyBbiTElwi6qWM0ye6rWvZnNdc4xlkWC5SxLivElwi6qlteY3h0eZKZZoC/q4uUW6/onqIMEtqtbWqRl2xVN4nQ4e6m7F6ZB2bVEdJLhFVZrIFXh+wO76d3dHEzGfp8wVCTF3JLhF1TEsiyd7hylZmktjEdbWRctdkhBzSoJbVJ0jV7OJed3c1SlXsxHVR4JbVJUPEineHo/jUIqHe1rxOp3lLkmIOSfBLapGqlji2f5RAD7V1iBXsxFVS4JbVAVLa57uGyFrmCyPBLm2sbbcJQkxbyS4RVV4Y2ya3lSWoFzNRiwBEtyi4g2mc2wcngTg/u4WQjKkXVQ5CW5R0fKGyVN9w1hac21jjBXRULlLEmLeSXCLiqW15rmBMeKFEi0BH7e21Ze7JCEWhAS3qFjbp5O8P53E43DwcE8rLrmajVgi5J0uKtJUvshvD9tD2u/qbKJOhrSLJUSCW1Qc09I82TtM0bK4JBbh0lik3CUJsaAkuEXF+cPwBCPZPDUeN3d3NEnXP7HkSHCLinJgJs0bY9NHh7T7XDKkXSw9EtyiYqRLBs/0jQBwc0s97SEZ0i6WJgluURG01vyyb4SMYdITDnBdc6zcJQlRNhLcoiK8OR7nQDJDwOXkwe4WHNKuLZYwCW6x6A1ncvxhaAKA+7qaCXvcZa5IiPKS4BaLWsE0ebJ3GFNrrm6s5cKacLlLEqLsJLjFovb8wDjThRJNfi+3tclV2oUACW6xiO2YmmHb1Axuh5Ih7UIcRz4JYlGKF4o8N3uV9jvam2jwe8tckRCLhwS3WHSODGkvmBara8JcUS9XaRfieBLcYtF5aWSCoUyeqMfNZ7qaZUi7EB8hwS0WlUPJDK+PTqOAh3pa8MuQdiFOIMEtFo1MyeDpvhE0cFNrPZ2hQLlLEmJRkuAWi4LWmmf7R0mXDLpCfq5vrit3SUIsWhLcYlF4ZyLOvpk0fqeTB3taZUi7EKcgwS3KbjSbZ8OgPaT93q5mojKkXYhTkuAWZVU0raND2tfV17CqVoa0C3E6EtyirH43OMZkvkiDz8OnOxrLXY4QFUGCW5TNrukk703O4Jq9mo1bhrQLcUbkkyLKIlEo8evDowB8uqORpoCvzBUJUTkkuMWCs7Tmqd5h8qbFRTUh1tXXlLskISqKBLdYcJtGJhnI5Ii4XdwrQ9qFOGsS3GJB9aWyvDIyhQIe7Gkl4HKVuyQhKo4Et1gwWcPg6d5hNHB9cx3dYRnSLsS5kOAWC0Jrza/6R0mWDDqCfm5qrS93SUJULAlusSA2TybYm0jjczp4SIa0C3FeJLjFvBvPFXhhYByAezqbqfExwFYcAAATWUlEQVTKkHYhzocEt5hXJcse0m5ozRX1UdbEIuUuSYiKJ8Et5tWGwXHGcwXqfR7uaG8qdzlCVAUJbjFv9sRTvDuRwKkUD/W04nHK202IuSCfJDEvZoolftVvD2m/ra2BFhnSLsSckeAWc87Smqd7h8mZJiujIa5urC13SUJUFQluMedeHZ2iP50j5HZxnwxpF2LOSXCLOXU4neXl4Ul7SHt3C0G3DGkXYq5JcIs5kzdMnuq1r9J+XXOMZZFguUsSoipJcIs5obXm14dHmSmWaAv6uLmlodwlCVG1JLjFnNg6NcOueAqv08FD3a04HdKuLcR8keAW520iV+D5gTEA7u5oIubzlLkiIaqbBLc4L4Zl8VTvMCVLc2kswtq6aLlLEqLqSXCL8/L7oQlGcwViXjd3dcqQdiEWggS3OGf7EmneHo/jmL1Ku9fpLHdJQiwJEtzinKSKJX7ZPwLAp9oaaA36y1yREEuHBLc4a5bWPN03QtYwWR4Jcq0MaRdiQUlwi7P2xtg0vaksQZeTB7pbZEi7EAtMgluclcF0jo3DkwDc391CSIa0C7HgJLjFGcsbJk/1DWNpzbWNMVZEQ+UuSYglSYJbnBGtNc8NjBEvlGgJ+Li1Ta7SLkS5SHCLM7J9Osn700k8DgcP97TicshbR4hykU+fOK2pfJHfHraHtN/V2USdDGkXoqwkuMUpmZbmyd5hipbFJbEIl8pV2oUoOwlucUp/GJ5gJJunxuPm7o4m6fonxCIgwS0+1oGZNG+MTR8d0u5zyZB2IRYDCW5xUumSwTN99pD2m1vqaQ/JkHYhFgsJbnECrTW/7BshY5j0hANc1xwrd0lCiONIcIsTvDke50AyQ8Dl5MHuFhzSri3EoiLBLT5kOJPjD0MTANzX1UzY4y5zRUKIj5LgFkcVTJMne4cxtebqxlourAmXuyQhxElIcIujnh8YZ7pQosnv5bY2uUq7EIuVBLcA4P3pJNumZnA7lAxpF2KRk0+nIF4o8pvDowDc0d5Eg99b5oqEEKciwb3EHRnSXjAtVteEuaJertIuxGInwb3EvTQywVAmT9Tj5jNdzTKkXYgKIMG9hB1KZnh9dBoFPNTTgl+GtAtRESS4l6jM7JB2DdzUWk9nKFDukoQQZ0iCewnSWvNs/yipkkFXyM/1zXXlLkkIcRYkuJegdybi7JtJ43c6ebCnVYa0C1FhJLiXmNFsng2D9pD2e7uaicqQdiEqjgT3ElI0raND2tfV17CqVoa0C1GJJLiXkN8NjjGZL9Lo9/LpjsZylyOEOEcS3EvErukk703O4FKKh7pbcMuQdiEqlnx6l4BEocSvZ4e0f7qjkaaAr8wVCSHOhwR3lbO05qneYfKmxUU1IdbV15S7JCHEeZLgrnKbRiYZyOSIuF3cK0PahagKEtxVrC+V5ZWRKRTwYE8rAZer3CUJIeaABHeVyhoGT/cOo4Hrm+voDsuQdiGqhQR3FdJa86v+UZIlg46gn5ta68tdkhBiDklwV6Etkwn2JtL4nA4ekiHtQlQdCe4qM54r8LuBcQDu6WymxitD2oWoNhLcVaRk2UPaDa25oj7Kmlik3CUJIeaBBHcV2TA4zniuQL3Pwx3tTeUuRwgxTyS4q8TeRIp3JxI4lX2Vdo9TXlohqpV8uqvATLHEs332kPbb2hpoliHtQlQ1Ce4KZ2nN073D5EyTldEQVzfWlrskIcQ8k+CucK+OTtGfzhFyu7hPhrQLsSRIcFeww+ksm44Mae9uIeiWIe1CLAUS3BUqb5g81TuCpTXXNcdYFgmWuyQhxAKR4K5AWmt+fXiUmWKJtqCPm1sayl2SEGIBSXBXoK1TM+yKp/A6HTzU3YrTIe3aQiwlEtwVZiJX4PmBMQDu7mgi5vOUuSIhxEKT4K4ghmXxVO8wJUtzaSzC2rpouUsSQpSBBHcF+f3QBKO5AjGvm7s6ZUi7EEuVBHeF2JdI8/Z4/OiQdq/TWe6ShBBlIsFdAVLFEr/sHwHg1rYGWoP+MlckhCgnCe5FTmvNM30jZA2T5ZEg18qQdiGWPAnuRe71sWkOpbIEXU4e6G6RIe1CCAnuxWwwnWPj8CQA93e3EJIh7UIIJLgXrbxh8lTfMJbWXNsYY0U0VO6ShBCLhAT3IqS15rmBMeKFEi0BH7e2yVXahRDHSHAvQjumk7w/ncTjcPBwTysuh7xMQohjJBEWmal8kecO20Pa7+psok6GtAshPkKCexFJFUv84tAQRcvikliES+Uq7UKIk5BuCovAeK7Am2PT7JhOYmpNrdfN3R1N0vVPCHFSEtxlYmnNvpk0704kOJjMAKCA1TVhPtXWgM8lQ9qFECcnwb3AxrJ5dkwn2TGdJFUyAHA7FJfX1XBNY61M0yqEOC0J7gWQKpZ4P55ix9QMo7nC0fvrvB6ubKjh8roofjnDFkKcIQnueaC1ZqpQZG8izd5EmsFM7ugyv9PJxbVh1tZF6Aj6pR1bCHHWJLjnSNYw6E1mOZjKcDCZZaZYOrrMpRQroiHWxiKsiAalX7YQ4rxIcJ+jVLHEQCbH4XSO/nSO0WwefdzygMvJBZEgq2rCLI8E8TglrIUQc0OC+zS01qRLBqO5AqPZAsPZPEOZHMnZLxaPcCpFZ8jP8kiQ5ZEgzX6vNIMIIeaFBPcsS2sShRLThSJThSKT+SLjuQIT+SJZwzxhfa/TQVvAR2coQGfIT0fIj1uaQIQQC2DRBvd0boodY9tJ5OPU+GrpiHQykDx89PbapkuJ+etO/0DjM9AYxbAsUsPTJKMBkiWDmWKJxIERPgjleXPvWwxk4hh+N63hdkKeEKWxKQaMUcZ0gkiiiNtfoMnr4oJQPTf7l9HT1cW+l17iHccoz7QrrnCtINbVA1oxU0gwlZ2iJqVIhDV1gTp6Sg2sWn0NMX/d0WPLDo0Q9dWwuuFiaju7j9UL0CgXAhZCnJzSWp9+JaXuAP5fwAn8f1rrb51q/XXr1unNmzefU0GGZfDM3ifZ1LcRFDiVk4PxA4ymh2kOtrIsthxLW6Dhxu5beOCih3E5Tv7/z1t9o7y39QCpthg504KhKWirA48LI5nm4CsbeTvUTySlAMVgdIqiI05NpsAduwNM+Mc4VD/EQ9vqGKzJUXCbeItOVsajTAYLXH24nsGaLD+/ZoS2aQ/vd2ew/B6UUvgLiosOeflgWYGQJ8L1Y22kL2vGE42QKWXw5RUX7tI48gaG20HTPbdw78r7cb36AaDgxtUQ8p3TcyiEqDxKqS1a63Vnsu5pz7iVUk7gvwO3AYPAu0qpX2mtd59fmSf3zN4n2di7gY5IF06Hk10TO5nJJ2gMNDFTSJDIJ7i4YQ2mZbKxdwMAn139+ZM+Vj6eZjyVg1QOh1KEcwYRCyK1YV564ycs37UPGibpSilcxRQbVg7SF05y0VCYT+7v5nBtgWDBQ/d0iJqsm5mAQTTrojXlp3PMQ0/cR6joYu1IDm0Y+GaSjDqTeJ0eapMeWouNjGTiRAwXzniOsYED7JuaYHXjWtboi+ieSeFJGhTCDl55/0U8iSL3HG4DDUylJLiFECd1Jk0l64EDWutDAEqpnwP3AXMe3NO5KTb1bTwa2rlSlv5EL1FvFKUcRLxR+hO9LKtZht8doCPSxab+jdzac5vdbGKYMJ22H8wwufxwnIuCQcJpi0CphEoU4L0BhiP7yf7qTSiVWDviJukr4Uu76RjrYlfNJNcerqc266Eh5eGavhihvItwPkLJYeEvuVGAz3CiFfgKLr74ahMHGjOsPRji3e4phusNLpwMY/ot1oxGMSwDh9tB+5YE7R1tOPeM06Ua8KcchCaKmAkHt2ZCJF99mUzsBoI1tbC1FzJ5+1hqgnbTiQzSEUJwZsHdBgwcd3sQuHo+itkxtt1uHnHYATWWsac3Vcr+0s8x+3MsM0Z3Tc/R9XaMbeem7lvAtODgKAzHwesieniSaL4ElgX5IuwcgKkUxdwQ60YCpD0GRTc04KB92ou/6OQmVy0+wwVoHJbCZzjxGk5cFijU0ekULTSWA6y8pinlZcVYmJLL4KKJEHubM7iVk1RE43S6KWkDLwadMyHaDnmpn3Hj88cptEcIpCxCfQU6MxY5h0Gm6wOCa9fYtW45BH4PXHcR1IUluIUQwJkF98n6tJ3QMK6UehR4FKCzs/Ociknk4x9qr84buaNhfWw/DvJG/uhtl3Ixk0/YN7xuuOoC2DcCfeNwZY8d4oNT4HXBg+thcJr4hn7e70kQLLjwmW6cRc3LKyYJ5Jx0zvgYCcVxaCddU0FqCi7QUMIiWvLgKoFHOzHQmC5NyWmRd5lMRQxeXD5C2m9ScmkOtuRoz0boTLswTM1QLMtQV4p1Uy6ak15i7iANXicTMRdFv4PAdImUyyKzJkrjJZ2QyEBDBD6xEla1S2gLIY46k+AeBDqOu90ODH90Ja31D4Efgv3l5LkUU+OrxbCO9Y/2ufz2F5Ef2o+Fz3Ws7dfQBlFfzbEVXE5Y3Q7x2SaTjjqYTIJSsKwJuhowt0c4ZBVpmdJ0JFwUXZr9LVkwTeKTfjQWJaemNufGgSYeLKK0Qs0oMpESdRkfPlORd1rkPSY5j8WWrgRvrUigtUahcCs347UlmvIWpsNiMmZhFRS9nSb9Ksd1CSexJDhLkKl1kqtxkiyU6GishZJhn2mv6YBLus7lqRRCVLEzCe53gRVKqR5gCPg88MX5KGZt06U8ufvnmJaJ0+GkKdjE7omdaG2hlONoiDcFmwAwLfPodh+SK0IyZ4e1acFkGtAwlYZskZVGM+1JP7EZRSzpxnJoWhJeXJaDthkPhkNjKYuGjI9w3oXXdIBW1OY9BItOwgUXTsuB16nxGRY+w6Ij7qcx6cF0ahw4MDyKkOWjJu3CtKC9GGIkX6ARLxQsWrJBQtMGWoHD1DiKFkprOnIhmMrY7fUTKftY/DJjoBDimNMGt9baUEp9DXgBuzvgj7TWu+ajmJi/jhu7bznaq8TvDtBV00Nv/CBhT5hUMUVP7XL87gCmZTKQ7OeWnttP7M89nQbDgovbweWAkbjduLOmA7b2EVEBIiuWs3l8K93jJcKWhz0dOVb1B4jm3GxrmWZr2wy+opOOuJ/RUJ6i08Q5FaWoLUyg6LAYj+RxOpyESi5K2iLnsmjI+XA5nLzTNoPf56Iu7sTv8vFBaxYfMaydGVY4O9BBN0NNbuKtbhr3ZgjvnKahsY3g7evtLyI37YFkxv7LwR+bj6dbCFGhzqgf99mak37c/RsB+wvJQ/GDjKZHaA62sKx2ORb2mfeNXR/Tj3siaQd2bejY7/aD293s8kXyLovvvvsdfjnxAu6UgdtQ+HKa+riL/c0ZMC38JQehgpuajIuM18RjOmifDOBwOzjYlMGDBw1YDotwysGh9gJulweU4nBDnkDCxO3xEPREaPE3EW3voGvCC/EsuYiTYtRNKmhRP2xyI5dw/dq7cHU3223b8TRMpuweJQ1yCTMhqt3Z9ONedMF9xJHRhTP5BFFfzdGRk0dun/HIydM4OH2AX+x6nD2Tu/G5vKxrvZrmUAsjqSHeGnqdvkQfmVKGoDtIc7iZ5dGVrG+/htZwO68d3sRAsp+SWeLihjW0RFpBK5KFGSazk/hdPnJGnvpAPW2RjqM1f/TY5upYhBCVqyqCWwghlpKzCW6ZFUkIISqMBLcQQlQYCW4hhKgwEtxCCFFhJLiFEKLCSHALIUSFkeAWQogKMy/9uJVSE0D/nD/w6dUDk2XYb7kt1eOGpXvsctzVp0tr3XAmK85LcJeLUmrzmXZgryZL9bhh6R67HPfSJk0lQghRYSS4hRCiwlRbcP+w3AWUyVI9bli6xy7HvYRVVRu3EEIsBdV2xi2EEFVPglsIISpM1QW3UuoxpdRepdQOpdQzSqma029V+ZRSn1VK7VJKWUqpqu8upZS6Qyn1gVLqgFLqr8tdz0JRSv1IKTWulNpZ7loWilKqQyn1klJqz+x7/M/KXVO5VV1wA78H1mit1wL7gG+UuZ6FshN4EHil3IXMN6WUE/jvwJ3AauALSqnV5a1qwfwYuKPcRSwwA/g/tdargGuA/30Jvd4nVXXBrbXeoLU2Zm++BbSXs56ForXeo7X+oNx1LJD1wAGt9SGtdRH4OXBfmWtaEFrrV4DpctexkLTWI1rr92Z/TwF7gLbyVlVeVRfcH/EV4PlyFyHmXBswcNztQZb4B3mpUEp1A5cDb5e3kvJynX6VxUcp9SLQfJJFf6u1fnZ2nb/F/hPrpwtZ23w6k+NeItRJ7pN+rVVOKRUCngL+D611stz1lFNFBrfW+lOnWq6U+jJwD3CrrqKO6qc77iVkEOg47nY7MFymWsQCUEq5sUP7p1rrp8tdT7lVXVOJUuoO4K+Ae7XW2XLXI+bFu8AKpVSPUsoDfB74VZlrEvNEKaWA/wHs0Vp/t9z1LAZVF9zA94Ew8Hul1Dal1A/KXdBCUEo9oJQaBK4FnlNKvVDumubL7JfPXwNewP6i6gmt9a7yVrUwlFKPA28CFyqlBpVS/0u5a1oA1wH/Bbhl9jO9TSl1V7mLKicZ8i6EEBWmGs+4hRCiqklwCyFEhZHgFkKICiPBLYQQFUaCWwghKowEtxBCVBgJbiGEqDD/P3yzq2BzJ5zWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#sklearn logit always score 100% accuracy\n",
    "logit(y,x)"
   ]
  },
  {
   "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.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}