Skip to content

Instantly share code, notes, and snippets.

@ssugiyama
Created March 27, 2018 01:06
Show Gist options
  • Save ssugiyama/75bc2c8792a983560169890c7873fe0c to your computer and use it in GitHub Desktop.
Save ssugiyama/75bc2c8792a983560169890c7873fe0c to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment