neshi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.linear_model import LinearRegression, Ridge, ElasticNet\n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics.pairwise import rbf_kernel, polynomial_kernel\n",
    "import statsmodels.api as sm\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('shojisu_kaiten2.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = df.iloc[:, 0:8].values\n",
    "y = df.iloc[:, 8].values.reshape(-1, 1)\n",
    "x = np.delete(x, 22, axis=0)\n",
    "y = np.delete(y, 22, axis=0)\n",
    "logx = np.hstack((x[:, 0:1], np.log(x[:,1:]+1)))\n",
    "logy = np.log(y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = x[x[:, 0] == 0]\n",
    "y0 = y[x[:, 0] == 0]\n",
    "x1 = x[x[:, 0] == 1]\n",
    "y1 = y[x[:, 0] == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GLMWrapper():\n",
    "    def __init__(self, family, fit_intercept=True, **glm_params):\n",
    "        self.family = family\n",
    "        self.fit_intercept = fit_intercept \n",
    "        self.glm_params = glm_params\n",
    "        \n",
    "    def fit(self, x, y, sample_weights=None):\n",
    "        _x = sm.add_constant(x) if self.fit_intercept else x\n",
    "        self.glm = sm.GLM(y, _x, family=self.family, \n",
    "                          freq_weights=sample_weights, **self.glm_params)\n",
    "        self.results = self.glm.fit()\n",
    "        return self\n",
    "\n",
    "    def score(self, x, y, sample_weights=None):\n",
    "        self.fit(x, y, sample_weights)\n",
    "        return self.glm.score(self.results.params)\n",
    "    \n",
    "    def predict(self, x):\n",
    "        _x = sm.add_constant(x) if self.fit_intercept else x\n",
    "        return self.glm.predict(self.results.params, _x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold\n",
    "kfold = KFold(5, shuffle=True)\n",
    "split = kfold.split(x)\n",
    "split = list(split)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(42, 8)"
      ]
     },
     "execution_count": 281,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_dataset(i, _x, _y):\n",
    "    return _x[split[i][0]], _y[split[i][0]], _x[split[i][1]], _y[split[i][1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.6/site-packages/scipy/linalg/basic.py:223: RuntimeWarning: scipy.linalg.solve\n",
      "Ill-conditioned matrix detected. Result is not guaranteed to be accurate.\n",
      "Reciprocal condition number: 9.221851574618809e-20\n",
      "  ' condition number: {}'.format(rcond), RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "KernelRidge(alpha=0, coef0=1, degree=3, gamma=1e-07, kernel='linear',\n",
       "      kernel_params=None)"
      ]
     },
     "execution_count": 275,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.kernel_ridge import KernelRidge\n",
    "kr = KernelRidge(alpha=0, kernel='linear', gamma=0.0000001)\n",
    "kr.fit(train_x, train_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_y = kr.predict(test_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "147994343.75911847"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean((test_y-pred_y)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6018182.34437\n",
      "12386836.578\n",
      "30947399.7867\n",
      "30877312.7392\n",
      "1318014.75806\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "13950382.090074688"
      ]
     },
     "execution_count": 227,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred_ys = []\n",
    "for i in range(5):\n",
    "    train_x, train_y, test_x, test_y = get_dataset(i, logx, logy)\n",
    "    k_train_x = rbf_kernel(train_x, gamma=0.22)\n",
    "    k_test_x = rbf_kernel(test_x, train_x, gamma=0.22)\n",
    "    kx = rbf_kernel(logx, train_x, gamma=0.22)\n",
    "    #kr = KernelRidge(alpha=0.001, kernel='rbf', degree=5, gamma=0.22)\n",
    "    #kr = GLMWrapper(sm.families.Gaussian())\n",
    "    kr = LinearRegression(normalize=True) \n",
    "    kr.fit(train_x, train_y)\n",
    "    pred_y = kr.predict(test_x)\n",
    "    pred_y[pred_y < 0] = 0\n",
    "    print(np.mean((np.exp(test_y)-np.exp(pred_y))**2))\n",
    "    pred_y = kr.predict(logx)\n",
    "    pred_y[pred_y < 0] = 0\n",
    "    pred_ys.append(pred_y)\n",
    "np.mean((np.exp(np.mean(pred_ys, axis=0)) - y)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.78022085517680417"
      ]
     },
     "execution_count": 229,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kr.score(train_x, train_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   1.,   0.,\n",
       "          0.,   1.,   0.,   0.,   0.,   0.,   1.,   0.,   0.,   0.,   0.,\n",
       "          1.,   0.,   0.,   0.,   0.,   1.,   0.,   2.,   1.,   1.,   1.,\n",
       "          4.,   5.,   5.,   5.,   9.,  17.,  27.]),\n",
       " array([-38155.11364212, -37203.23290137, -36251.35216062, -35299.47141987,\n",
       "        -34347.59067912, -33395.70993837, -32443.82919763, -31491.94845688,\n",
       "        -30540.06771613, -29588.18697538, -28636.30623463, -27684.42549388,\n",
       "        -26732.54475313, -25780.66401238, -24828.78327163, -23876.90253088,\n",
       "        -22925.02179013, -21973.14104938, -21021.26030863, -20069.37956788,\n",
       "        -19117.49882713, -18165.61808638, -17213.73734563, -16261.85660488,\n",
       "        -15309.97586413, -14358.09512338, -13406.21438263, -12454.33364188,\n",
       "        -11502.45290113, -10550.57216038,  -9598.69141963,  -8646.81067888,\n",
       "         -7694.92993813,  -6743.04919738,  -5791.16845663,  -4839.28771588,\n",
       "         -3887.40697513,  -2935.52623438,  -1983.64549363,  -1031.76475288,\n",
       "           -79.88401213]),\n",
       " <a list of 40 Patch objects>)"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAADt1JREFUeJzt3X+MHOV9x/H3tzhQqUkbE5+cK4EeRDSqaRWSnNzQtBUVaWpMK4L6Q1gVcptUjhrcJhX9w8ZSQKoqkaQkUtSUYAQKqmgSGkKgwSlxUGgUqSU9qDEG4mKoKbaMfSRqoFKVyvDtH/scXl9vb3/M7N3x+P2SVjf7zMzO52bvPreemV1HZiJJeu37seUOIElqh4UuSZWw0CWpEha6JFXCQpekSljoklQJC12SKmGhS1IlLHRJqsSqpdzYmjVrcmpqaik3KUmveQ8//PALmTnRb7klLfSpqSlmZmaWcpOS9JoXEc8OspyHXCSpEha6JFXCQpekSljoklQJC12SKmGhS1IlLHRJqoSFLkmVsNAlqRJL+k5RSarV1Lb7Fp1/8IbLxp7BV+iSVAkLXZIqYaFLUiUsdEmqhIUuSZWw0CWpEha6JFXCQpekSljoklQJC12SKtG30CPi7Ij4VkQ8ERGPR8RHyvj1EXE4IvaU28bxx5Uk9TLIZ7kcB67JzEci4g3AwxGxu8z7dGb+1fjiSZIG1bfQM/MIcKRMvxQRTwJnjTuYJGk4Qx1Dj4gp4B3AQ2Voa0TsjYjbImJ1y9kkSUMYuNAj4vXAXcBHM/NF4CbgrcCFdF7B39hjvS0RMRMRM7Ozsy1EliQtZKBCj4jX0SnzOzLzKwCZeTQzX87MV4BbgPULrZuZOzNzOjOnJyYm2sotSZpnkKtcArgVeDIzP9U1Ptm12BXAvvbjSZIGNchVLu8BrgIei4g9ZexaYFNEXAgkcBD40FgSSpIGMshVLt8BYoFZu9qPI0kale8UlaRKWOiSVAkLXZIqYaFLUiUsdEmqhIUuSZWw0CWpEha6JFXCQpekSljoklQJC12SKmGhS1IlLHRJqoSFLkmVsNAlqRIWuiRVwkKXpEpY6JJUCQtdkiphoUtSJSx0SaqEhS5JlbDQJakSFrokVcJCl6RKWOiSVAkLXZIqYaFLUiUsdEmqhIUuSZXoW+gRcXZEfCsinoiIxyPiI2X8zIjYHRFPla+rxx9XktTLIK/QjwPXZOY64N3A1RGxDtgGPJCZ5wMPlPuSpGXSt9Az80hmPlKmXwKeBM4CLgduL4vdDrx/XCElSf0NdQw9IqaAdwAPAWsz80iZ9Tywtsc6WyJiJiJmZmdnG0SVJC1m4EKPiNcDdwEfzcwXu+dlZgK50HqZuTMzpzNzemJiolFYSVJvAxV6RLyOTpnfkZlfKcNHI2KyzJ8Ejo0noiRpEINc5RLArcCTmfmprln3ApvL9GbgnvbjSZIGtWqAZd4DXAU8FhF7yti1wA3AnRHxQeBZ4PfGE1GSNIi+hZ6Z3wGix+xL2o0jSRqV7xSVpEpY6JJUCQtdkiphoUtSJSx0SaqEhS5JlbDQJakSFrokVcJCl6RKWOiSVAkLXZIqYaFLUiUsdEmqhIUuSZWw0CWpEha6JFXCQpekSgzyX9BJ0ilvatt9yx2hL1+hS1IlLHRJqoSFLkmVsNAlqRIWuiRVwkKXpEpY6JJUCQtdkiphoUtSJSx0SaqEhS5Jlehb6BFxW0Qci4h9XWPXR8ThiNhTbhvHG1OS1M8gr9A/D2xYYPzTmXlhue1qN5YkaVh9Cz0zvw38YAmySJIaaHIMfWtE7C2HZFa3lkiSNJJRC/0m4K3AhcAR4MZeC0bEloiYiYiZ2dnZETcnSepnpELPzKOZ+XJmvgLcAqxfZNmdmTmdmdMTExOj5pQk9TFSoUfEZNfdK4B9vZaVJC2Nvv8FXUR8AbgYWBMRh4DrgIsj4kIggYPAh8aYUZI0gL6FnpmbFhi+dQxZJEkN+E5RSaqEhS5JlbDQJakSFrokVcJCl6RKWOiSVAkLXZIqYaFLUiUsdEmqhIUuSZWw0CWpEha6JFXCQpekSljoklQJC12SKmGhS1IlLHRJqoSFLkmVsNAlqRIWuiRVwkKXpEpY6JJUCQtdkiphoUtSJSx0SaqEhS5JlbDQJakSFrokVcJCl6RK9C30iLgtIo5FxL6usTMjYndEPFW+rh5vTElSP4O8Qv88sGHe2Dbggcw8H3ig3JckLaO+hZ6Z3wZ+MG/4cuD2Mn078P6Wc0mShjTqMfS1mXmkTD8PrG0pjyRpRI1PimZmAtlrfkRsiYiZiJiZnZ1tujlJUg+jFvrRiJgEKF+P9VowM3dm5nRmTk9MTIy4OUlSP6MW+r3A5jK9GbinnTiSpFENctniF4B/Bt4WEYci4oPADcCvR8RTwHvLfUnSMlrVb4HM3NRj1iUtZ5EkNeA7RSWpEha6JFXCQpekSvQ9hi5Jp4qpbfctd4RGfIUuSZWw0CWpEha6JFXCQpekSljoklQJC12SKmGhS1IlLHRJqoSFLkmVsNAlqRIWuiRVwkKXpEpY6JJUCQtdkiphoUtSJSx0SaqEhS5JlbDQJakSFrokVcJCl6RKWOiSVAkLXZIqYaFLUiUsdEmqhIUuSZWw0CWpEquarBwRB4GXgJeB45k53UYoSdLwGhV68WuZ+UILjyNJasBDLpJUiaav0BP4RkQkcHNm7py/QERsAbYAnHPOOQ03J+lUN7XtvuWOsGI1fYX+y5n5TuBS4OqI+NX5C2TmzsyczszpiYmJhpuTJPXSqNAz83D5egy4G1jfRihJ0vBGLvSI+ImIeMPcNPA+YF9bwSRJw2lyDH0tcHdEzD3O32XmP7aSSpI0tJELPTOfAd7eYhZJUgNetihJlbDQJakSbbxTVJJa43Xmo/MVuiRVwkKXpEpY6JJUCQtdkiphoUtSJSx0SaqEhS5JlbDQJakSFrokVcJCl6RKWOiSVAkLXZIqYaFLUiUsdEmqhB+fK2lBi32M7cEbLlvCJBqUr9AlqRIWuiRVwkKXpEpY6JJUCQtdkiphoUtSJSx0SarEa+Y69MWuiYVT87rYcV0nvJL39Uq9NvpUy9XvZ0TLw1foklQJC12SKmGhS1IlGhV6RGyIiP0RcSAitrUVSpI0vJELPSJOAz4LXAqsAzZFxLq2gkmShtPkFfp64EBmPpOZ/wt8Ebi8nViSpGE1KfSzgOe67h8qY5KkZRCZOdqKEb8DbMjMPyr3rwJ+MTO3zltuC7Cl3P15YN/occduDfDCcodYhPmaMV8z5mumSb6fycyJfgs1eWPRYeDsrvtvKWMnycydwE6AiJjJzOkG2xwr8zVjvmbM14z5mh1y+Vfg/Ig4NyJOB64E7m0nliRpWCO/Qs/M4xGxFbgfOA24LTMfby2ZJGkojT7LJTN3AbuGWGVnk+0tAfM1Y75mzNfMKZ9v5JOikqSVxbf+S1ItMrOVG3ANkMCacj+AzwAHgL3AO7uW3Qw8VW6bu8bfBTxW1vkMJ/4FcSawuyy/G1g9RK6/KNvfA3wD+OkyfjHwwzK+B/hY1zobgP0lx7au8XOBh8r4l4DTy/gZ5f6BMn+qhXzLvv+ATwLfK9u/G3hjGZ8C/qdr331u1AyLfZ9NMpZ528vj7gd+Y6mfX+B3gceBV4DprvEVsf965VsJ+65H3uvpXEk3t982tp13XLdeOVrfTkthz6ZzcvRZThT6RuDr5Qfu3cBDXT+Yz5Svq8v03A/nd8uyUda9tIx/Ym4nANuAjw+R7Se7pv907peHTqF/bYHlTwOeBs4DTgceBdaVeXcCV5bpzwF/XKY/3PW4VwJfaiHfsu8/4H3AqjL98bn16BTSvh7rDJWh1/c5xP7rlXFdee7OoPPL+3R5bpfs+QV+Dngb8CD/v9CXff8tkm/Z912PvNcDf77AeGt5x3FbLEfr22op8JeBtwMHOVHoNwObupbZD0wCm4Cbu8ZvLmOTwPe6xl9dbm7dMj0J7B8x53bgpjJ9MQsX+kXA/fPW2V5+YV7gRHm8uhydP2YXlelVZblomG9F7T/gCuCOMj3FAoU0SoZe3+eIz293xu3A9q5595fnbMmfXwYs9OXafwvkWzH7bl7O61m40FvLO45brxzj2FbjY+gRcTlwODMfnTer10cDLDZ+aIFxgLWZeaRMPw+sHTLjX0bEc8DvAx/rmnVRRDwaEV+PiAv65H4T8F+ZeXyBfK+uU+b/sCzfJN+K2X/FB+i8EpxzbkT8W0T8U0T8SlfmYTO0+RES3RmH3X9je357WIn7b85K3ndbI2JvRNwWEavHkHccluxjUga6bDEivgm8eYFZO4Br6fyzd0lkZkZEdo8tli8z78nMHcCOiNgObAWuAx6h83ba/46IjcBXgfPHkXnEfGMxf//1y1aW2QEcB+4o844A52Tm9yPiXcBXu/4gDp2hnxEzLolBsi1gyfbfiPmWTZ+uuYnOOacsX2+k80dcxUCFnpnvXWg8In6BzjGrRyMCOm//fyQi1tP7owEO0znc0T3+YBl/ywLLAxyNiMnMPBIRk8CxQfIt4A46181fl5kvdq2/KyL+JiLWLJL7+8AbI2JV+cvenW9unUMRsQr4qbL8yPkWydHq/uuXLSL+APhN4JIs/17MzB8BPyrTD0fE08DPjpih70dIjJKxz+O29vwO8dx2fz9Ltv9GydfnMVv93Zhv0LwRcQvwtTHkHYeBPialFS0fKzrIiWPol3HyyZrvlvEzgf+gc0JvdZk+s8ybf0JoYxn/JCefEPrEEJnO75r+E+DLZfrNnLiCYD3wn2W7q+icaDyXEycwLijL/T0nn0j5cJm+mpNP/NzZQr5l3390zsw/AUzMG58ATivT59H54RwpQ6/vc4j91yvjBZx8ouwZOienlvT5Les8yMnHqFfM/uuRb8Xsu3k5J7um/wz4Ytt5x3FbLEfr22o5+EFOvmzxs3TO7j427wfmA3Qu3zkA/GHX+DSdT2N8GvhrThTum4AH6Fyy9c25H/4BM91VHnMv8A/AWWV8K51Lth4F/gX4pa51NgL/XnLs6Bo/r/zCHSg/EGeU8R8v9w+U+ee1kG/Z9195/OeYd3kd8Ntl3+2hc+jqt0bNsNj32SRjmbejPO5+ytUiS/n80jlJe4jOq/GjnDhRuCL2X698K2Hf9cj7t+V73Evnc6Mm2847rluvHG3ffKeoJFXCd4pKUiUsdEmqhIUuSZWw0CWpEha6JFXCQpekSljoklQJC12SKvF/gZ+AMMyMuMIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11fc88828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist((np.mean(pred_ys, axis=0)-y).flatten(), bins=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1.,  5.,  2.,  2.,  5.,  6.,  5.,  6.,  8.,  2.,  2.,  2.,  4.,\n",
       "         4.,  3.,  4.,  2.,  1.,  2.,  3.,  2.,  1.,  0.,  2.,  2.,  2.,\n",
       "         0.,  1.,  0.,  1.,  0.,  1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]),\n",
       " array([ 0.14917904,  0.23525483,  0.32133061,  0.4074064 ,  0.49348219,\n",
       "         0.57955798,  0.66563377,  0.75170956,  0.83778534,  0.92386113,\n",
       "         1.00993692,  1.09601271,  1.1820885 ,  1.26816429,  1.35424007,\n",
       "         1.44031586,  1.52639165,  1.61246744,  1.69854323,  1.78461902,\n",
       "         1.87069481,  1.95677059,  2.04284638,  2.12892217,  2.21499796,\n",
       "         2.30107375,  2.38714954,  2.47322532,  2.55930111,  2.6453769 ,\n",
       "         2.73145269,  2.81752848,  2.90360427,  2.98968005,  3.07575584,\n",
       "         3.16183163,  3.24790742,  3.33398321,  3.420059  ,  3.50613479,\n",
       "         3.59221057]),\n",
       " <a list of 40 Patch objects>)"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAADc5JREFUeJzt3W2spOVdx/Hvr7uLQCHFdCdKgOPBxJAgkYeeIBVDKogB1sALebEkRWlsTqLWUmNiVl9IamKyL0zjY6wnLYoW6QMFg2ypJSm1aWK37q7QLiwYiksLoru04ama4pK/L84sbA8zZ+7dc+4zc5HvJ5kwc+5rZn7cGX7c55r7OneqCklSO9427QCSpONjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5Ias7mPF926dWvNz8/38dKS9Ja0d+/e56tq0GVsL8U9Pz/Pnj17+nhpSXpLSvJ017FOlUhSYyxuSWqMxS1JjbG4JakxFrckNaZTcSf5rSSPJtmf5K4kJ/cdTJI02sTiTnIW8EFgoaouADYB2/sOJkkaretUyWbglCSbgVOB/+wvkiRpNROLu6qeBf4I+BbwHPBiVX2h72CSpNEmrpxM8sPADcC5wAvAZ5K8t6o+sWLcIrAIMDc310PUt675HbvGbju4c9sGJpHUgi5TJT8P/EdVHa6q/wPuAX5m5aCqWqqqhapaGAw6LbeXJJ2ALsX9LeCyJKcmCXAVcKDfWJKkcbrMce8G7gb2Ad8YPmep51ySpDE6/XXAqroNuK3nLJKkDlw5KUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2ZWNxJzkvy8DG3l5J8aCPCSZLebOKly6rqCeAigCSbgGeBe3vOJUka43inSq4CvllVT/cRRpI02fEW93bgrj6CSJK66VzcSU4Crgc+M2b7YpI9SfYcPnx4vfJJklY4niPua4F9VfXfozZW1VJVLVTVwmAwWJ90kqQ3OZ7ivgmnSSRp6joVd5K3A1cD9/QbR5I0ycTTAQGq6nvAO3vOIknqwJWTktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1Jiuly47I8ndSR5PciDJu/sOJkkardOly4A/AT5fVTcmOQk4tcdMkqRVTCzuJO8ArgBuAaiqV4FX+40lSRqnyxH3ucBh4K+TXAjsBW4dXkD4dUkWgUWAubm59c45dfM7do3ddnDnthN+riQdry5z3JuBS4C/rKqLge8BO1YOqqqlqlqoqoXBYLDOMSVJR3Up7meAZ6pq9/Dx3SwXuSRpCiYWd1X9F/DtJOcNf3QV8FivqSRJY3U9q+Q3gTuHZ5Q8Bbyvv0iSpNV0Ku6qehhY6DmLJKkDV05KUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYzpdASfJQeBl4DXgSFV5NRxJmpKu15wE+Lmqer63JJKkTpwqkaTGdD3iLuALSQr4q6paWjkgySKwCDA3N7d+CY8xv2PXqtsP7tzW22v39VxJOl5dj7h/tqouAa4FfiPJFSsHVNVSVS1U1cJgMFjXkJKkN3Qq7qp6dvjPQ8C9wKV9hpIkjTexuJO8PcnpR+8DvwDs7zuYJGm0LnPcPwLcm+To+L+vqs/3mkqSNNbE4q6qp4ALNyCLJKkDTweUpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxnQu7iSbkvxbkvv7DCRJWt3xHHHfChzoK4gkqZtOxZ3kbGAb8LF+40iSJulylXeAPwZ+Bzh93IAki8AiwNzc3NqTqXfzO3ad8HMP7ty2ptee9HxJ40084k7yi8Chqtq72riqWqqqhapaGAwG6xZQkvSDukyVXA5cn+Qg8EngyiSf6DWVJGmsicVdVb9bVWdX1TywHfhiVb2392SSpJE8j1uSGtP1y0kAqupLwJd6SSJJ6sQjbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSY7pcLPjkJF9L8kiSR5N8eCOCSZJG63IFnO8DV1bVK0m2AF9J8kBVfbXnbJKkESYWd1UV8Mrw4ZbhrfoMJUkar9Mcd5JNSR4GDgEPVtXufmNJksbpdLHgqnoNuCjJGcC9SS6oqv3HjkmyCCwCzM3NrXtQzZb5HbumHWGkSbkO7ty2QUmk/hzXWSVV9QLwEHDNiG1LVbVQVQuDwWC98kmSVuhyVslgeKRNklOAq4HH+w4mSRqty1TJmcAdSTaxXPSfrqr7+40lSRqny1klXwcu3oAskqQOXDkpSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjelyzclzkjyU5LEkjya5dSOCSZJG63LNySPAb1fVviSnA3uTPFhVj/WcTZI0wsQj7qp6rqr2De+/DBwAzuo7mCRptOOa404yz/KFg3f3EUaSNFmqqtvA5DTgn4E/rKp7RmxfBBYB5ubm3vX000+vZ04A5nfsWvfX1Ow5uHPbqtv7/BxMeu/VzGoutSHJ3qpa6DK20xF3ki3AZ4E7R5U2QFUtVdVCVS0MBoPuaSVJx6XLWSUBPg4cqKqP9B9JkrSaLkfclwM3A1cmeXh4u67nXJKkMSaeDlhVXwGyAVkkSR24clKSGmNxS1JjLG5JaozFLUmNsbglqTEWtyQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4JakxFrckNcbilqTGWNyS1BiLW5Ia0+Wak7cnOZRk/0YEkiStrssR998A1/ScQ5LU0cTirqovA9/dgCySpA6c45akxky8yntXSRaBRYC5ubkTfp35HbvWK5IaNc3PwGrvfXDntg1M0t2k/bVa7rU8t8vz1/Las2oWPiPrdsRdVUtVtVBVC4PBYL1eVpK0glMlktSYLqcD3gX8C3BekmeS/Gr/sSRJ40yc466qmzYiiCSpG6dKJKkxFrckNcbilqTGWNyS1BiLW5IaY3FLUmMsbklqjMUtSY2xuCWpMRa3JDXG4pakxljcktQYi1uSGmNxS1JjLG5JaozFLUmNsbglqTGdijvJNUmeSPJkkh19h5IkjdflmpObgL8ArgXOB25Kcn7fwSRJo3U54r4UeLKqnqqqV4FPAjf0G0uSNE6X4j4L+PYxj58Z/kySNAWpqtUHJDcC11TV+4ePbwZ+uqo+sGLcIrA4fHge8ASwFXh+vUP3yLz9Mm//Wsts3jf8WFUNugzc3GHMs8A5xzw+e/izH1BVS8DSsT9LsqeqFroEmQXm7Zd5+9daZvOemC5TJf8K/ESSc5OcBGwH7us3liRpnIlH3FV1JMkHgH8CNgG3V9WjvSeTJI3UZaqEqvoc8LkTeP2lyUNminn7Zd7+tZbZvCdg4peTkqTZ4pJ3SWrMmot70nL4JD+U5FPD7buTzK/1PdeiQ95bkhxO8vDw9v5p5Dwmz+1JDiXZP2Z7kvzp8N/n60ku2eiMK/JMyvueJC8es39/f6MzrshzTpKHkjyW5NEkt44YMzP7uGPeWdvHJyf5WpJHhpk/PGLMzPREx7zT7YmqOuEby19WfhP4ceAk4BHg/BVjfh346PD+duBTa3nPDch7C/Dn08o4IvMVwCXA/jHbrwMeAAJcBuye8bzvAe6f9n49Js+ZwCXD+6cD/z7iMzEz+7hj3lnbxwFOG97fAuwGLlsxZpZ6okveqfbEWo+4uyyHvwG4Y3j/buCqJFnj+56o5pbvV9WXge+uMuQG4G9r2VeBM5KcuTHp3qxD3plSVc9V1b7h/ZeBA7x5ZfDM7OOOeWfKcL+9Mny4ZXhb+eXazPREx7xTtdbi7rIc/vUxVXUEeBF45xrf90R1Xb7/S8Nfie9Ocs6I7bOkxT9J8O7hr6EPJPnJaYc5avjr+cUsH2Edayb38Sp5Ycb2cZJNSR4GDgEPVtXYfTwDPdElL0yxJ/xy8s3+EZivqp8CHuSNowCtj30sL+29EPgz4B+mnAeAJKcBnwU+VFUvTTvPJBPyztw+rqrXquoilldeX5rkgmlnWk2HvFPtibUWd5fl8K+PSbIZeAfwnTW+74mamLeqvlNV3x8+/Bjwrg3KdqI6/UmCWVFVLx39NbSW1wdsSbJ1mpmSbGG5BO+sqntGDJmpfTwp7yzu46Oq6gXgIeCaFZtmqSdeNy7vtHtircXdZTn8fcCvDO/fCHyxhrP7UzAx74q5y+tZnkOcZfcBvzw88+Ey4MWqem7aocZJ8qNH5y6TXMryZ3Bq/4EOs3wcOFBVHxkzbGb2cZe8M7iPB0nOGN4/BbgaeHzFsJnpiS55p90TnVZOjlNjlsMn+QNgT1Xdx/KH7O+SPMnyl1bb1xq657wfTHI9cGSY95Zp5QVIchfLZwlsTfIMcBvLX5ZQVR9leUXrdcCTwP8A75tO0mUd8t4I/FqSI8D/Atun+D9ygMuBm4FvDOc0AX4PmIOZ3Mdd8s7aPj4TuCPLF2V5G/Dpqrp/VnuCbnmn2hOunJSkxvjlpCQ1xuKWpMZY3JLUGItbkhpjcUtSYyxuSWqMxS1JjbG4Jakx/w/3qavry7LBgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11fc19ac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist((y/(np.mean(np.exp(pred_ys), axis=0))).flatten(), bins=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (83,1) (42,1) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-127-ed027b09e54a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpred_ys\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0my0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (83,1) (42,1) "
     ]
    }
   ],
   "source": [
    "(((np.mean(pred_ys, axis=0))-y0).flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.delete(y0, 39)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.delete(x0, 39, axis=0)\n",
    "y0 = np.delete(y0, 39, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(42, 8)"
      ]
     },
     "execution_count": 287,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[statsmodels.genmod.families.links.log,\n",
       " statsmodels.genmod.families.links.cloglog,\n",
       " statsmodels.genmod.families.links.identity,\n",
       " statsmodels.genmod.families.links.nbinom,\n",
       " statsmodels.genmod.families.links.Power]"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sm.families.family.NegativeBinomial.links"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}