Class1_GLM_Poisson_and_Gradient_Descent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson Regression with boston data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm # for GLM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',\n",
       "       'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "boston = datasets.load_boston()\n",
    "X_train = boston['data']\n",
    "y_train = boston['target'] #number of crime\n",
    "boston.feature_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.320e-03, 1.800e+01, 2.310e+00, 0.000e+00, 5.380e-01, 6.575e+00,\n",
       "       6.520e+01, 4.090e+00, 1.000e+00, 2.960e+02, 1.530e+01, 3.969e+02,\n",
       "       4.980e+00])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train[0, ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "y = np.round_(y_train, decimals=0, out=None) #make into integers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(506, 13) (506,)\n"
     ]
    }
   ],
   "source": [
    "print(X_train.shape, y.shape) # dimension\n",
    "# of obs = 506, # of exp v = 13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[6.3200e-03, 1.8000e+01, 2.3100e+00, 0.0000e+00, 5.3800e-01,\n",
       "        6.5750e+00, 6.5200e+01, 4.0900e+00, 1.0000e+00, 2.9600e+02,\n",
       "        1.5300e+01, 3.9690e+02, 4.9800e+00],\n",
       "       [2.7310e-02, 0.0000e+00, 7.0700e+00, 0.0000e+00, 4.6900e-01,\n",
       "        6.4210e+00, 7.8900e+01, 4.9671e+00, 2.0000e+00, 2.4200e+02,\n",
       "        1.7800e+01, 3.9690e+02, 9.1400e+00]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train[0:2,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* X_train does not contain the column of 1's. In the following, we add the column of 1's into the model matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(506, 1) (506, 13)\n"
     ]
    }
   ],
   "source": [
    "my_ones = np.ones([506,1]) # 1로만 이루어진 열 추가\n",
    "print(my_ones.shape, X_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(506, 14)\n"
     ]
    }
   ],
   "source": [
    "my_ones = np.ones([506,1])\n",
    "X = np.concatenate((my_ones, X_train), axis=1) # axis=1 : columns\n",
    "# combine two matrices so that axis=1 is increased\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.00632"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "6.3200e-03 # 6.3200* 10**(-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.0000e+00, 6.3200e-03, 1.8000e+01, 2.3100e+00, 0.0000e+00,\n",
       "        5.3800e-01, 6.5750e+00, 6.5200e+01, 4.0900e+00, 1.0000e+00,\n",
       "        2.9600e+02, 1.5300e+01, 3.9690e+02, 4.9800e+00],\n",
       "       [1.0000e+00, 2.7310e-02, 0.0000e+00, 7.0700e+00, 0.0000e+00,\n",
       "        4.6900e-01, 6.4210e+00, 7.8900e+01, 4.9671e+00, 2.0000e+00,\n",
       "        2.4200e+02, 1.7800e+01, 3.9690e+02, 9.1400e+00],\n",
       "       [1.0000e+00, 2.7290e-02, 0.0000e+00, 7.0700e+00, 0.0000e+00,\n",
       "        4.6900e-01, 7.1850e+00, 6.1100e+01, 4.9671e+00, 2.0000e+00,\n",
       "        2.4200e+02, 1.7800e+01, 3.9283e+02, 4.0300e+00]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[0:3, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(506,)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# poisson regression model -> distribution\n",
    "# X=X1+ X2+ ... +X13 (설명 변수)\n",
    "poisson_model = sm.GLM(y, X, family=sm.families.Poisson(), link='log') \n",
    "# Find W based on MLE\n",
    "# y~pois(lambda), log(lambda) = b0+b1*x1+ ... + b13*x13\n",
    "# automatically link function implied link='log' -> 따로 안써도됨\n",
    "# family : poisson, logistic, gamma regression ...\n",
    "\n",
    "poisson_results = poisson_model.fit(); # optimizing function \n",
    "# -> find beta_0, ... beta_13, maximize likelihood function of this obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                  506\n",
      "Model:                            GLM   Df Residuals:                      492\n",
      "Model Family:                 Poisson   Df Model:                           13\n",
      "Link Function:                    log   Scale:                          1.0000\n",
      "Method:                          IRLS   Log-Likelihood:                -1416.0\n",
      "Date:                Wed, 16 Dec 2020   Deviance:                       362.69\n",
      "Time:                        12:58:48   Pearson chi2:                     389.\n",
      "No. Iterations:                     5                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          3.7780      0.232     16.268      0.000       3.323       4.233\n",
      "x1            -0.0103      0.002     -4.699      0.000      -0.015      -0.006\n",
      "x2             0.0010      0.001      1.900      0.057   -3.29e-05       0.002\n",
      "x3             0.0022      0.003      0.829      0.407      -0.003       0.008\n",
      "x4             0.0922      0.035      2.664      0.008       0.024       0.160\n",
      "x5            -0.6932      0.182     -3.812      0.000      -1.050      -0.337\n",
      "x6             0.1225      0.018      6.735      0.000       0.087       0.158\n",
      "x7             0.0004      0.001      0.627      0.531      -0.001       0.002\n",
      "x8            -0.0490      0.009     -5.614      0.000      -0.066      -0.032\n",
      "x9             0.0154      0.003      4.827      0.000       0.009       0.022\n",
      "x10           -0.0005      0.000     -3.114      0.002      -0.001      -0.000\n",
      "x11           -0.0359      0.006     -6.281      0.000      -0.047      -0.025\n",
      "x12            0.0005      0.000      3.195      0.001       0.000       0.001\n",
      "x13           -0.0314      0.003    -12.235      0.000      -0.036      -0.026\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "print(poisson_results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 여기 있는 coef가 각각 b0 ... b13\n",
    "# 목표는 lambda를 찾는 것 ... mean of y\n",
    "# lambda hat : beta hat을 사용했으므로\n",
    "# lambda : always positive -> log 함수 안에 있음이 이를 보여줌\n",
    "# -> exponential (X(i)*beta hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = poisson_results.mu # prediction for each observation\n",
    "# yhat[0] ... 29.9xx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [],
   "source": [
    "# yhat = poisson_results.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's calculate yhat[0] using hands\n",
    "# X[0,]\n",
    "# yhat0=exp(3.7780+beta1*x1+beta2*x2+... +beta13*x13)\n",
    "# ==  yhat0=np.exp(np.dot(X[0, ],beta_hat)) ...29.9xx same"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.7256638587818"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((y_train-yhat)**2)  # Calculation of MSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Calculate MSE again!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = poisson_results.params # Estimated parameters\n",
    "# ... (14,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_hat = poisson_results.params\n",
    "beta_hat.shape \n",
    "beta_hat = np.reshape(poisson_results.params, [14, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(506,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "16.7256638587818"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_yhat = np.exp(np.dot(X, W) )\n",
    "print(my_yhat.shape) \n",
    "# yhat = np.exp(np.dot(X, beta_hat)) # real y.hat\n",
    "# print(yhat.shape) # (506, 1)\n",
    "\n",
    "# y_train : (506,)\n",
    "np.mean((y_train-my_yhat)**2)\n",
    "# [506,]-[506, 1] -> [1, 506] - [506, 1] -> [506, 506] - [506, 1] ... ugly\n",
    "# Y = np.reshape(y_train, [506,1])\n",
    "# Y.shape -> (506, 1) ... 일치\n",
    "# MSE = np.mean((Y-yhat)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Draw plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2641d43eec8>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0EUlEQVR4nO3deXhTVfrA8e/bEKBFoaCoWMEqKggiFHEb1BFUcLfjOu7jgIjrCPzQMqCgg4IwLrjhoOLKKLhVBFkUcFewWHZxFERqQChLQaFAac/vj5uUNM1NbtKkSdr38zw+pbfJvecUeXv6nnPeI8YYlFJKpZ60RDdAKaVUdDSAK6VUitIArpRSKUoDuFJKpSgN4EoplaIa1ObDDjzwQJOdnV2bj1RKqZS3cOHCTcaYloHXazWAZ2dnU1BQUJuPVEqplCcivwS7rikUpZRKURrAlVIqRWkAV0qpFKUBXCmlUpQGcKWUSlGOAriIZIrI2yKyUkS+F5FTRaSFiHwkIj96PzaPd2OVUkrt43QEPg6YaYxpD3QGvgfygDnGmKOBOd7PlVJK1ZKwAVxEmgJnAC8CGGP2GGNKgEuAV7wvewXIjU8TlVIqhW3eDHffDdu2xfzWTkbgRwLFwEsiUigiL4hIE+BgY8x6AO/Hg4K9WUT6iUiBiBQUFxfHrOFKKZXUjIG33oIOHeCZZ+Czz2L+CCcBvAHQFRhvjMkBdhBBusQYM8EY080Y061ly2o7QZVSqu5Zvx4uvRSuvBJat4aFC+Gii2L+GCcB/FfgV2PMfO/nb2MF9A0i0grA+3FjzFunlFKpxBiYOBGOPRZmzoQxY+Cbb+D44+PyuLAB3BjzG1AkIu28l84CVgBTgRu9124E3o9LC5VSKhWsXg3nnAN9+kDnzrBkCQweDA3iV3LK6Z3vBCaJSENgNXATVvCfIiJ9gLXAFfFpolJKJbHycnjqKRg6FFwuGD8e+vWDtPhvs3EUwI0xi4BuQb50Vkxbo5RSqWTFCmvE/c03cP758NxzVs67luhOTKWUitSePfCvf0FODvz4I7z+OkybVqvBG2q5HrhSSqW8ggJr1L1kCfz1rzBuHBwUdBV13OkIXCmlnNi5E+65B04+GTZtgvffhzfeSFjwBh2BK6VUeJ9+Cn37wk8/wc03w9ix0KxZolulI3CllLK1fTvceiuceSZUVMCcOTBhQlIEb9AArpRSwU2fDh07WgF74EBYuhR69kx0q6rQAK6UUv42bYLrroMLL7RG2l99BY8+ChkZiW5ZNRrAlVIKrG3wb75pbYOfMgVGjIDvvrMmLZOUTmIqpZTHY+W6P/gATjoJXnwRjjsu0a0KS0fgSqn6yxh4/nmr5OvHH1upkq++SongDToCV0rVV6tWWUsC582DHj2sQN62baJbFREdgSul6pfycnjsMejUyarTPWGCtTwwxYI36AhcKVWfLFtmbYNfsMA6YGH8eMjKSnSroqYjcKVU3bdnj7WqpGtX+Plna7XJ+++ndPAGHYErpeq6BQvg73+H5cvh2mvhiSfgwAMT3aqY0BG4Uqpu2rkTBg2CU0+1ToSfNs0q+1pHgjfoCFwpVRfNm2cVn1q9Gvr3h0cegaZNE92qmNMRuFKq7ti2zTrOrGdP60izTz6xJirrYPAGDeBKqbpi6lRrQ86LL1qHCS9eDH/+c6JbFVcawJVSqW3jRutknEsugQMOgPnzYcyYpCw+FWsawJVSqckYmDTJGnW/9551RmVBAXQLdv563aSTmEqp1FNUZBWfmj4dTjnFSpt06JDoVtU6HYErpVJHRQU895x10MK8edaa7i++qJfBG3QErpRKFT/+aBWf+vRTOPtsq4bJEUckulUJ5SiAi8ga4HegHNhrjOkmIi2AyUA2sAa40hizNT7NVPVFfqGHsbN+YF1JKYdmpjO4dzuAatdyc7KCvjY3J/m2RkfTzvxCDyOmLqektAyA5hluhl/U0fZ9gc/o0b4l81YWh3ym3ff6gQ+Ws3Wn9dzMdDcjLrZ/bqj7+Lc/TaDCQFaQtoT9/uzdC48/DvffD40aWemSm24CkZDfw/pAjDHhX2QF8G7GmE1+18YAW4wxo0UkD2hujLk31H26detmCgoKathkVVflF3oY8u5SSsvKK6+50wQEysr3/X+a7nZx2QlZvLPQU+W16W4Xoy7tlFRBPFifwrUzv9DD4LcWU1ZR9d+m2yWMvbxz0EAc+IxAgc8M+r12CeUVhoDH4k4Txl5R/bm290kTyk31+wRrS9jvz+LFVvGphQshNxeeeQYOPdS2n3WViCw0xlSbna1JDvwS4BXvn18BcmtwL6UYO+uHakGorMJUCd4ApWXlvDG/qNprS8vKGTvrh7i3MxLB+hSunWNn/VAteIP1QyzY+4I9I1DgM4N+r8uDB92yiuDPtb1PkB8Cdm2x+/48MX0p3HeftaKkqMg64uzdd+tl8A7FaQ7cALNFxAD/McZMAA42xqwHMMasF5GDgr1RRPoB/QDatGkTgyarumpdSanj15bb/OYYyT1qg117QrUz0q857bP/6yL9PkXTDyf3C/b+rp7veWTGk7C5CG64wardfcABUT2nrnM6Au9ujOkKnAfcLiJnOH2AMWaCMaabMaZby5Yto2qkqh8OzUx3/FqXTf4zknvUBrv2hGpnpF9z2mf/10X6fYqmH07u5//+9D27uP/jCbz9+j3sv3c3fPghvPKKBu8QHAVwY8w678eNwHvAScAGEWkF4P24MV6NVPXD4N7tSHe7qlxzpwluV9Vgne52cfXJrau9Nt3tqpxASxbB+hSunYN7t7Ny/wHcLgn6vmDPCBT4zKDfa5f9pGCP9i3JL/TQffRcjsibTvfRc61cvc3fWZDmB22L7/3d1yxi9sTb+fvCqbzR7UIKpn0G550Xsk/KQQAXkSYisr/vz0AvYBkwFbjR+7Ibgffj1UhVP+TmZDHq0k5kZaYjWCsWxl7RmbGXd65ybdSlnRiZ26naa5NtAhOC9ylcO3Nzshh7RWcy092V15pnuINOYPo/w//1TRq6yEx32z4z6Pf68s40z3BXuz/A9CXrGfLuUjwlpRjAU1LKkHeXAgT9O3vsyi5V2uML6IFtyc3OYObSl5k0eRhlrgbc1u9xmkwYz4WnJdcP4mQVdhWKiByJNeoGK2f+X2PMQyJyADAFaAOsBa4wxmwJdS9dhaJUfESz2iWYI/KmE35d2j5Zmel8mdczgnf4ee89uO02KC62ik8NHw6NG0d3rzrObhVK2ElMY8xqoHOQ65uBs2LTPKVUTYRa7RJJAD80Mx1PBBOTUU1ibtgAd94Jb70FXbpY2+G7do38Pkq30itVF8RqlYhdzt4/HeIvoklMY+DVV+HYY63zKB96yDruTIN31HQrvVJ1gN3IOdJVIr7RerCdlcFSNI4njdeuhVtugZkz4U9/snZTtm8fUdtUdRrAlUoy0Wy9H9y7Xc0CrJ/cnCzb50VcuqCiwjoRJy/PGoE/9ZSV907TX/5jQQO4UkkkcDLSf7VHuJUrEEWAjUCowB7UDz9Y51J+8QX06gX/+Q9kZ8esPUoDuFJJpSaTkREH2Cg4+u2grAwefRRGjLBOxXn5ZWtHpRafijkN4EolkVhvWY8lR78dFBZaxacKC+Gyy+Dpp+GQQxLV5DpPE1FKJZFYb1kPJdjOylBCFubatQuGDoUTT4R16+Dtt63/NHjHlQZwpZJINFvvo+EbTQfurAwVxO1+C2i1bKG1nvvhh+H662HFCmv0reJOUyhKJZF4Tkb656/TRKpVdAyXaw9cqthk904Gf/YqNxROhzZtYNYsa7JS1RoN4EolmXhMRgbmr6Mpx+u/VPGM1Qt5eNbTHLp9Ez9ffRNtJ4yD/faLaZtVeBrAlaoHnBz6AKFz7bk5Wbi3bUUGDeL872azpmUbvpj4LmfclBvDlqpIaABXKoVEew6ok1UsYXPt77zDBbffDps2wdChZA8bRrYWn0ooDeBKJVAkATnaTT5gv9XeJUKFMaGfvX493HGHdaRZTo61Hb5Ll8g6quJCV6EolSCRrgSJ5nxNH7vVLY9e2ZnHr+oCwIDJi6ouJzTG2oTToYNVMXD0aKv4lAbvpKEjcKUSxMmuS/8Rul2dbifpEadFqnw/RDLWFdFr3P3w0Udw+unwwgtwzDFR9lTFiwZwpYg+t1wT4XZdBjukIRinm3yCrW7pPnpulfunVZRz1ddTOW3Mq9CwATzzDPTvr8WnkpQGcFXv1SS3XBPhSsA6WTlS000+/j9E2m4q4pGZT9LN8z2fHHECZ37yrrW+WyUt/bGq6r2a5JZrItSuy/xCT8iTcWJ1Duihmek0KN/L7V9N5sOX76Tt5l8ZcMFAht78iAbvFKAjcFXvJaqAVLi8tJ0anUMZYGSbPRzyxECO3bCaae1PZ8TZ/diReSCjztXDFlKBBnBV78XqNJtoOMlL+4tFXZT8Qg/jpi3hyukTufnbd9nZ7ADyrnuAyVkncGhmOsNqIf+vYkMDuKr3YnmaTaBoJkdDjfxrmjLJL/TwzuP/5cVpT3Dk1nW8eXwvHut1M/+85lRGa9BOORrAVb0XrwJS0U6O2v1GkJWZXrM2bd9O+W238do3U1nb7GCuuWokX2V3AYj49HqVHDSAK0V8CkhFe7pOXH4jmDEDbrmFvxT9yovdLuHfp19PacN92+CT4cAIFTkN4ErFSbSTozH9jWDzZhgwAF57DTp0oP+tTzK76RHVXlYb+X4Ve44DuIi4gALAY4y5UERaAJOBbGANcKUxZms8GqlUKqrJ5GiNfyMwBt56y6phsnUr3HcfDB3K+Ss28Xmc8v2q9kWyDvwfwPd+n+cBc4wxRwNzvJ8rpbyCrfMWoEf7lvF98Lp18Je/wFVXweGHw8KF8OCD0KgRuTlZjLq0E1mZ6TFbS64Sx9EIXEQOAy4AHgIGei9fApzp/fMrwCfAvbFtnlKpKzcni4JftjDpm7WVdUwM8M5CD90ObxH7oGkMTJwIgwbB7t0wdizcfTc0qPrPvDZOr1e1w+kI/AngHqDC79rBxpj1AN6PBwV7o4j0E5ECESkoLi6uSVuVSjnzVhZXK0IVl12eq1fD2WdD375WtcClS+H//q9a8FZ1S9gALiIXAhuNMQujeYAxZoIxppsxplvLlnH+1VGpJBPNRGZEp8WXl8MTT0CnTvDtt/DcczB3Lhx1VA1brlKBkx/P3YGLReR8oDHQVEReBzaISCtjzHoRaQVsjGdDlUpFkU5kRrR2fPly6NMH5s+HCy6wgvdhh8W2AyqphR2BG2OGGGMOM8ZkA38F5hpjrgOmAjd6X3Yj8H7cWqlUigpVsCoYR4W19uyxJiVzcuCnn2DSJPjgAw3e9VBNEmSjgSki0gdYC1wRmyYpVXf4r+n2lJTiEqkSkANH1WFTLt9+a426ly6Fq6+GceNAU5P1VkQB3BjzCdZqE4wxm4GzYt8kpeoWX5B2khqxS7kckSEweDA89hi0agVTp8JFFzl6fiIOq1C1Q6eolYqRwEDZo31L5q0sZl1JKWkilJuq61GCbasPto3+DM9ynv10PBStgX79YMwYaNbMcZsScViFqh0awJWi5qPUYIHy9W/WVn49MHj7BKZM/FMuv2/YxINfvUbugmnQtq21uqRHj4j6FW09FpUaNICrei8Wo1Qnx58FE2w1Sm5OFrmeQnh6IKxfb23MefBByMiI+P6JOqxC1Q49Uk3Ve7E4Ui3agLhzz96q67yLi+Gaa6z8dvPm8PXX8O9/RxW8wX65ohavqht0BK5SWrjUh5PUSLSjVP97I1Bty2UQItaOd5+tO8sY8u5SCtZsxjV5CndNfYqme3byU/9BHDvuYWjYMPxNQ4jnYRUq8TSAq5QVLvXhNDUSTdXAwHs7Cd7pbheNGqRRUlpW5XqzzRvoMfh+zlr1LYtaHcM9591F0YFtGbW8uMZ56ngdVqGSgwZwlbLCTdA5ncALNkoNVzUwkpy3QGXgHDB50b7rpoK/Lp7NkHkTcVeU86+efXnphIuoSHNBDCcatXhV3aUBXKWscKkPp6mRaKoGOs15B54gP2LqckpKyzh86zpGz3yKU9cu5as2x5N37p2sbd6qynuD/VaglD8N4CplhUt92H29Wbq72rVQVQODBXC7e/tzu4Qdu/eSnTcdl3cduKuinL4F7zPo80mUpbm499w7mXx8Lys5HkCwUjU6elZ2dBWKSlnh6owM7t0Od1r1wLgjcOUH9iNqT0lp0MqAwZ7tShN8T0sTKC83lfnucmNoV7yGd17/P4bNm8gX2V04p++zTO7cO2jwBuu3gJiXnVV1io7AVcoKN0GXm5PFAx8sZ+vOqpOGZeWm2sjabkQt7EtlBJsE9T07M8PNH7v24suKV/gN5xvuLeP2r6dw2zdT2NZ4P+64+B6mtT/dNnD70/XaKhQN4CqlhZugKwkI3j6BgdFuIjNUWsX/2d1Hz632gwKgy7ofeGTGONptWsu7HXvwr5592ZpRfRu83SpEXa+tQtEAruo0p0sEg43m7XLcwUbFgdfS9+xi0Oev8feCqfy2/wHcdPlw5rU9Mej9mme4ueD4Vryz0KPrtVVENICrOi3cRpZQG326j57reH14s3R3Zb771F8WM3rmUxxe8huv5ZzPI3/+G3802reT0jfazvIWvJq+ZH1l3RT/r+l6bRWOTmKqhIno6LAohTqF3bcZx1NSimFfjjvURKUvJ+7f3vxCD9t3ldF01x+MmvEkb7w5lAoRrrp6FPf3vo0/GmXg8ua7szLTufaUNjTPcFcWvPJPvRjAnSYavJUjYmyqpMVDt27dTEFBQa09TyWvajsZsUbGvuAaz+f6RtzBSrxC1bXbvtd7Skqr5anT3S4uOyGLdxZ6OG3Fl4yc/SwH7ijh+ZP+wuPdr2G3uxEC/Dz6girPH/z2YsrKQ/+7C1w/ruo3EVlojOkWeF1TKCohElHmNPCHhl2JV/+0iW+isssDs6ttgS8tK2f23KWM+eg5Llr5Od+3zKbvpfextNXRla8JTLeMnfVD2OANVk5dD2JQ4WgAVwmRiDKnTre/uwKW9+UXeqoFb4whd8UnDP94Ahllpfz79Ot47uTL2eva908qcBIyv9DjeHdls3S3HsSgwtIArhIimgJSNeX0h0PgyDxwM02r7cU8NOsZeq4u4LtD23HPef/gpwPbVHmNS4RRl3YC9k2Ghl/1bXGnCSLoQQwqLA3gKiESUebUyfZ3sPLP/nyBX0wF1y6ayb2fvITLVPDAWTfzStcLreJTftwu4aoTW1fbRORktikz3c2IiztWKXoVrC1KgQZwlSDxLnMaLH88uHc77rYJjD7Blhga4IgtHkbPeJKTf13O54d3Yci5d/Br5iFB71FWbpj8bZGjXDfsC9r+ffdNnAbSjT3Kn65CUXVOqBUuwbbW+/jWZc9bWVyZ8kirKKfvt+8x4Iv/stvlZmTPvrzV6WxH2+CdsFttkqhVOio56SoUVW888MFy2/zx8Is6VguMbpfQpGEDPCWlVUrKtt+4mjEfjqPThlXMPOZU7jvnVor3axGzdoZKGelBDMoJDeCqTskv9NiOsNeVlNoWofKtMjFYxafu+OpNbp3/NiWN9+fWS/KY0a57TEbdkey01IMYVDhhA7iINAY+Axp5X/+2MWa4iLQAJgPZwBrgSmPM1vg1VdVnTtdEhyq/6ssf+wfGLg/MpsyvdGBXz/c8MuNJjt5cxDvH9eRfPftSkt608usZ7jR2llWEbGsa4HJJtRx4ZrqbCzu3Yt7KYtaVlFa2VYO0ipaTEfhuoKcx5g8RcQNfiMgM4FJgjjFmtIjkAXnAvXFsq6ojIj2IuEf7llUKPXlKShkweREFv2xhZG6nKvcOtUojMF3hv747Y08p//fZa/xt4Qesa3ogN17xAJ8eeUK1e+wqqwh5frFvQhKqpz8AXdutYipsADfWLOcf3k/d3v8McAlwpvf6K8AnaABXYURzELF/XtrHAJO+WVvtyDO7pYKZ6e5qQdI3Aj7t50JGzXqa1ts28HLXCxl7xg3s8Cs+5c9u7B1sgjHwed1Hz9W13SqmHBWzEhGXiCwCNgIfGWPmAwcbY9YDeD8eZPPefiJSICIFxcXFMWq2SlWhttDbfd1utBvsxBq7U3p8o2Kf/EIPv/9WzJgPn+D1Kfexx+Xm8msfYcQ5/W2Dtx3fpp1wQTgRu09V3eZoEtMYUw50EZFM4D0ROc7pA4wxE4AJYC0jjKaRqu6I9iBiO4GjbV8Q9V8u2KhB1XFKfqGHOQ8/x8cznqHFzm08c8oVPNn9anY3aBjRs30qjHE0gk7E7lNVt0W0CsUYUyIinwDnAhtEpJUxZr2ItMIanSsVUrQHEduxO/h3l99EY0lpGQMmL+KtgrVsX+Oh/ztP8NQPX7L8oCO56fLhLD/kqOg6E9D2cLn9ROw+VXVb2BSKiLT0jrwRkXTgbGAlMBW40fuyG4H349RGVYeEO4i4R/uW1WqGpLtddG8bfP11sDRK0DSMMRyS/xavPXYTZ/+0gDFn3MAlNzxW4+At3jaHqy0OoWuTKxUNJyPwVsArIuLCCvhTjDHTRORrYIqI9AHWAlfEsZ0qCcSivGmoDSr5hR7eWeipkvMW4LITshiZ24nsvOlB7xmYdgn8PGvbRh6e9TR//vk7CrKO5d7z7mLVAa0jarcdA7yz0MO0xesdTVDq2m4VS05WoSwBcoJc3wycFY9GqfgKFYjtvhZu9UgkAoOY72SeYKkTA8xbaU1+ZznMIfvSMGIquP676dz76SsA3H/2LbzW9QKMxPYgqtKyctsytTpBqeJJd2LWM6ECMdivU47XAQzBan4E8gVBpznkwb3b8dSz0xg940lO9Kzg0yO68s/ed+BpFnShVBWh1nhHQycoVTxpAK9nwi3js/tavJbAOTlkwX8Hpe89tmmcsjKOnvg0H770KKXuRgw6fwDvHNfT0TZ43xFpvmJWkWie4eaP3Xur7L50u0QnKFVcaQCvZ0IFYruRp6ek1HH6wp+TnHm4HwCBI+yQOeTCQkquuYGOK5cxvV13Rpzdn+L9mgd9qX/lwUhOpLezu6yc8sDysbpoVsWZBvB6JtQyvt+27Qp6TqRLJKIlcPmFHkZMXV7lGDK7nHmoZYO+IDt21g8MmLzIfuJ01y548EEYM4a96U25JfefzGr3J9vvwXWntKm2BT/Q4N7tGDB5UdAYnJnuRoQqRbOC1UcpqzC6y1LFVWxnc1TSC7WMz+6Q33LvRhUnS+B8Oe1qZ0hSNVUTqj1gpSR8NVBCLc3jiy+gc2cYNQpuuIGefcaHDN6wb1I0lNycLK49pU3QJY0jLu5IRkNnYx+dxFTxpAG8ngkViAOPEvPJ8stBf5nXk59HX8CXeT1tqwE6mZAMbE9murvK9a07y5j0zVr7fP3vv8Mdd8Dpp8OePTBrFkycyPbG+4X9HjgNqiNzO/H4VV2Cfq+c3kMnMVU8aQqlHrLLI8dip2C4wBYsoOXmZDF21g/VRu12KeSjvvsCjrsJiorgrrvgoYdgv/CBO1Qb7Nh9r5zsGNVdlireNICrSrE4BSZUYAsV0JyMaJuV/s79c5/nsmVz+T37KPb/4gvy0w9n7NMLKtvbpKGLHXvsfwPwtaGmm5KC/bDzneyzrbRMT9BRtULPxFQxZbeuu0lDF25XWpXgBvt+WIhARZD/FQVrG/x5P3zJgx89R+au33nu5Mt5+k9XYRo1wkDVpXtpQgVQHuRmWZnpZB+QzlertxD4v30kJ+X491WPPFO1we5MTA3gyjGnAWtY/lLemF9EuTG4RDjlyOYsWLO1SqBNE3ClVT+1JlDr0q0Mnfks5/7va5Ye3JZ7z/sHKw4+MuR7MtxplJZVVEnBpLtddG3TjC9XbQnbTz08WCUbDeCqRuxOSfdtfLE7PSdqxnDF0o8ZNvcFGpWX8fhp1/DCiX+hPK36ipV4sDstXqlE0FPpVY3Y7eD0Py3HU1LK69+srfGzDiv5jVEzn+b0XxYx/7CO5J13Fz+3qN3RsC7/U6lAA3iSSPZ8ql1Ai+Xvb2kV5dzw3XTu+ewVKiSNYb1uY1KXc2NefMpRW0Q4Im96Uv5dKOWjATwJxLLSXyTPjOQHRqQHLUSq7aYixswYxwnrVjLvyBMY2vt21jUNX3wqXnybmvTgYZXMNIAngUgr/UU7Wve9z1NSWqXqnpMgFWzZXCw0KN9L//lvc+dXb7LTnc7dFw4iv8OZjopP1RY9eFglKw3gSSCSSn/RjtYD3xeY+ggXpHJzsij4ZQv/nb826HK/aBz320+M/fAJji1ew7T2pzP87FvY3CQzNjePMc2Jq2SkW+mTgN3OwGDXw5WDteOkbKunpLRqnRE/vtNyYhG8G5XtJu+Tl3j/1YG0KN1Ov78M5Y5L7k3a4A26JV4lJw3gSSDcOZH+oq3L7XQEWa1YlNcDHyyPSfrkpKJlzHjpTvrPf4cpnc7mnD7PMvuYU2t831jp3raF478LpRJNUyhJIJIt7OFOdbfjdBIyWColv9BTpXRqNPbbvZN7P32Z6ws/ZG2zg7nmqpF8ld2lRveMh0k3n5r0K4KCScU2q5rTAJ4knB52G23BqUgmIQNH6+HSM+GcuepbHp71DIf8vpkXul3Co6dfT2nDxjW6Zzz4V11MpeCXiFVMKjloAE8x0RacCnxfmoht/e/A0Xy0E3jNd27j/jnP85cVn/C/A9pw2XV5FGa1j+pe8ZbKaZJ4nVeqkp8G8BQU7QjR/31H5E23fZ2npJTuo+dW/mCIeA24MVy48nNGfPwfmu36g3F/uppnTr2SPQ3c4d9bC4KVAEjllEO8zitVyU8DeD0VLih7SkoZMHkRBb9sYXDvdtw9eZGj+x70+2Yemv0s5/w0n8WHHM11V41k5UFHxKjV0anrZV6jnRdRqU8DeD3Vo33LKnVMgjHApG/WMn/15vA3NIarlsxm6LyJNCwvY2SPv/NSt0viXnwqTaBRA6v6IFhHsV1wfKs6M7p2IhYHcajUFDaAi0hr4FXgEKACmGCMGSciLYDJQDawBrjSGLM1fk1V0fJfodAs3U1ZeUXIQw/8GeDHjTtCvqZ1yW+Mnvkk3X9Zwjetj+Pe8+7il+aHxqDl4VUYaNGkUb2uHBiLgzhUanIyAt8LDDLGfCci+wMLReQj4G/AHGPMaBHJA/KAe+PXVGUn1BKywBUKwQ4bjlZaRTk3LfyA//vsNfampTGk9x282blXzIpPCVSWqJ23stg25aO53tRbOaNiI2wAN8asB9Z7//y7iHwPZAGXAGd6X/YK8AkawGtdqCVkAIOmLLZdbVITxxSvYcyMJ+my/n/MaXsiQ3vdzm9ND4zZ/YPV4+4+eq7mepXyE1EOXESygRxgPnCwN7hjjFkvIkFLx4lIP6AfQJs2bWrUWFWd3RKyuycvqlKwKlbc5WXc9vVb3P71FH5vlMFdFw1m6rFnxLT4VBoEzd9qrlepqhwHcBHZD3gHuNsYs10c/oM1xkwAJoB1Ik80jVT2QqUPYv3NPn79/xjz4Tjab/qF/A5/5sGz+rElo1mMnwLNMtxB0wGa61WqKkcBXETcWMF7kjHmXe/lDSLSyjv6bgVsjFcjlb141+kGaFy2i4GfT6JPwftsbNKcPpfdx5yjTo7b87buLKuyDt2f5nqV2ifsbJNYQ+0Xge+NMY/5fWkqcKP3zzcC78e+eSqcYIWwYumUtUuYOfFO+n37Hm927kWvvs/GNXj7+HL5dtURlVLORuDdgeuBpSKyyHvtn8BoYIqI9AHWAlfEpYVJJtZFg2p6P99rYz1Zuf/uHQyZ9xLXLJ7JmsxWXP3Xh/n68ONjdn8ndDu4UqE5WYXyBdaKrmDOim1zklusiwbF6n65OVkMcLhT0omzfprPQ7OeoeWOEv5z0qU8fto17HInpviULhFUyp7uxIxArIsGRXM/uxF7s3R3jdd4t9i5jeEfT+CS7z/l+5bZ9Lt0GEtaHVOje/pLd7sYdWknAIa+t9TRZiJdIqiUPQ3gEYh10aBI72c3Yi/4ZQs79uyNqg0AGMPF33/KiI8nsN/unTx22rWMP+Vyylw1Kz6VJlg/WHZWr0GSm5NVbYfojj17KSvflwbSJYJKhaYBPAKxLhoU6f3sRuxvzC+KOv99yPZNjJz9DGev+pbCVu2457y7+LHl4VHdK1CFgYyGDSi8v1fQrweuKNFDCZSKjJg47NKz061bN1NQUFBrz4u1wBEw7EsLxCIHDlRuvnF563Vn+QWyI/Kmx2xtt5gKrl48iyHzJtKgooJ/n3E9L51wERVxKD6VpcFYqRoRkYXGmG6B13UEHoFYbyTxv5+npLTKzknfiNp/YtNuxO4KcThDMNlbPIye+RSnFC3jy8OPJ+/cuyjKPCSqPjihJ8QoFR86Ak8SdnU+fHyj2GDHoqW709hbYarkj4NxVZTz92/fZ9AXr7MnrQEP9ezD5ON7xXQbfCjB6psopcLTEXiSCzcR6ikprRy9PvDB8iqHDPtqYfs0apBGuttVZVVK+40/88iMJ+n82498dNTJDOt1Kxv2j13xKSd0SaBSsRWbup+qxpxMhA7LX0puThbhfmnavbeCjofuD0DDvWUM+HwSH7xyN1nbN3L7xfdy86XDaj14gy4JVCrWdASeJJycGj/pm7WAs5reX67aQo5nJY/MeJJjNq/l3Y49ePCsmylJbxqzNttJd6dV+60ArFOAlFKxowE8hiJZBud7raektHISMjPdTWN3WpX0iD8DvDG/KGw70vfsYtDnr/H3gqn8tv8B/O3y4XzS9sSadA3AUXlad5rQ2O0KGsDnrSyucRuUUvtoAI+RSLbFB77Wt4KkpLSMdLeL5hlu2yAebrXJn9YsYvTMp2izbQOv5ZzPI3/+G380yqhR33zCBe/MdDcjLu5ou61fc+BKxZYG8AiEGmGHOlhh7Kwfwr7W/z2h0ihpYm2QCdR01x/8c95E/rpkNqubH8qV14xmQevjouxpdc0z3Gwv3Rv0B0jg6hLfbxaBNAeuVGxpAHco3Ag71OgykteG4nYJ5UGWCp7z4zeMnP0sB+woYfzJl/NE96vZ7W4U1TPs/LErePAOtt1dT85RqnboKhSHQhWegvCjy0he68/lXaMtAmXlBv/M8oE7tvL0+4/w/Lsj2ZzRjNwbHuORM/8W8+AtQFmQYb9LJOgu1NycLEZd2omszHQEa4Qe7W5VpZQ9HYE7FK7wlJNVJL7X9mjfkte9K0rCKTcGt0uqbtIxhr8sn8f9c54no6yUsadfz39Ovoy9rtj/dYaauKwwxjYo68k5SsWfBnCHwhWeCtwWb3cPiGw1hkuqBu9Dt2/koVnP0GP1QhYe2p57zvsHqw5s7fh+kQo1cak5baUSS1MoDgU7uiwwr5ubk8WXeT154qouIV/rNAee7nZV5p3FVHDdd9OZ/eLtnFy0jBFn9eOKax+Ja/AO1zbNaSuVWDoCdyiSQlbhXms3ms9Md9OkUYPK+tgiVu78iC0eRs94kpN/Xc7nh3dhyLl38Gsci0+Fo9UFlUoOWswqAfILPQycvKjKhGQa8NhVXcjNyWJY/lImfbOWtIpybl7wHgO+mMSuBg0Z2bMvb3U6u9aKTwUTbUEqrfWtVPS0mFUSKfhlC4H7FCu818HaMt9+42rGfDiOThtWMfOYU7nvnFsp3q9FjZ5rt4bcqWjTJrE+S1QpZdEAngB22+HfmF/El0s9DPzsNfrPf5uS9P3pnzuEme26x+S5TRtHfm6mS4QKY2o0ao71WaJKKYsG8ASw2w7fuWgFY2aM46gtv/L2cWfxr5592Za+f8yeW1JaFtHhDzU5bchfrM8SVUpZNIDHSaicb2AQzdhTyuDPXuXGhdPY0KwlN1zxAJ8deUK1e3Zv24Ll636v0enzdrspfafFxyNPHeuzRJVSFp3EjIP8Qg8Dpyyqlm9unuFm+EUdKfhlS+VGntN//o5RM5/m0O3FzD//rwzqfCXrKoKfBu+kGqBTsUiNOBXrs0SVqm+insQUkYnAhcBGY8xx3mstgMlANrAGuNIYszWWDU5l/3x3SdDJwq07yxjy7lJGXdqJxn9sp/2YEVy+9GNWtTiMFx96iZv/eSPr86bb3tf/sOOaqjCGn0dfUOP7OBHrs0SVUhYnKZSXgaeBV/2u5QFzjDGjRSTP+/m9sW9eatoZpBa2T2lZOQsefYGH5zwHxcUwZAht77+fto0bk1/oIS1MgC43JiYj8dpOX+jWeqViL2wAN8Z8JiLZAZcvAc70/vkV4BM0gIfV8o+tPPDReM7/31fQpQt8+CHk5AD70gzhRtfNM9x0aLU/X67aEnU7dBelUnVDtJOYBxtj1gMYY9aLyEF2LxSRfkA/gDZt2kT5uNQiQtVzK43h8mVzGDb3BdLLdvNc7z70/2A8uPflukPVCPf3x669rFj/e9Rt012UStUdcV+FYoyZAEwAaxIz3s+LJ7uVJYHX/3Rki8oR8mHbNvDwzKc5Y00h32Z1IO+8O1l1QGtee/TzKoHU6ZK6sgpje1pPKDppqFTdE20A3yAirbyj71bAxlg2KhnZ7SYs+GULkxcUVdbL9pSU4ikppVGa4epvpzH401cwItx3Tn8m5ZxPhaRVeT9Y+WG7pXaRaJ7hJqNh1VoqJTvLyMxwYwwMCHI6kFIqdTlaRujNgU/zW4UyFtjsN4nZwhhzT7j7pPIywu6j5wYNsMEmFNtuLmL0jKc40bOC/+X8iWPe+y/d31gV9P2+5XyZGW7+2LW3ysEJkUxW2o2wdQmfUqnPbhlh2HKyIvIG8DXQTkR+FZE+wGjgHBH5ETjH+3mdZpfi8A+wDcr3ctvXU/jwpTs5evNaBl4wgN7nDCF/SwPb95cbg8FaYohYFQl9p9hce0ob3K7whatCnXgT7iQhpVTqcrIK5WqbL50V47YktXApjo6//cSYGU/SceNqprfrzvBz+rOpSXPACqJOUiRl5YYmjRqwaHivyrx6WZAzMP0JhKwOqNvYlaq79EAHh+wOdGghe7nn05d5/9WBtNyxlVty/8ntuUMqgzdYwTL7AGfrrteVlFamPZzkxMOt57b7um5jVyr11ftaKE7rVAfbTTjqoG3kPPh/7P/LaqZ0OpuRPfuyvfF+1d7bsEGa43XbzdLdjpcUOlnPbXf+Zo/2LR21RymVvOp1AI+0TnXlbsLff4chQ2DIM5CdzZfP/pdx2w5mu82Iefde+52ZgXbs2euoWJWvrkq4iUi78zcjOZdTKZWc6nUAj6pO9cyZcMstUFQE//gHjBxJ9/3240vgiLzpNd7iXlZuQh68EOlGHM2BK1V31esAHlFw27wZBg6EV1+FY4+FL7+EU0+t8pJYrOUGK3i7XVVPo4926Z+WclWq7qrXk5iOJviMgbfegg4d4L//hWHDoLCwWvAG+4nO7m2DH4XWpKEr6HXraw3IykyvXFIY7bptuzZpLRSlUl+9HoEP7t0u6CaXyuC2fj3cdhvk58MJJ8Ds2dC5s+39QpVNHZa/lDfmF1FuDC4Rrj65Nd0Ob8HdkxcFvVdJaRmLhveqcR+1lKtSdVe9P9Ah6CqULofCSy9ZKZPdu+GBB6w/N6j+8y6/0MMDHyyvrE+Sme5mxMXhJxd97PLmLhFWjTq/Jl1TStUReiq9jWp1qn/+GXr1go8/hjPOgOefh2OOCfre/EIPg99eXCVXXVJaxuC3FlfeOxy7H5+xOLRBKVW31esceBXl5TBuHBx3HMyfD+PHw7x5tsEbsN0pWVZhHG9Vz7LJw9tdV0opHw3gACtWwGmnwd13w5//DMuXQ//+kBb62xNqKZ7TZXo6yaiUilb9DuB79sC//mWdivPjj/D66zB9OrRu7ejtoZbiOV2ml5uTxahLO8VkxYlSqn6pvznwggLo0weWLIGrroInn4SDbA8WCmpw73bVcuAA7jSJaASt50UqpaJR/0bgpaVwzz1w8smwaZO1RPDNNyMO3mAF3rGXd6Z5xr6j0TLT3Yy9orMGZKVU3NWvEfinn0LfvvDTT3DzzTBmDGRm1uiWOnpWSiVK/RiBb98Ot94KZ54JFRUwZw5MmFDj4K2UUolU9wP49OnQsaMVsAcOtHLePe0PQFBKqVRRdwP4pk1w3XVw4YXQtCl89RU8+ig0aZLolimlVEzUvQBujDUpeeyxMHkyDB8O331nTVoqpVQdUrcmMT0eq/jU1Klw4onw4ovQqVOiW6WUUnFRN0bgxlg1Szp0gI8+gn//G77+WoO3UqpOS/0R+KpV1pLAefOsVSbPPw9HHZXoVimlVNyl7gi8vBwee8waZS9cCP/5j7U8UIO3UqqeSM0R+LJl1jb4BQusVSbjx8NhhyW6VUopVatqNAIXkXNF5AcR+UlE8mLVKFt79liHK3TtCqtXW0ecTZ2qwVspVS9FPQIXERfwDHAO8CvwrYhMNcasiFXjqliwwBp1L1sG11wDTzwBLVvG5VFKKZUKajICPwn4yRiz2hizB3gTuCQ2zQowcqR1iPDWrfDBBzBpkgZvpVS9V5MAngUU+X3+q/daFSLST0QKRKSguLg4uie1bWutNFm+3Mp5K6WUqlEAlyDXqp0vZoyZYIzpZozp1jLaUfPVV8Nzz0GzZtG9Xyml6qCaBPBfAf+jaw4D1tWsOUoppZyqSQD/FjhaRI4QkYbAX4GpsWmWUkqpcKJehWKM2SsidwCzABcw0RizPGYtU0opFVKNNvIYYz4EPoxRW5RSSkUgdbfSK6VUPacBXCmlUpQGcKWUSlEawJVSKkWJMdX23sTvYSLFwC9Rvv1AYFMMm5NI2pfkU1f6AdqXZFWTvhxujKm2E7JWA3hNiEiBMaZbotsRC9qX5FNX+gHal2QVj75oCkUppVKUBnCllEpRqRTAJyS6ATGkfUk+daUfoH1JVjHvS8rkwJVSSlWVSiNwpZRSfjSAK6VUikqJAF7rhyfHkIhMFJGNIrLM71oLEflIRH70fmyeyDY6ISKtRWSeiHwvIstF5B/e66nYl8YiskBEFnv78oD3esr1BazzaUWkUESmeT9P1X6sEZGlIrJIRAq811K1L5ki8raIrPT+mzk1Hn1J+gDud3jyeUAH4GoR6ZDYVkXkZeDcgGt5wBxjzNHAHO/nyW4vMMgYcyxwCnC79+8hFfuyG+hpjOkMdAHOFZFTSM2+APwD+N7v81TtB0APY0wXv/XSqdqXccBMY0x7oDPW30/s+2KMSer/gFOBWX6fDwGGJLpdEfYhG1jm9/kPQCvvn1sBPyS6jVH06X3gnFTvC5ABfAecnIp9wToJaw7QE5jmvZZy/fC2dQ1wYMC1lOsL0BT4Ge8ikXj2JelH4Dg8PDnFHGyMWQ/g/XhQgtsTERHJBnKA+aRoX7xph0XARuAjY0yq9uUJ4B6gwu9aKvYDrDN1Z4vIQhHp572Win05EigGXvKmtl4QkSbEoS+pEMAdHZ6saoeI7Ae8A9xtjNme6PZEyxhTbozpgjWCPUlEjktwkyImIhcCG40xCxPdlhjpbozpipUuvV1Ezkh0g6LUAOgKjDfG5AA7iFPqJxUCeF08PHmDiLQC8H7cmOD2OCIibqzgPckY8673ckr2xccYUwJ8gjVPkWp96Q5cLCJrgDeBniLyOqnXDwCMMeu8HzcC7wEnkZp9+RX41ftbHcDbWAE95n1JhQBeFw9Pngrc6P3zjVj55KQmIgK8CHxvjHnM70up2JeWIpLp/XM6cDawkhTrizFmiDHmMGNMNta/i7nGmOtIsX4AiEgTEdnf92egF7CMFOyLMeY3oEhE2nkvnQWsIB59SXTC3+GkwPnA/4BVwNBEtyfCtr8BrAfKsH4y9wEOwJp4+tH7sUWi2+mgH6dhpa6WAIu8/52fon05Hij09mUZcL/3esr1xa9PZ7JvEjPl+oGVN17s/W+57995KvbF2+4uQIH3/7F8oHk8+qJb6ZVSKkWlQgpFKaVUEBrAlVIqRWkAV0qpFKUBXCmlUpQGcKWUSlEawJVSKkVpAFdKqRT1/1URMjbbvkvQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(yhat,y_train, 'o')\n",
    "plt.plot([0, 60], [0,60], 'r-') # prediction looks like good"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Don't look at this part -- 2 weeks later\n",
    "## Same Solution - Optimization based on Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* This part will be revisited later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\anaconda3\\envs\\env\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in true_divide\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    }
   ],
   "source": [
    "XX=(X-X.mean(axis=0))/np.std(X, axis=0)\n",
    "# GD algorithm을 쓰려면 X standardized를 해줘야 함"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# XX = X -> standardization is very important for gradient descent working\n",
    "# 안하면 overflow ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "XX[:,0]=1 # 1st column = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(506, 1)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y=y.reshape([506,1])\n",
    "Y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(beta, X, Y): # 정의 # this is -1 (beta)\n",
    "    # -f(beta) : negative value of log likelihood function : maximizing\n",
    "    # = minimize (- log function)\n",
    "    temp = np.exp(np.dot(X,beta))-Y*np.dot(X, beta)+scipy.special.gammaln(Y+1)\n",
    "    return np.sum(temp)\n",
    "\n",
    "# df obtained for free once you define f\n",
    "# df(f) \n",
    "\n",
    "def df(beta, X, Y): # don't need to calculate ... automatical\n",
    "    # in actual data analysis, using tensorflow, df is freely obtained from f function, derivative(f)\n",
    "    temp = np.exp(np.dot(X,beta))-Y\n",
    "    temp2 = np.mean(X*temp, axis=0)\n",
    "    return temp2.reshape([14,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gradient descent ... to find minimize things\n",
    "beta0=np.random.normal(0,1,[14,1]) # 14 column ... dim of beta0 = [14,1]\n",
    "lr=0.001\n",
    "iter1=5000 # iteration : 5000 -> # of data in each iteration : 506\n",
    "# of total data -> 5000 * 506\n",
    "beta=beta0\n",
    "history=np.ones([iter1,])\n",
    "for i in range(0,iter1):\n",
    "    beta = beta - lr * df(beta, XX, Y) # y is moving\n",
    "    history[i]=f(beta, XX, Y)\n",
    "    #print(\"************************\",  i, \"**********************\",    beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x201575f1348>]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZGklEQVR4nO3df4xd5X3n8ffn3ju/PGMb/xi7w9hgQ7wtBiVOmFjuUlVsvBvctClUgpUjbbFWqI4Q6RJtpApaadv9w1KpmrCiK1iRQjFsAmHzQ9AWWhAkStJSnHEgMcYxDGDA8WAPtmPGY2bGM/PdP+5zx/eO7/xg7Jk7nvN5SVfn3O85597nIY4/fs5zzrmKCMzMzHK1boCZmc0NDgQzMwMcCGZmljgQzMwMcCCYmVlSqHUDpmv58uWxZs2aWjfDzOyCsnv37vcjorXatgs2ENasWUNnZ2etm2FmdkGR9PZ423zKyMzMAAeCmZklDgQzMwMcCGZmljgQzMwMcCCYmVniQDAzMyCDgfCTA8f46jP7GRwaqXVTzMzmlMwFwu63j/M3z3cxNOJAMDMrl7lAyKm49O8CmZlVylwgiGIijDgRzMwqZC8QSiOE2jbDzGzOyVwglHiAYGZWKXOBIA8RzMyqylwgjE4qOxHMzCpkLhBSHjDiPDAzq5C9QEinjMKTCGZmFTIYCMWl48DMrFIGA6E0QqhxQ8zM5pjsBUJa+pSRmVml7AWCTxmZmVWVvUDAp4zMzKrJXiD4PgQzs6oyFwh+2qmZWXWZCwQ/7dTMrLrMBQIeIZiZVZW5QNDku5iZZVLmAiHnG9PMzKqaNBAkrZb0fUn7JO2VdHuq/4WkX0p6Ob0+V3bMnZK6JO2XdF1Z/WpJe9K2e5RuG5bUIOlbqf6ipDUz0NfUhuLScwhmZpWmMkIYAr4SEVcAm4DbJK1P2+6OiA3p9RRA2rYVuBLYAtwrKZ/2vw/YDqxLry2pfgtwPCI+BtwN3HXuXavON6aZmVU3aSBERHdE/DSt9wL7gPYJDrkeeCwiBiLiLaAL2CipDVgUES9E8bkRDwM3lB2zM61/G9hcGj2cb2duTHMkmJmV+0hzCOlUzieBF1PpS5J+LulBSUtSrR14t+ywg6nWntbH1iuOiYgh4ASwrMr3b5fUKamzp6fnozS97DOKS8eBmVmlKQeCpBbgO8CXI+IDiqd/Lgc2AN3AV0u7Vjk8JqhPdExlIeL+iOiIiI7W1tapNr2Cn3ZqZlbdlAJBUh3FMPhGRHwXICIOR8RwRIwAXwc2pt0PAqvLDl8FHEr1VVXqFcdIKgCLgWPT6dCkfUlLnzIyM6s0lauMBDwA7IuIr5XV28p2+wPglbT+JLA1XTm0luLk8a6I6AZ6JW1Kn3kz8ETZMdvS+o3A8zFDf2P7lJGZWXWFKexzDfCHwB5JL6fanwJfkLSB4t+tB4AvAkTEXkmPA69SvELptogYTsfdCjwENAFPpxcUA+cRSV0URwZbz6VTE/HTTs3Mqps0ECLix1Q/x//UBMfsAHZUqXcCV1Wp9wM3TdaW8yHnp52amVWVuTuVR29MG6ltO8zM5prMBUJpsOMRgplZpcwFgvy0UzOzqrIXCLVugJnZHJW5QPDTTs3MqstcIPhpp2Zm1WU2EBwHZmaVshcIftqpmVlV2QsEjxDMzKrKYCB4hGBmVk32AiEtnQdmZpWyFwg+ZWRmVlX2AsFPOzUzqypzgTD6tFMngplZhcwFAqM3ptW2GWZmc03mAkF+2qmZWVXZC4TRy4xq2gwzszknc4Ew+nC7GrfDzGyuyVwg+OF2ZmbVZS8Q0tJ5YGZWKXuB4BvTzMyqylwg4KedmplVlblAyKc70zyHYGZWKXuBkM4ZDY/UuCFmZnNM9gIhVwoEjxDMzMo5EMzMDMhyIHgOwcysQnYDYcSTCGZm5SYNBEmrJX1f0j5JeyXdnupLJT0r6fW0XFJ2zJ2SuiTtl3RdWf1qSXvStnuUfs9SUoOkb6X6i5LWzEBfASikQBga9gjBzKzcVEYIQ8BXIuIKYBNwm6T1wB3AcxGxDnguvSdt2wpcCWwB7pWUT591H7AdWJdeW1L9FuB4RHwMuBu46zz0rSpfdmpmVt2kgRAR3RHx07TeC+wD2oHrgZ1pt53ADWn9euCxiBiIiLeALmCjpDZgUUS8EMW7wh4ec0zps74NbC6NHs63UiAMeVLZzKzCR5pDSKdyPgm8CKyMiG4ohgawIu3WDrxbdtjBVGtP62PrFcdExBBwAlhW5fu3S+qU1NnT0/NRmj5qdITgQDAzqzDlQJDUAnwH+HJEfDDRrlVqMUF9omMqCxH3R0RHRHS0trZO1uSqSjemeYRgZlZpSoEgqY5iGHwjIr6byofTaSDS8kiqHwRWlx2+CjiU6quq1CuOkVQAFgPHPmpnpiKf930IZmbVTOUqIwEPAPsi4mtlm54EtqX1bcATZfWt6cqhtRQnj3el00q9kjalz7x5zDGlz7oReD5m6OlzZx5d4UAwMytXmMI+1wB/COyR9HKq/Snwl8Djkm4B3gFuAoiIvZIeB16leIXSbRExnI67FXgIaAKeTi8oBs4jkroojgy2nlu3xucb08zMqps0ECLix1Q/xw+weZxjdgA7qtQ7gauq1PtJgTLTRgPB9yGYmVXI3p3KnlQ2M6sqc4GQy4mcfGOamdlYmQsEKJ428gjBzKxSZgPBN6aZmVXKZiDIIwQzs7GyGQg5+T4EM7MxHAhmZgZkNhByvjHNzGyMjAaCb0wzMxsrk4FQ8AjBzOwsmQyEXM4PtzMzGyuTgVDI5XzZqZnZGJkMhLq8OD00UutmmJnNKZkMhPpCjtPDDgQzs3KZDIS6fI5BB4KZWYXsBoJPGZmZVchkIDQUPEIwMxsrk4FQl/ccgpnZWJkMhPp8jtNDvuzUzKxcJgOhzqeMzMzOkslAqPekspnZWbIZCAV5hGBmNkYmA8GTymZmZ8tkIBQnlR0IZmblMhkInlQ2MztbJgOhPp/j9HAQ/k0EM7NR2QyEQrHbHiWYmZ2RyUCoywuA0/4ZTTOzUZMGgqQHJR2R9EpZ7S8k/VLSy+n1ubJtd0rqkrRf0nVl9asl7Unb7pGkVG+Q9K1Uf1HSmvPcx7PU59MIwRPLZmajpjJCeAjYUqV+d0RsSK+nACStB7YCV6Zj7pWUT/vfB2wH1qVX6TNvAY5HxMeAu4G7ptmXKatLp4x86amZ2RmTBkJE/BA4NsXPux54LCIGIuItoAvYKKkNWBQRL0RxJvdh4IayY3am9W8Dm0ujh5niEYKZ2dnOZQ7hS5J+nk4pLUm1duDdsn0Oplp7Wh9brzgmIoaAE8Cyal8oabukTkmdPT090264J5XNzM423UC4D7gc2AB0A19N9Wr/so8J6hMdc3Yx4v6I6IiIjtbW1o/U4HKlEYJPGZmZnTGtQIiIwxExHBEjwNeBjWnTQWB12a6rgEOpvqpKveIYSQVgMVM/RTUtdSkQBk47EMzMSqYVCGlOoOQPgNIVSE8CW9OVQ2spTh7viohuoFfSpjQ/cDPwRNkx29L6jcDzMcN3jDXVF+e5+08Pz+TXmJldUAqT7SDpUeBaYLmkg8CfA9dK2kDx1M4B4IsAEbFX0uPAq8AQcFtElP7WvZXiFUtNwNPpBfAA8IikLoojg63noV8Taqwr5mC/J5XNzEZNGggR8YUq5Qcm2H8HsKNKvRO4qkq9H7hpsnacT411xRHCh4MeIZiZlWTyTuVSIAwMORDMzEoyGQhNHiGYmZ0lk4FQGiF4UtnM7IyMBoInlc3MxspmIBR8ysjMbKxMBkIuJ+oLOfo9qWxmNiqTgQDFieV+jxDMzEZlNhAa63L0+9EVZmajMhsITXV5PvRVRmZmozIbCI11eV92amZWJtOB4BGCmdkZGQ6EnB9/bWZWJrOB4DkEM7NKmQ0EzyGYmVXKbCB4hGBmVimzgdBQl/d9CGZmZTIbCAvq83w4OFTrZpiZzRmZDYTm+jx9g8OMjMzozzebmV0wshsIDcVfD/U8gplZUWYDYUEKhL4BnzYyM4MMB0JLQ/E3Efr8xFMzMyDDgdBc7xGCmVm57AZCOmV00oFgZgY4EDjlS0/NzIAMB0JpDuHkgOcQzMwgw4GwIM0hnPIpIzMzIMOB4DkEM7NK2Q2E+nTZqU8ZmZkBUwgESQ9KOiLplbLaUknPSno9LZeUbbtTUpek/ZKuK6tfLWlP2naPJKV6g6RvpfqLktac5z5WVcjnaCjkPKlsZpZMZYTwELBlTO0O4LmIWAc8l94jaT2wFbgyHXOvpHw65j5gO7AuvUqfeQtwPCI+BtwN3DXdznxULQ0FnzIyM0smDYSI+CFwbEz5emBnWt8J3FBWfywiBiLiLaAL2CipDVgUES9ERAAPjzmm9FnfBjaXRg8zrbmh4BvTzMyS6c4hrIyIboC0XJHq7cC7ZfsdTLX2tD62XnFMRAwBJ4Bl1b5U0nZJnZI6e3p6ptn0MxakJ56amdn5n1Su9i/7mKA+0TFnFyPuj4iOiOhobW2dZhPPaPEIwcxs1HQD4XA6DURaHkn1g8Dqsv1WAYdSfVWVesUxkgrAYs4+RTUjmj2HYGY2arqB8CSwLa1vA54oq29NVw6tpTh5vCudVuqVtCnND9w85pjSZ90IPJ/mGWbcoqY6evsdCGZmAIXJdpD0KHAtsFzSQeDPgb8EHpd0C/AOcBNAROyV9DjwKjAE3BYRpZP0t1K8YqkJeDq9AB4AHpHURXFksPW89GwKFjcVOPHh6dn6OjOzOW3SQIiIL4yzafM4++8AdlSpdwJXVan3kwJlti1qrOPEh6eJCGbpwiYzszkrs3cqAyxuqmN4JDjlK43MzLIdCIua6gB82sjMjIwHwuIUCB/0OxDMzDIdCIsa0wjhlAPBzCzTgXBmhOBLT83MMh0Ii5qKF1l5DsHMLOOBMDpCcCCYmWU7EBY2+iojM7OSTAdCPicWNhR8lZGZGRkPBCjei+ARgpmZA6EYCL7s1MzMgbCsuZ5jpwZr3Qwzs5pzILTUc/SkA8HMzIHQ3MDRkwO1boaZWc05EFrq6Rsc5kM/8dTMMi7zgbC8pR6Ao30eJZhZtmU+EJY1NwB4HsHMMi/zgbDUIwQzM8CBwHKPEMzMAAcCy0ZHCA4EM8u2zAfCgvo8jXU5enp9ysjMsi3zgSCJixc38d6J/lo3xcyspjIfCABtFzXyy199WOtmmJnVlAMBuHhxE90nHAhmlm0OBKDtoiaO9A5wenik1k0xM6sZBwLQflEjEXgewcwyzYEAtC1uAqDbgWBmGeZAANqXFAPh4PFTNW6JmVntnFMgSDogaY+klyV1ptpSSc9Kej0tl5Ttf6ekLkn7JV1XVr86fU6XpHsk6Vza9VFdsnQBhZx4o+fkbH6tmdmccj5GCP8hIjZEREd6fwfwXESsA55L75G0HtgKXAlsAe6VlE/H3AdsB9al15bz0K4pq8vnuHTZAt440jebX2tmNqfMxCmj64GdaX0ncENZ/bGIGIiIt4AuYKOkNmBRRLwQEQE8XHbMrLm8tcUjBDPLtHMNhACekbRb0vZUWxkR3QBpuSLV24F3y449mGrtaX1s/SyStkvqlNTZ09Nzjk2vdPmKFg4c7WPIl56aWUadayBcExGfAn4HuE3Sb0+wb7V5gZigfnYx4v6I6IiIjtbW1o/e2gmsW9HC6eHgzfd92sjMsumcAiEiDqXlEeB7wEbgcDoNRFoeSbsfBFaXHb4KOJTqq6rUZ9XHV10EwMvv/mq2v9rMbE6YdiBIapa0sLQOfBZ4BXgS2JZ22wY8kdafBLZKapC0luLk8a50WqlX0qZ0ddHNZcfMmsuWN7OoscBL7/xqtr/azGxOKJzDsSuB76UrRAvANyPinyT9BHhc0i3AO8BNABGxV9LjwKvAEHBbRJR+2f5W4CGgCXg6vWZVLic2XLKEl945PttfbWY2J0w7ECLiTeATVepHgc3jHLMD2FGl3glcNd22nC8b1yzhr595jSO9/axY2Fjr5piZzSrfqVzmM7+xEoDv/+LIJHuamc0/DoQyV7QtpP2iJp7Ze7jWTTEzm3UOhDKS+P0NF/OD13r85FMzyxwHwhhf+PQlDI8E33zx7Vo3xcxsVjkQxrhk2QKuu3IlD/7LAY73Dda6OWZms8aBUMVXPvvr9A0O8bVnX6t1U8zMZo0DoYp/t3Ih//Xfr+WRf3ubf+l6v9bNMTObFQ6EcfzJll/nstZm/tujL/HuMf9wjpnNfw6EcTTW5fn6zR0MjQTb/m4Xhz/wVUdmNr85ECZweWsLX7+5g/dO9HPj//lXDvhJqGY2jzkQJrFx7VK++Ueb6O0f4vN/82Oe3tNd6yaZmc0IB8IUbFh9EX//pd/ishUt3PqNn/LHj77kU0hmNu84EKZo9dIF/L8v/ia3b17HP+99j8/89Q+4659+wfsnB2rdNDOz80LFnzG+8HR0dERnZ2dNvvvto3381T/v56k93TQUcnz+4xfznz+9mo5Ll5AeB25mNidJ2h0RHVW3ORCmr+vISf72R2/y9z87RN/gMJcsXcB/Wr+SzVes4NNrllKX9wDMzOYWB8IM6xsY4h/3dPOPP+/mhTeOMjg8QktDgU9duoRPX7qEjjVL+cTqxSyoP5ffIzIzO3cOhFnUNzDEj15/nx+93sPut4+z/3AvESDBmmXNXNG2kPVti7iibRGXtbawakmTRxJmNmsmCgT/k/U8a24osOWqX2PLVb8GwIlTp9n9zjH2HPyAfd0f8MovP+CpPe+N7p/PidVLmrh0WTNrlzdz6bIFtC1uom1xI20XNbK8uYFczvMSZjbzHAgzbPGCOj7zGytHf40NoLf/NK8d7uXNnj4OHO3jwNFTHHi/j84Dx+gbHK44vpATKxc10ra4kZWLG1nWXM/S5vq0bCiutxRrFzXVUfBow8ymyYFQAwsb67j60qVcfenSinpEcLRvkPdO9NN9op/3TnyYlsX3+w59wNG+QU58eHrcz15Qn6eloUBLY4GFo8s6WhoLtDQUaG7I01DI01iXq7psKC0LOQp5UciJfC6Xlipb5sjnz7zPSx7JmF3gHAhziCSWtzSwvKWBq9oXj7vf6eERjp8a5FjfIMdODnK0r7h+/NQgfQND9PYP0TswxMn+IU4ODNHTe5KTqXZqcJjhkZmZN5IgJ6G0LgSCXFov1or9FEDZ+5zO1JU2ntm//HOnFzoTHTbhNsbfOPFxE33fBJ85wXETbZzu99mF6fbN6/j8Jy4+75/rQLgA1eVzrFjYyIqFjdM6fmh4hP6hEQZOD59Znh5hYOjMcmBohJGRYGgkGB5djpx5P1y9PhJBBASkZfFNtXrpeoaIIDh7H9I+pf0nyrGJro0ofvq4G6eziYkuxpj4uNn9vok32oVqcVPdjHyuAyGDCvkcLfkcLQ3+n9/MzvAMpJmZAQ4EMzNLHAhmZgY4EMzMLHEgmJkZ4EAwM7PEgWBmZoADwczMkgv28deSeoC3p3n4cuD989icC4H7nA3uczacS58vjYjWahsu2EA4F5I6x3se+HzlPmeD+5wNM9VnnzIyMzPAgWBmZklWA+H+WjegBtznbHCfs2FG+pzJOQQzMztbVkcIZmY2hgPBzMyADAaCpC2S9kvqknRHrdtzLiQ9KOmIpFfKakslPSvp9bRcUrbtztTv/ZKuK6tfLWlP2naP5uhvLkpaLen7kvZJ2ivp9lSfz31ulLRL0s9Sn/9nqs/bPpdIykt6SdI/pPfzus+SDqS2viypM9Vmt88RkZkXkAfeAC4D6oGfAetr3a5z6M9vA58CXimr/RVwR1q/A7grra9P/W0A1qb/Dvm0bRfwmxR/mvdp4Hdq3bdx+tsGfCqtLwReS/2az30W0JLW64AXgU3zuc9lff/vwDeBf5jvf7ZTWw8Ay8fUZrXPWRshbAS6IuLNiBgEHgOur3Gbpi0ifggcG1O+HtiZ1ncCN5TVH4uIgYh4C+gCNkpqAxZFxAtR/NP0cNkxc0pEdEfET9N6L7APaGd+9zki4mR6W5dewTzuM4CkVcDvAn9bVp7XfR7HrPY5a4HQDrxb9v5gqs0nKyOiG4p/gQIrUn28vren9bH1OU3SGuCTFP/FPK/7nE6dvAwcAZ6NiHnfZ+B/AX8CjJTV5nufA3hG0m5J21NtVvuctV9Zr3YuLSvX3Y7X9wvuv4mkFuA7wJcj4oMJTpHOiz5HxDCwQdJFwPckXTXB7hd8nyX9HnAkInZLunYqh1SpXVB9Tq6JiEOSVgDPSvrFBPvOSJ+zNkI4CKwue78KOFSjtsyUw2nYSFoeSfXx+n4wrY+tz0mS6iiGwTci4rupPK/7XBIRvwJ+AGxhfvf5GuD3JR2geFr3M5L+L/O7z0TEobQ8AnyP4inuWe1z1gLhJ8A6SWsl1QNbgSdr3Kbz7UlgW1rfBjxRVt8qqUHSWmAdsCsNQ3slbUpXI9xcdsycktr3ALAvIr5Wtmk+97k1jQyQ1AT8R+AXzOM+R8SdEbEqItZQ/P/o8xHxX5jHfZbULGlhaR34LPAKs93nWs+sz/YL+BzFq1PeAP6s1u05x748CnQDpyn+y+AWYBnwHPB6Wi4t2//PUr/3U3blAdCR/vC9Afxv0h3sc+0F/BbF4e/PgZfT63PzvM8fB15KfX4F+B+pPm/7PKb/13LmKqN522eKVz7+LL32lv5umu0++9EVZmYGZO+UkZmZjcOBYGZmgAPBzMwSB4KZmQEOBDMzSxwIZmYGOBDMzCz5/4dWHXqcNC6NAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(history[:]) # learning is over\n",
    "# negative function의 값이 점점 작아진다는 것은 곧 maximize를 의미함\n",
    "# learning curve - flat : learning is over when 1000\n",
    "# x축 : iteration number\n",
    "# negative라서 값이 작을 수록 좋음"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1eedf236908>]"
      ]
     },
     "execution_count": 205,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAf30lEQVR4nO3dfZAd1Xnn8e9v7rzqdfQykoWkWMKWMbLXlvGs0BapxAZbCJyNcC3egqSCykuVEhtq7XW8MSSpJX6hymwqxiHBTuGgRWw5FhjbQcuKyFrA6/WuERqMEAgZayxkIyRLA3p/GWlenv2jz0hXozszV6OZuaPp36eq63Y/fbrvOarRPHNOn+5WRGBmZvlWVekKmJlZ5TkZmJmZk4GZmTkZmJkZTgZmZgZUV7oCgzV9+vSYN29epathZnZRef7559+MiKbe8Ys2GcybN4+WlpZKV8PM7KIi6Vel4h4mMjOz8pOBpIKkFyQ9kbbnS9ooabukRyTVpnhd2m5N++cVnePOFH9V0rVF8WUp1irpjqFrnpmZleN8egafAbYVbd8D3BsRC4ADwK0pfitwICLeCdybyiFpIXAT8B5gGfCNlGAKwP3AdcBC4OZU1szMRkhZyUDSHOBjwD+mbQFXA4+lIquBG9L68rRN2n9NKr8cWBMRJyPiNaAVWJyW1ojYERGngDWprJmZjZByewZfB/4M6E7b04CDEdGZtncBs9P6bOB1gLT/UCp/Ot7rmL7i55C0UlKLpJa2trYyq25mZgMZMBlI+j1gX0Q8XxwuUTQG2He+8XODEQ9ERHNENDc1nTMzyszMBqmcqaVXAb8v6XqgHphE1lNolFSd/vqfA+xO5XcBc4FdkqqBycD+oniP4mP6ipuZ2QgYsGcQEXdGxJyImEd2AfjpiPhD4BngxlRsBfB4Wl+btkn7n47sOdlrgZvSbKP5wALgOWATsCDNTqpN37F2SFpXwkP/9zX+x4vONWZmxS7kPoMvAJ+T1Ep2TeDBFH8QmJbinwPuAIiIrcCjwCvAvwC3RURX6lncDqwnm630aCo7LP7puV+z7qU9w3V6M7OL0nndgRwRPwJ+lNZ3kM0E6l2mHfhEH8ffDdxdIr4OWHc+dRmsmkIVHV3dAxc0M8uR3N2BXF2o4lSX3+5mZlYsd8mgtiA6Ot0zMDMrlrtkUFOoorPbycDMrFguk4GHiczMzpbDZOBhIjOz3nKYDDybyMyst1wmg85uDxOZmRXLXTKoLohTHiYyMztL7pJBrYeJzMzOkbtk4GEiM7Nz5TIZeDaRmdnZcpgMxEkPE5mZnSV3yaCupsCpzm66PVRkZnZa7pLBuNoCAO2dXRWuiZnZ6JG7ZNBQkyWD46ecDMzMeuQvGaSewQknAzOz03KXDHqGiU50OBmYmfUYMBlIqpf0nKQXJW2V9MUUf0jSa5I2p2VRikvSfZJaJW2RdEXRuVZI2p6WFUXxD0p6KR1znyQNR2PhTDLwMJGZ2RnlvPbyJHB1RByVVAP8RNKTad9/jojHepW/juxl9wuAK4FvAldKmgrcBTQDATwvaW1EHEhlVgLPkr3+chnwJMOgoSZrsoeJzMzOGLBnEJmjabMmLf3Ny1wOPJyOexZolDQLuBbYEBH7UwLYACxL+yZFxE8jIoCHgRsuoE39On3NoKNzuL7CzOyiU9Y1A0kFSZuBfWS/0DemXXenoaB7JdWl2Gzg9aLDd6VYf/FdJeKl6rFSUouklra2tnKqfg4PE5mZnausZBARXRGxCJgDLJb0XuBO4N3AvwamAl9IxUuN98cg4qXq8UBENEdEc1NTUzlVP4enlpqZneu8ZhNFxEHgR8CyiNiThoJOAv8NWJyK7QLmFh02B9g9QHxOifiwOH3TmWcTmZmdVs5soiZJjWm9AfgI8PM01k+a+XMD8HI6ZC1wS5pVtAQ4FBF7gPXAUklTJE0BlgLr074jkpakc90CPD60zTyjwcNEZmbnKGc20SxgtaQCWfJ4NCKekPS0pCayYZ7NwJ+k8uuA64FW4DjwSYCI2C/py8CmVO5LEbE/rX8KeAhoIJtFNCwziQDqq50MzMx6GzAZRMQW4AMl4lf3UT6A2/rYtwpYVSLeArx3oLoMhaoq0VBT4MQpzyYyM+uRuzuQIbtu4DuQzczOyGUyqK8peJjIzKxILpPBhLpqjp30MJGZWY9cJoPxdQWOnXTPwMysR06TQTVH3TMwMzstl8nAw0RmZmfLZTIY72RgZnaWXCaDCR4mMjM7Sy6Twfi6AsdOdZHdH2dmZjlNBtV0dQcnO7srXRUzs1Ehl8lgQl32FI4j7R4qMjODnCaD8bVZMvBFZDOzTC6TwYT6LBn4IrKZWSafyaDOPQMzs2K5TAbje5KBH2NtZgbkNBlMqMtecHPUzycyMwNymgzGe5jIzOws5bwDuV7Sc5JelLRV0hdTfL6kjZK2S3pEUm2K16Xt1rR/XtG57kzxVyVdWxRflmKtku4Y+maezcnAzOxs5fQMTgJXR8T7gUXAsvSi+3uAeyNiAXAAuDWVvxU4EBHvBO5N5ZC0ELgJeA+wDPiGpEJ6t/L9wHXAQuDmVHbY9Ewt9X0GZmaZAZNBZI6mzZq0BHA18FiKrwZuSOvL0zZp/zWSlOJrIuJkRLwGtAKL09IaETsi4hSwJpUdNoX0HmT3DMzMMmVdM0h/wW8G9gEbgF8CByOi57fpLmB2Wp8NvA6Q9h8CphXHex3TV3xYja+r9mwiM7OkrGQQEV0RsQiYQ/aX/OWliqVP9bHvfOPnkLRSUouklra2toEr3o8JdQXPJjIzS85rNlFEHAR+BCwBGiVVp11zgN1pfRcwFyDtnwzsL473OqaveKnvfyAimiOiuamp6Xyqfg6/08DM7IxyZhM1SWpM6w3AR4BtwDPAjanYCuDxtL42bZP2Px3Zs6LXAjel2UbzgQXAc8AmYEGanVRLdpF57VA0rj9+9aWZ2RnVAxdhFrA6zfqpAh6NiCckvQKskfQV4AXgwVT+QeC/S2ol6xHcBBARWyU9CrwCdAK3RUQXgKTbgfVAAVgVEVuHrIV9mFBXzd7D7cP9NWZmF4UBk0FEbAE+UCK+g+z6Qe94O/CJPs51N3B3ifg6YF0Z9R0yE+qq2eGegZkZkNM7kKFnmMgXkM3MIMfJIJtN1FHpapiZjQq5TQbj66pp7+ims8uvvjQzy20yOPNOAw8VmZnlNhlMqq8B4HC7h4rMzHKbDCbW+2F1ZmY9cpwMsp7BEfcMzMzynAzcMzAz6+Fk4OmlZmZ5TgY9w0TuGZiZ5TgZeJjIzKxHbpNBfU2B2kKVp5aamZHjZABZ78A9AzMzJwMnAzMzcp8ManyfgZkZuU8G7hmYmYGTgXsGZmaU9w7kuZKekbRN0lZJn0nxv5L0hqTNabm+6Jg7JbVKelXStUXxZSnWKumOovh8SRslbZf0SHoX8rDLhoncMzAzK6dn0An8aURcDiwBbpO0MO27NyIWpWUdQNp3E/AeYBnwDUmF9A7l+4HrgIXAzUXnuSedawFwALh1iNrXr0lOBmZmQBnJICL2RMTP0voRYBswu59DlgNrIuJkRLwGtJK9K3kx0BoROyLiFLAGWC5JwNXAY+n41cANg23Q+ZhYX83Rk510dcdIfJ2Z2ah1XtcMJM0DPgBsTKHbJW2RtErSlBSbDbxedNiuFOsrPg04GBGdveLDrucu5KMn3Tsws3wrOxlImgB8D/hsRBwGvgm8A1gE7AH+pqdoicNjEPFSdVgpqUVSS1tbW7lV79MkP8bazAwoMxlIqiFLBN+OiO8DRMTeiOiKiG7gW2TDQJD9ZT+36PA5wO5+4m8CjZKqe8XPEREPRERzRDQ3NTWVU/V++flEZmaZcmYTCXgQ2BYRXyuKzyoq9nHg5bS+FrhJUp2k+cAC4DlgE7AgzRyqJbvIvDYiAngGuDEdvwJ4/MKaVR4/udTMLFM9cBGuAv4IeEnS5hT7c7LZQIvIhnR2An8MEBFbJT0KvEI2E+m2iOgCkHQ7sB4oAKsiYms63xeANZK+ArxAlnyGXU/P4PAJDxOZWb4NmAwi4ieUHtdf188xdwN3l4ivK3VcROzgzDDTiPELbszMMjm/A9nDRGZmkPtk4AvIZmaQ82TgF9yYmWVynQzATy41MwMnAyY31Hg2kZnlXu6TwaSGGg45GZhZzuU+GTSOczIwM8t9MpjsnoGZmZNBY0MNB487GZhZvuU+GUxuqOFwewfdfqeBmeWYk8G4WiJ845mZ5ZuTQUP2SIqDJ05VuCZmZpWT+2TQmJKBLyKbWZ45GYxLPQNfRDazHMt9MpjsnoGZmZPB5J6egZOBmeWYk0FPz+C4LyCbWX7lPhnUVRdoqCl4mMjMcm3AZCBprqRnJG2TtFXSZ1J8qqQNkranzykpLkn3SWqVtEXSFUXnWpHKb5e0oij+QUkvpWPuk1TqNZvDZrLvQjaznCunZ9AJ/GlEXA4sAW6TtBC4A3gqIhYAT6VtgOuABWlZCXwTsuQB3AVcSfa+47t6Ekgqs7LouGUX3rTy+WF1ZpZ3AyaDiNgTET9L60eAbcBsYDmwOhVbDdyQ1pcDD0fmWaBR0izgWmBDROyPiAPABmBZ2jcpIn4aEQE8XHSuETG5ocYXkM0s187rmoGkecAHgI3AzIjYA1nCAGakYrOB14sO25Vi/cV3lYiX+v6VkloktbS1tZ1P1fvlF9yYWd6VnQwkTQC+B3w2Ig73V7RELAYRPzcY8UBENEdEc1NT00BVLlvjOF8zMLN8KysZSKohSwTfjojvp/DeNMRD+tyX4ruAuUWHzwF2DxCfUyI+YrJhIk8tNbP8Kmc2kYAHgW0R8bWiXWuBnhlBK4DHi+K3pFlFS4BDaRhpPbBU0pR04XgpsD7tOyJpSfquW4rONSIax9XS3tFNe0fXSH6tmdmoUV1GmauAPwJekrQ5xf4c+CrwqKRbgV8Dn0j71gHXA63AceCTABGxX9KXgU2p3JciYn9a/xTwENAAPJmWEVP8SIr6msJIfrWZ2agwYDKIiJ9Qelwf4JoS5QO4rY9zrQJWlYi3AO8dqC7DZer4WgD2HzvFzEn1laqGmVnF5P4OZDg7GZiZ5ZGTATAtJYO3nAzMLKecDCjqGRw9WeGamJlVhpMB2WwiCfb7XgMzyyknA6BQJRobath/zD0DM8snJ4Nk6vhaX0A2s9xyMkimja/jraNOBmaWT04GiXsGZpZnTgbJ1Am1HPCrL80sp5wMkqnjajlwvIPu7pIPTDUzG9OcDJKp42vp6g6/8czMcsnJIJk2wXchm1l+ORkkfj6RmeWZk0HiZGBmeeZkkEwbXwc4GZhZPjkZJD09g7YjfiSFmeWPk0FSW13F1PG17DvSXumqmJmNuHLegbxK0j5JLxfF/krSG5I2p+X6on13SmqV9Kqka4viy1KsVdIdRfH5kjZK2i7pEUm1Q9nA8zFjYh373DMwsxwqp2fwELCsRPzeiFiUlnUAkhYCNwHvScd8Q1JBUgG4H7gOWAjcnMoC3JPOtQA4ANx6IQ26EE1OBmaWUwMmg4j4MbB/oHLJcmBNRJyMiNeAVmBxWlojYkdEnALWAMslCbgaeCwdvxq44TzbMGRmTKyn7bCHicwsfy7kmsHtkrakYaQpKTYbeL2ozK4U6ys+DTgYEZ294iVJWimpRVJLW1vbBVS9tBmT6mg7epIIP5LCzPJlsMngm8A7gEXAHuBvUlwlysYg4iVFxAMR0RwRzU1NTedX4zLMmFhHR1dwwG88M7OcGVQyiIi9EdEVEd3At8iGgSD7y35uUdE5wO5+4m8CjZKqe8Uromlidq+BZxSZWd4MKhlImlW0+XGgZ6bRWuAmSXWS5gMLgOeATcCCNHOoluwi89rIxmOeAW5Mx68AHh9MnYbCjIn1gO81MLP8qR6ogKTvAB8CpkvaBdwFfEjSIrIhnZ3AHwNExFZJjwKvAJ3AbRHRlc5zO7AeKACrImJr+oovAGskfQV4AXhwyFp3nmb09AwOOxmYWb4MmAwi4uYS4T5/YUfE3cDdJeLrgHUl4js4M8xUUTMm9QwTORmYWb74DuQi42qrmVBX7WsGZpY7Tga9+C5kM8sjJ4Ne3ja5nt8ccs/AzPLFyaCXSxob2H3wRKWrYWY2opwMermksYG9h9vp6OqudFXMzEaMk0Evsxvr6Q7Y62cUmVmOOBn0ckljAwC7DzoZmFl+OBn00pMM3jh4vMI1MTMbOU4GvVwy2T0DM8sfJ4NeGmoLTB1fyxueUWRmOeJkUMIljfWeXmpmueJkUMIlk32vgZnli5NBCZc0NvDGgRN+45mZ5YaTQQm/NXUcx051sf/YqUpXxcxsRDgZlDBv+jgAdr51rMI1MTMbGU4GJcybNh6AnW/6XgMzywcngxLmTBlHleBX7hmYWU4MmAwkrZK0T9LLRbGpkjZI2p4+p6S4JN0nqVXSFklXFB2zIpXfLmlFUfyDkl5Kx9wnSUPdyPNVW13FnCnjeO0t9wzMLB/K6Rk8BCzrFbsDeCoiFgBPpW2A64AFaVkJfBOy5EH27uQryV5xeVdPAkllVhYd1/u7KuLt08a5Z2BmuTFgMoiIHwP7e4WXA6vT+mrghqL4w5F5FmiUNAu4FtgQEfsj4gCwAViW9k2KiJ9GNo/z4aJzVdT86eN57c1jnl5qZrkw2GsGMyNiD0D6nJHis4HXi8rtSrH+4rtKxEuStFJSi6SWtra2QVa9PG+fNp4j7Z0cON4xrN9jZjYaDPUF5FLj/TGIeEkR8UBENEdEc1NT0yCrWJ5Lp2czina0HR3W7zEzGw0Gmwz2piEe0ue+FN8FzC0qNwfYPUB8Tol4xS2YOQGAX+x1MjCzsW+wyWAt0DMjaAXweFH8ljSraAlwKA0jrQeWSpqSLhwvBdanfUckLUmziG4pOldFzW5sYHxtgV/sPVLpqpiZDbvqgQpI+g7wIWC6pF1ks4K+Cjwq6Vbg18AnUvF1wPVAK3Ac+CRAROyX9GVgUyr3pYjouSj9KbIZSw3Ak2mpOEm8620TefU3TgZmNvYNmAwi4uY+dl1TomwAt/VxnlXAqhLxFuC9A9WjEi6bOZENr+ytdDXMzIad70Dux7tmTuStY6d48+jJSlfFzGxYORn047K3TQTgFx4qMrMxzsmgH++amSWDV/YcrnBNzMyGl5NBP5om1jFrcj1bdh2qdFXMzIaVk8EA3jdnMlt2Hax0NczMhpWTwQDeN6eRnW8d5+Bxv/XMzMYuJ4MBLJrbCOChIjMb05wMBvDe2ZMBPFRkZmOak8EAJjfU8I6m8Tz/qwOVroqZ2bBxMijDlZdOY9POA3R2dVe6KmZmw8LJoAz/5tJpHD3ZyUtv+LqBmY1NTgZlWHLpNAB+uuOtCtfEzGx4OBmUoWliHe+aOYGf/tLJwMzGJieDMl31zulsfG0/x052VroqZmZDzsmgTB9dOJNTnd38n+3D++5lM7NKcDIo0+J5U5ncUMMP/X4DMxuDnAzKVF2o4pp3z+Cpbfvo8BRTMxtjLigZSNop6SVJmyW1pNhUSRskbU+fU1Jcku6T1Cppi6Qris6zIpXfLmlFX99Xadf9q1kcOtHB/37VQ0VmNrYMRc/gwxGxKCKa0/YdwFMRsQB4Km0DXAcsSMtK4JuQJQ+y9ypfCSwG7upJIKPNhy5rYvqEWr77/OuVroqZ2ZAajmGi5cDqtL4auKEo/nBkngUaJc0CrgU2RMT+iDgAbACWDUO9LlhNoYqPf2A2T23b51dhmtmYcqHJIIAfSnpe0soUmxkRewDS54wUnw0U/0m9K8X6ip9D0kpJLZJa2toqM1Tz75vn0tkdfGfjryvy/WZmw+FCk8FVEXEF2RDQbZJ+p5+yKhGLfuLnBiMeiIjmiGhuamo6/9oOgQUzJ/Lhy5p46P/tpL2jqyJ1MDMbaheUDCJid/rcB/yAbMx/bxr+IX3uS8V3AXOLDp8D7O4nPmr9ye++g7eOneKRTb52YGZjw6CTgaTxkib2rANLgZeBtUDPjKAVwONpfS1wS5pVtAQ4lIaR1gNLJU1JF46XptiotXj+VBbPn8rfPb2dI+0dla6OmdkFu5CewUzgJ5JeBJ4D/mdE/AvwVeCjkrYDH03bAOuAHUAr8C3g0wARsR/4MrApLV9KsVFLEn/5sct58+gp7n/ml5WujpnZBVNEyeH5Ua+5uTlaWloqWofPf/dFHt/8Bj/49FWn34hmZjaaSXq+6FaA03wH8gX4y49dzpRxtfynRzb7YrKZXdScDC5A47ha/voT72f7vqPc8b0tXKy9LDMzJ4ML9LvvauLzS9/FP2/ezdf/1/ZKV8fMbFCqK12BseC2D7+TnW8d52+fypLBZz+yAKnU7RNmZqOTk8EQkMQ9/+59CPjbp7bzm0PtfHH5e6ivKVS6amZmZXEyGCKFqiwhzJxUz98/08rWPYf46xvfz+WzJlW6amZmA/I1gyFUVSU+f+1lfOuWZvYcbOff/t1P+MoTr/CWH2pnZqOc7zMYJgeOneKrT/6cR59/nfrqAjcv/i3+4Mq5vHPGxEpXzcxyrK/7DJwMhlnrvqP8/dPbeWLLHjq7g0VzG1n6npl8+LIZvPttE32h2cxGlJNBhb159CQ/+Nkb/PPmN9i6+zAA0yfUsmhuI4vmNvLe2ZN5R9MEZjc2UFXlBGFmw8PJYBTZe7idZ36+j+d27mfz6wfZ0Xbs9L666irmTx/PnCnjmDmpjrdNqmfmpHpmTKpjyrhaJjfUMLmhhkkNNRScNMzsPDkZjGKHjnfw898cZsebx9jRdpQdbcd44+AJ9h5u58Dxvp+KOrGumkkNNUysr6ahtkBDTbbUp/Vx6bO+pkBdTRU1VVXUFER1oYraQhXVBVFTyGI1hSqqi9ZrClVUV4lClaiSKFRlU2gLymISp/dlCynecwyn9/Vse0jMrPL6SgaeWjoKTB5Xw5WXTuPKS6eds6+9o4u2IyfZe7idg8c7OHTi7OXwiQ6OnOykvaOL9o4uDp3o4ERHF+2nujjRkS3tHd0VaNW5lBKE0rp63mskzonpdExn3n5UHNOZtyKp6Jycc3zpc3K6fK+6DKJN533MoL5pcMl00Ol3kAcO5rDB/pGQ5z8tnviPv01d9dDex+RkMMrV1xSYO3Ucc6eOG/Q5uruDU13ddHR109kVdHR109EddHR209ndzanOLNaz3tmdle3oCiKCrm7oip71oDuyc3ZH0BW9truDCFI8UpwUj3SeM6+yy9azwJnYmTKn96f10/t7Hd/7nPTEeo6h9DmLv/d8DaZXPfjvGsQxg/6uwR05qKMGWckYdOvGhsH+8dIfJ4McqKoS9VUF3xFtZn3yTWdmZuZkYGZmoygZSFom6VVJrZLuqHR9zMzyZFQkA0kF4H7gOmAhcLOkhZWtlZlZfoyKZAAsBlojYkdEnALWAMsrXCczs9wYLclgNvB60fauFDuLpJWSWiS1tLW1jVjlzMzGutGSDEpNmj1nInFEPBARzRHR3NTUNALVMjPLh9GSDHYBc4u25wC7K1QXM7PcGRXPJpJUDfwCuAZ4A9gE/EFEbO3nmDbgV4P8yunAm4M89mLlNudD3tqct/bChbf57RFxztDKqLgDOSI6Jd0OrAcKwKr+EkE6ZtDjRJJaSj2oaSxzm/Mhb23OW3th+No8KpIBQESsA9ZVuh5mZnk0Wq4ZmJlZBeU1GTxQ6QpUgNucD3lrc97aC8PU5lFxAdnMzCorrz0DMzMr4mRgZmb5SgZj6cmoklZJ2ifp5aLYVEkbJG1Pn1NSXJLuS+3eIumKomNWpPLbJa2oRFvKJWmupGckbZO0VdJnUnzMtltSvaTnJL2Y2vzFFJ8vaWOq/yOSalO8Lm23pv3zis51Z4q/KunayrSoPJIKkl6Q9ETaHtPtBZC0U9JLkjZLakmxkfvZjvQ6w7G+kN2/8EvgUqAWeBFYWOl6XUB7fge4Ani5KPZfgTvS+h3APWn9euBJssd+LAE2pvhUYEf6nJLWp1S6bf20eRZwRVqfSHaj4sKx3O5U9wlpvQbYmNryKHBTiv8D8Km0/mngH9L6TcAjaX1h+pmvA+an/wuFSrevn3Z/Dvgn4Im0Pabbm+q8E5jeKzZiP9t56hmMqSejRsSPgf29wsuB1Wl9NXBDUfzhyDwLNEqaBVwLbIiI/RFxANgALBv+2g9OROyJiJ+l9SPANrIHGo7Zdqe6H02bNWkJ4GrgsRTv3eaef4vHgGuUvXF+ObAmIk5GxGtAK9n/iVFH0hzgY8A/pm0xhts7gBH72c5TMijryagXuZkRsQeyX5zAjBTvq+0X7b9JGg74ANlfymO63WnIZDOwj+w/9y+BgxHRmYoU1/9029L+Q8A0Lq42fx34M6A7bU9jbLe3RwA/lPS8pJUpNmI/26PmDuQRUNaTUceovtp+Uf6bSJoAfA/4bEQczv4QLF20ROyia3dEdAGLJDUCPwAuL1UsfV7UbZb0e8C+iHhe0od6wiWKjon29nJVROyWNAPYIOnn/ZQd8nbnqWeQhyej7k1dRdLnvhTvq+0X3b+JpBqyRPDtiPh+Co/5dgNExEHgR2RjxI3KHvAIZ9f/dNvS/slkw4kXS5uvAn5f0k6yodyryXoKY7W9p0XE7vS5jyzpL2YEf7bzlAw2AQvSrIRasotNaytcp6G2FuiZPbACeLwofkuagbAEOJS6nOuBpZKmpFkKS1NsVEpjwQ8C2yLia0W7xmy7JTWlHgGSGoCPkF0reQa4MRXr3eaef4sbgacju7K4Frgpzb6ZDywAnhuZVpQvIu6MiDkRMY/s/+jTEfGHjNH29pA0XtLEnnWyn8mXGcmf7UpfQR/JhewK/C/Ixlz/otL1ucC2fAfYA3SQ/TVwK9lY6VPA9vQ5NZUV2Tumfwm8BDQXnec/kF1cawU+Wel2DdDm3ybr8m4BNqfl+rHcbuB9wAupzS8D/yXFLyX75dYKfBeoS/H6tN2a9l9adK6/SP8WrwLXVbptZbT9Q5yZTTSm25va92Jatvb8fhrJn20/jsLMzHI1TGRmZn1wMjAzMycDMzNzMjAzM5wMzMwMJwMzM8PJwMzMgP8PExPWWdEi2UoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(history[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.04864465],\n",
       "       [-0.08872099],\n",
       "       [ 0.02445404],\n",
       "       [ 0.01577132],\n",
       "       [ 0.02333752],\n",
       "       [-0.08034107],\n",
       "       [ 0.08589046],\n",
       "       [ 0.01040775],\n",
       "       [-0.10300676],\n",
       "       [ 0.13523935],\n",
       "       [-0.09143892],\n",
       "       [-0.07765697],\n",
       "       [ 0.04368575],\n",
       "       [-0.22417392]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta # 위에서 사용한 x와 약간 다름 : standardized 여부"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1416.0202939098758"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(beta, XX, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_yhat = np.exp(np.dot(X, W) )\n",
    "my_yhat2 = np.exp(np.dot(XX, beta) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([24.69482883, 31.57108125, 29.83295113, 28.29189   , 26.02959171,\n",
       "       22.14341111, 18.13403117, 11.96775387, 18.19758204])"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_yhat[1:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[24.69994292],\n",
       "       [31.5711263 ],\n",
       "       [29.83283516],\n",
       "       [28.29114655],\n",
       "       [26.03383602],\n",
       "       [22.15077532],\n",
       "       [18.14084703],\n",
       "       [11.97271604],\n",
       "       [18.20521123]])"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_yhat2[1:10]"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}