Empirical bayes - Chapter 8.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import scipy\n",
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "from patsy import dmatrix\n",
    "sns.set(style=\"darkgrid\")\n",
    "pd.set_option('display.max_rows', 500)\n",
    "pd.set_option('display.max_columns', 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create batting average dataset (sum over years)\n",
    "batting = pd.read_csv(\"~/baseball/core/Batting.csv\")\n",
    "pitching = pd.read_csv(\"~/baseball/core/Pitching.csv\")\n",
    "# create batting average dataset (sum over years)\n",
    "batting = pd.read_csv(\"~/baseball/core/Batting.csv\")\n",
    "pitching = pd.read_csv(\"~/baseball/core/Pitching.csv\")\n",
    "pitching = pitching.groupby(\"playerID\").agg({\"G\": \"sum\"}).reset_index()\n",
    "pitching = pitching.loc[pitching[\"G\"]>3]\n",
    "# filter pitchers with 4+ games\n",
    "batting_pitchers = batting.playerID.isin(pitching.playerID)\n",
    "batting = batting[~batting_pitchers]\n",
    "# sum over seasons\n",
    "batting = batting.groupby(\"playerID\").agg({\"AB\": \"sum\", \"H\": \"sum\", \"yearID\":\"mean\"}).reset_index()\n",
    "batting = batting.loc[batting[\"AB\"] > 0]\n",
    "batting[\"average\"] = batting[\"H\"]/batting[\"AB\"]\n",
    "# add actual player name\n",
    "people = pd.read_csv(\"~/baseball/core/People.csv\")\n",
    "people[\"name\"] = people[\"nameFirst\"] + \" \" + people[\"nameLast\"]\n",
    "batting = batting.merge(people[[\"playerID\", \"name\", \"bats\"]], on=\"playerID\")\n",
    "# I am using the book values here, not the values that we got in our last notebook. I think the difference comes\n",
    "# from the fact that we are using more data (data from more recent seasons).\n",
    "alpha_0 = 101.4\n",
    "beta_0 = 287.3\n",
    "# add empirical Bayes estimate\n",
    "# this is called a point estimate\n",
    "batting[\"eb_estimate\"] = (batting[\"H\"] + alpha_0) / (batting[\"AB\"] + alpha_0 + beta_0)\n",
    "# sort by eb_estimate\n",
    "batting = batting.sort_values(\"eb_estimate\", ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>playerID</th>\n",
       "      <th>AB</th>\n",
       "      <th>H</th>\n",
       "      <th>yearID</th>\n",
       "      <th>average</th>\n",
       "      <th>name</th>\n",
       "      <th>bats</th>\n",
       "      <th>eb_estimate</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1791</th>\n",
       "      <td>cobbty01</td>\n",
       "      <td>11435</td>\n",
       "      <td>4189</td>\n",
       "      <td>1916.500000</td>\n",
       "      <td>0.366331</td>\n",
       "      <td>Ty Cobb</td>\n",
       "      <td>L</td>\n",
       "      <td>0.362864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4512</th>\n",
       "      <td>hornsro01</td>\n",
       "      <td>8173</td>\n",
       "      <td>2930</td>\n",
       "      <td>1926.291667</td>\n",
       "      <td>0.358497</td>\n",
       "      <td>Rogers Hornsby</td>\n",
       "      <td>R</td>\n",
       "      <td>0.354065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4720</th>\n",
       "      <td>jacksjo01</td>\n",
       "      <td>4981</td>\n",
       "      <td>1772</td>\n",
       "      <td>1914.071429</td>\n",
       "      <td>0.355752</td>\n",
       "      <td>Shoeless Joe Jackson</td>\n",
       "      <td>L</td>\n",
       "      <td>0.348884</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>465</th>\n",
       "      <td>barnero01</td>\n",
       "      <td>2391</td>\n",
       "      <td>860</td>\n",
       "      <td>1875.333333</td>\n",
       "      <td>0.359682</td>\n",
       "      <td>Ross Barnes</td>\n",
       "      <td>R</td>\n",
       "      <td>0.345865</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2377</th>\n",
       "      <td>delahed01</td>\n",
       "      <td>7510</td>\n",
       "      <td>2597</td>\n",
       "      <td>1895.500000</td>\n",
       "      <td>0.345806</td>\n",
       "      <td>Ed Delahanty</td>\n",
       "      <td>R</td>\n",
       "      <td>0.341626</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       playerID     AB     H       yearID   average                  name  \\\n",
       "1791   cobbty01  11435  4189  1916.500000  0.366331               Ty Cobb   \n",
       "4512  hornsro01   8173  2930  1926.291667  0.358497        Rogers Hornsby   \n",
       "4720  jacksjo01   4981  1772  1914.071429  0.355752  Shoeless Joe Jackson   \n",
       "465   barnero01   2391   860  1875.333333  0.359682           Ross Barnes   \n",
       "2377  delahed01   7510  2597  1895.500000  0.345806          Ed Delahanty   \n",
       "\n",
       "     bats  eb_estimate  \n",
       "1791    L     0.362864  \n",
       "4512    R     0.354065  \n",
       "4720    L     0.348884  \n",
       "465     R     0.345865  \n",
       "2377    R     0.341626  "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "batting.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use MCMC values for our beta binomial model:\n",
    "\n",
    "```\n",
    "Mu: 0.14255250930015972\n",
    "Sigma: 0.001828644970651173\n",
    "mu_ab: 0.015265514081094674\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 0.14255250930015972\n",
    "sigma = 0.001828644970651173\n",
    "mu_ab = 0.015265514081094674"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_0 = mu/sigma\n",
    "alpha_ab = mu_ab/sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "batting[\"alpha_0\"] = alpha_0 + alpha_ab * np.log(batting[\"AB\"])\n",
    "batting[\"beta_0\"]= 1/sigma - batting[\"alpha_0\"]\n",
    "batting[\"alpha_1\"] = batting[\"alpha_0\"] + batting[\"H\"]\n",
    "batting[\"beta_1\"] = batting[\"beta_0\"] + batting[\"AB\"] - batting[\"H\"]\n",
    "batting[\"new_eb\"] = batting[\"alpha_1\"] / (batting[\"alpha_1\"] + batting[\"beta_1\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>playerID</th>\n",
       "      <th>AB</th>\n",
       "      <th>H</th>\n",
       "      <th>yearID</th>\n",
       "      <th>average</th>\n",
       "      <th>name</th>\n",
       "      <th>bats</th>\n",
       "      <th>eb_estimate</th>\n",
       "      <th>alpha_0</th>\n",
       "      <th>beta_0</th>\n",
       "      <th>alpha_1</th>\n",
       "      <th>beta_1</th>\n",
       "      <th>new_eb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1791</th>\n",
       "      <td>cobbty01</td>\n",
       "      <td>11435</td>\n",
       "      <td>4189</td>\n",
       "      <td>1916.500000</td>\n",
       "      <td>0.366331</td>\n",
       "      <td>Ty Cobb</td>\n",
       "      <td>L</td>\n",
       "      <td>0.362864</td>\n",
       "      <td>155.962532</td>\n",
       "      <td>390.890474</td>\n",
       "      <td>4344.962532</td>\n",
       "      <td>7636.890474</td>\n",
       "      <td>0.362629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4512</th>\n",
       "      <td>hornsro01</td>\n",
       "      <td>8173</td>\n",
       "      <td>2930</td>\n",
       "      <td>1926.291667</td>\n",
       "      <td>0.358497</td>\n",
       "      <td>Rogers Hornsby</td>\n",
       "      <td>R</td>\n",
       "      <td>0.354065</td>\n",
       "      <td>153.158919</td>\n",
       "      <td>393.694088</td>\n",
       "      <td>3083.158919</td>\n",
       "      <td>5636.694088</td>\n",
       "      <td>0.353579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4720</th>\n",
       "      <td>jacksjo01</td>\n",
       "      <td>4981</td>\n",
       "      <td>1772</td>\n",
       "      <td>1914.071429</td>\n",
       "      <td>0.355752</td>\n",
       "      <td>Shoeless Joe Jackson</td>\n",
       "      <td>L</td>\n",
       "      <td>0.348884</td>\n",
       "      <td>149.024948</td>\n",
       "      <td>397.828058</td>\n",
       "      <td>1921.024948</td>\n",
       "      <td>3606.828058</td>\n",
       "      <td>0.347517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>465</th>\n",
       "      <td>barnero01</td>\n",
       "      <td>2391</td>\n",
       "      <td>860</td>\n",
       "      <td>1875.333333</td>\n",
       "      <td>0.359682</td>\n",
       "      <td>Ross Barnes</td>\n",
       "      <td>R</td>\n",
       "      <td>0.345865</td>\n",
       "      <td>142.898198</td>\n",
       "      <td>403.954808</td>\n",
       "      <td>1002.898198</td>\n",
       "      <td>1934.954808</td>\n",
       "      <td>0.341371</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2377</th>\n",
       "      <td>delahed01</td>\n",
       "      <td>7510</td>\n",
       "      <td>2597</td>\n",
       "      <td>1895.500000</td>\n",
       "      <td>0.345806</td>\n",
       "      <td>Ed Delahanty</td>\n",
       "      <td>R</td>\n",
       "      <td>0.341626</td>\n",
       "      <td>152.452674</td>\n",
       "      <td>394.400332</td>\n",
       "      <td>2749.452674</td>\n",
       "      <td>5307.400332</td>\n",
       "      <td>0.341256</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       playerID     AB     H       yearID   average                  name  \\\n",
       "1791   cobbty01  11435  4189  1916.500000  0.366331               Ty Cobb   \n",
       "4512  hornsro01   8173  2930  1926.291667  0.358497        Rogers Hornsby   \n",
       "4720  jacksjo01   4981  1772  1914.071429  0.355752  Shoeless Joe Jackson   \n",
       "465   barnero01   2391   860  1875.333333  0.359682           Ross Barnes   \n",
       "2377  delahed01   7510  2597  1895.500000  0.345806          Ed Delahanty   \n",
       "\n",
       "     bats  eb_estimate     alpha_0      beta_0      alpha_1       beta_1  \\\n",
       "1791    L     0.362864  155.962532  390.890474  4344.962532  7636.890474   \n",
       "4512    R     0.354065  153.158919  393.694088  3083.158919  5636.694088   \n",
       "4720    L     0.348884  149.024948  397.828058  1921.024948  3606.828058   \n",
       "465     R     0.345865  142.898198  403.954808  1002.898198  1934.954808   \n",
       "2377    R     0.341626  152.452674  394.400332  2749.452674  5307.400332   \n",
       "\n",
       "        new_eb  \n",
       "1791  0.362629  \n",
       "4512  0.353579  \n",
       "4720  0.347517  \n",
       "465   0.341371  \n",
       "2377  0.341256  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "batting.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "R    6103\n",
       "L    3004\n",
       "B     859\n",
       "Name: bats, dtype: int64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "batting.bats.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have some null values for handedness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 10800 entries, 1791 to 683\n",
      "Data columns (total 13 columns):\n",
      "playerID       10800 non-null object\n",
      "AB             10800 non-null int64\n",
      "H              10800 non-null int64\n",
      "yearID         10800 non-null float64\n",
      "average        10800 non-null float64\n",
      "name           10764 non-null object\n",
      "bats           9966 non-null object\n",
      "eb_estimate    10800 non-null float64\n",
      "alpha_0        10800 non-null float64\n",
      "beta_0         10800 non-null float64\n",
      "alpha_1        10800 non-null float64\n",
      "beta_1         10800 non-null float64\n",
      "new_eb         10800 non-null float64\n",
      "dtypes: float64(8), int64(2), object(3)\n",
      "memory usage: 1.2+ MB\n"
     ]
    }
   ],
   "source": [
    "batting.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use righthanded as default value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "batting[\"bats\"] = batting[\"bats\"].fillna(\"R\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One hot encode and drop indicator for right handed players:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "batting = pd.get_dummies(batting, columns=['bats'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `R` as the default value, hence we can drop that column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "batting.drop(\"bats_R\", inplace=True, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use good old `patsy` to create a spline and do some renaming:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transformed_batting = dmatrix(\"cr(train['yearID'],df = 3)\", {\"train\": batting[[\"yearID\"]]}, \n",
    "                              return_type='dataframe')\n",
    "transformed_batting.rename(columns={\"cr(train['yearID'], df=3)[0]\": \"year_1\",\n",
    "                                    \"cr(train['yearID'], df=3)[1]\": \"year_2\",\n",
    "                                    \"cr(train['yearID'], df=3)[2]\": \"year_3\"}, inplace=True)\n",
    "full_batting = batting.merge(transformed_batting[[\"year_1\", \"year_2\", \"year_3\"]], left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add some interactions using `patsy` again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "trans_batting = dmatrix(\"0+(year_1 + year_2 + year_3) * (bats_L + bats_B)\", \n",
    "                        full_batting, return_type='dataframe')\n",
    "new_batting = batting.merge(trans_batting[['year_1', 'year_2', 'year_3', 'year_1:bats_L',\n",
    "       'year_1:bats_B', 'year_2:bats_L', 'year_2:bats_B', 'year_3:bats_L',\n",
    "       'year_3:bats_B']], left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>playerID</th>\n",
       "      <th>AB</th>\n",
       "      <th>H</th>\n",
       "      <th>yearID</th>\n",
       "      <th>average</th>\n",
       "      <th>name</th>\n",
       "      <th>eb_estimate</th>\n",
       "      <th>alpha_0</th>\n",
       "      <th>beta_0</th>\n",
       "      <th>alpha_1</th>\n",
       "      <th>beta_1</th>\n",
       "      <th>new_eb</th>\n",
       "      <th>bats_B</th>\n",
       "      <th>bats_L</th>\n",
       "      <th>year_1</th>\n",
       "      <th>year_2</th>\n",
       "      <th>year_3</th>\n",
       "      <th>year_1:bats_L</th>\n",
       "      <th>year_1:bats_B</th>\n",
       "      <th>year_2:bats_L</th>\n",
       "      <th>year_2:bats_B</th>\n",
       "      <th>year_3:bats_L</th>\n",
       "      <th>year_3:bats_B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1791</th>\n",
       "      <td>cobbty01</td>\n",
       "      <td>11435</td>\n",
       "      <td>4189</td>\n",
       "      <td>1916.500000</td>\n",
       "      <td>0.366331</td>\n",
       "      <td>Ty Cobb</td>\n",
       "      <td>0.362864</td>\n",
       "      <td>155.962532</td>\n",
       "      <td>390.890474</td>\n",
       "      <td>4344.962532</td>\n",
       "      <td>7636.890474</td>\n",
       "      <td>0.362629</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.385753</td>\n",
       "      <td>0.805572</td>\n",
       "      <td>-0.191325</td>\n",
       "      <td>0.385753</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.805572</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.191325</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4512</th>\n",
       "      <td>hornsro01</td>\n",
       "      <td>8173</td>\n",
       "      <td>2930</td>\n",
       "      <td>1926.291667</td>\n",
       "      <td>0.358497</td>\n",
       "      <td>Rogers Hornsby</td>\n",
       "      <td>0.354065</td>\n",
       "      <td>153.158919</td>\n",
       "      <td>393.694088</td>\n",
       "      <td>3083.158919</td>\n",
       "      <td>5636.694088</td>\n",
       "      <td>0.353579</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.275787</td>\n",
       "      <td>0.920190</td>\n",
       "      <td>-0.195976</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4720</th>\n",
       "      <td>jacksjo01</td>\n",
       "      <td>4981</td>\n",
       "      <td>1772</td>\n",
       "      <td>1914.071429</td>\n",
       "      <td>0.355752</td>\n",
       "      <td>Shoeless Joe Jackson</td>\n",
       "      <td>0.348884</td>\n",
       "      <td>149.024948</td>\n",
       "      <td>397.828058</td>\n",
       "      <td>1921.024948</td>\n",
       "      <td>3606.828058</td>\n",
       "      <td>0.347517</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.414766</td>\n",
       "      <td>0.772548</td>\n",
       "      <td>-0.187314</td>\n",
       "      <td>0.414766</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.772548</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.187314</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>465</th>\n",
       "      <td>barnero01</td>\n",
       "      <td>2391</td>\n",
       "      <td>860</td>\n",
       "      <td>1875.333333</td>\n",
       "      <td>0.359682</td>\n",
       "      <td>Ross Barnes</td>\n",
       "      <td>0.345865</td>\n",
       "      <td>142.898198</td>\n",
       "      <td>403.954808</td>\n",
       "      <td>1002.898198</td>\n",
       "      <td>1934.954808</td>\n",
       "      <td>0.341371</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.937880</td>\n",
       "      <td>0.086291</td>\n",
       "      <td>-0.024171</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2377</th>\n",
       "      <td>delahed01</td>\n",
       "      <td>7510</td>\n",
       "      <td>2597</td>\n",
       "      <td>1895.500000</td>\n",
       "      <td>0.345806</td>\n",
       "      <td>Ed Delahanty</td>\n",
       "      <td>0.341626</td>\n",
       "      <td>152.452674</td>\n",
       "      <td>394.400332</td>\n",
       "      <td>2749.452674</td>\n",
       "      <td>5307.400332</td>\n",
       "      <td>0.341256</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.654586</td>\n",
       "      <td>0.472539</td>\n",
       "      <td>-0.127126</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>-0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       playerID     AB     H       yearID   average                  name  \\\n",
       "1791   cobbty01  11435  4189  1916.500000  0.366331               Ty Cobb   \n",
       "4512  hornsro01   8173  2930  1926.291667  0.358497        Rogers Hornsby   \n",
       "4720  jacksjo01   4981  1772  1914.071429  0.355752  Shoeless Joe Jackson   \n",
       "465   barnero01   2391   860  1875.333333  0.359682           Ross Barnes   \n",
       "2377  delahed01   7510  2597  1895.500000  0.345806          Ed Delahanty   \n",
       "\n",
       "      eb_estimate     alpha_0      beta_0      alpha_1       beta_1    new_eb  \\\n",
       "1791     0.362864  155.962532  390.890474  4344.962532  7636.890474  0.362629   \n",
       "4512     0.354065  153.158919  393.694088  3083.158919  5636.694088  0.353579   \n",
       "4720     0.348884  149.024948  397.828058  1921.024948  3606.828058  0.347517   \n",
       "465      0.345865  142.898198  403.954808  1002.898198  1934.954808  0.341371   \n",
       "2377     0.341626  152.452674  394.400332  2749.452674  5307.400332  0.341256   \n",
       "\n",
       "      bats_B  bats_L    year_1    year_2    year_3  year_1:bats_L  \\\n",
       "1791       0       1  0.385753  0.805572 -0.191325       0.385753   \n",
       "4512       0       0  0.275787  0.920190 -0.195976       0.000000   \n",
       "4720       0       1  0.414766  0.772548 -0.187314       0.414766   \n",
       "465        0       0  0.937880  0.086291 -0.024171       0.000000   \n",
       "2377       0       0  0.654586  0.472539 -0.127126       0.000000   \n",
       "\n",
       "      year_1:bats_B  year_2:bats_L  year_2:bats_B  year_3:bats_L  \\\n",
       "1791            0.0       0.805572            0.0      -0.191325   \n",
       "4512            0.0       0.000000            0.0      -0.000000   \n",
       "4720            0.0       0.772548            0.0      -0.187314   \n",
       "465             0.0       0.000000            0.0      -0.000000   \n",
       "2377            0.0       0.000000            0.0      -0.000000   \n",
       "\n",
       "      year_3:bats_B  \n",
       "1791           -0.0  \n",
       "4512           -0.0  \n",
       "4720           -0.0  \n",
       "465            -0.0  \n",
       "2377           -0.0  "
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_batting.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating named arrays\n",
    "year_1 = np.array(new_batting[\"year_1\"])\n",
    "year_2 = np.array(new_batting[\"year_2\"])\n",
    "year_3 = np.array(new_batting[\"year_3\"])\n",
    "y1l = np.array(new_batting[\"year_1:bats_L\"])\n",
    "y1b = np.array(new_batting[\"year_1:bats_B\"])\n",
    "y2l = np.array(new_batting[\"year_2:bats_L\"])\n",
    "y2b = np.array(new_batting[\"year_2:bats_B\"])\n",
    "y3l = np.array(new_batting[\"year_3:bats_L\"])\n",
    "y3b = np.array(new_batting[\"year_3:bats_B\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "# taken from https://stackoverflow.com/questions/54505173/finding-alpha-and-beta-of-beta-binomial-distribution-with-scipy-optimize-and-log\n",
    "from scipy.special import gammaln, logit, digamma\n",
    "\n",
    "def loglike_betabinom(params, *args):\n",
    "\n",
    "    k = args[0] # the OVERALL conversions\n",
    "    n = args[1] # the number of at-bats (AE)\n",
    "    b = args[2]\n",
    "    l = args[3]\n",
    "    y1 = args[4]\n",
    "    y2 = args[5]\n",
    "    y3 = args[6]\n",
    "    mu_y1 = params[5]\n",
    "    mu_y2 = params[6]\n",
    "    mu_y3 = params[7]\n",
    "    mu_y1b = params[8]\n",
    "    mu_y1l = params[9]\n",
    "    mu_y2b = params[10]\n",
    "    mu_y2l = params[11]\n",
    "    mu_y3b = params[12]\n",
    "    mu_y3l = params[13]\n",
    "\n",
    "    alpha = (params[0] + params[1] * np.log(n) + params[3] * b + params[4] * l\n",
    "             + mu_y1 * y1 + mu_y2 * y2 + mu_y1 * y3 + mu_y1b * b * y1 + mu_y1l * l * y1\n",
    "             + mu_y2b * b * y2 + mu_y2l * l * y2 + mu_y3b * b * y3 + mu_y3l * l * y3) / params[2], \n",
    "    beta = 1/params[2] - alpha\n",
    "\n",
    "    logpdf = gammaln(n+1) + gammaln(k+alpha) + gammaln(n-k+beta) + gammaln(alpha+beta) - \\\n",
    "     (gammaln(k+1) + gammaln(n-k+1) + gammaln(alpha) + gammaln(beta) + gammaln(n+alpha+beta))\n",
    "    \n",
    "    #return -np.sum(logpdf) \n",
    "    mask = np.isfinite(logpdf)\n",
    "    nll = -logpdf[mask].sum()\n",
    "    return nll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 37795.43338073921\n",
      " hess_inv: <14x14 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([ 21.27781045, 134.18029994, 117.78174667, -23.89206202,\n",
      "        28.14922482,   0.54715201,  20.78959369,   0.        ,\n",
      "        -1.22745405,  -1.38679752, -12.26944732,  25.2854079 ,\n",
      "       -10.35805326,   4.29354259])\n",
      "  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "     nfev: 1695\n",
      "      nit: 95\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([ 0.15414599,  0.01524164,  0.00162399,  0.00488026,  0.01199432,\n",
      "       -0.01458635, -0.01402161,  0.01      ,  0.01019966,  0.01746908,\n",
      "       -0.01034605, -0.00333706, -0.00384607, -0.01101001])\n"
     ]
    }
   ],
   "source": [
    "# the trick here is to set a non-zero lower bound for sigma and mu_ab\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "init_params = [0.2, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,]\n",
    "res = minimize(loglike_betabinom, x0=init_params,\n",
    "            args=(np.array(new_batting['H']), np.array(new_batting['AB']), np.array(new_batting[\"bats_B\"]), np.array(new_batting[\"bats_L\"]),\n",
    "                  year_1, year_2, year_3),\n",
    "            method='L-BFGS-B', options={'disp': True, 'maxiter': 500},\n",
    "              bounds=[(0.001, 0.4), (0.0001, 0.5), (0.001, 0.5), (-0.1, 0.1), (-0.1, 0.1),\n",
    "                      (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1),\n",
    "                      (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1), (-0.1, 0.1),])\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "(mu_0, mu_ab, sigma, mu_b, mu_l,\n",
    " mu_y1,mu_y2, mu_y3, mu_y1b, mu_y1l, mu_y2b,\n",
    " mu_y2l, mu_y3b, mu_y3l\n",
    ") = res.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_batting[\"alpha\"] = (mu_0 + mu_ab * np.log(new_batting[\"AB\"]) + mu_b * new_batting[\"bats_B\"] + mu_l * new_batting[\"bats_L\"]\n",
    "             + mu_y1 * year_1 + mu_y2 * year_2 + mu_y1 * year_3 + mu_y1b * new_batting[\"bats_B\"] * year_1 + mu_y1l * new_batting[\"bats_L\"] * year_1\n",
    "             + mu_y2b * new_batting[\"bats_B\"] * year_2 + mu_y2l * new_batting[\"bats_L\"] * year_2 + mu_y3b * new_batting[\"bats_B\"] * year_3\n",
    "                        + mu_y3l * new_batting[\"bats_L\"] * year_3) / sigma\n",
    "new_batting[\"beta\"] = 1/sigma - new_batting[\"alpha\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's calculate the new posterior batting averages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_batting[\"beta_2\"] = new_batting[\"beta\"] + new_batting[\"AB\"] - new_batting[\"H\"]\n",
    "new_batting[\"full_eb\"] = new_batting[\"alpha\"] / (new_batting[\"alpha\"] + new_batting[\"beta\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>playerID</th>\n",
       "      <th>AB</th>\n",
       "      <th>H</th>\n",
       "      <th>yearID</th>\n",
       "      <th>average</th>\n",
       "      <th>name</th>\n",
       "      <th>eb_estimate</th>\n",
       "      <th>alpha_0</th>\n",
       "      <th>beta_0</th>\n",
       "      <th>alpha_1</th>\n",
       "      <th>beta_1</th>\n",
       "      <th>new_eb</th>\n",
       "      <th>bats_B</th>\n",
       "      <th>bats_L</th>\n",
       "      <th>year_1</th>\n",
       "      <th>year_2</th>\n",
       "      <th>year_3</th>\n",
       "      <th>year_1:bats_L</th>\n",
       "      <th>year_1:bats_B</th>\n",
       "      <th>year_2:bats_L</th>\n",
       "      <th>year_2:bats_B</th>\n",
       "      <th>year_3:bats_L</th>\n",
       "      <th>year_3:bats_B</th>\n",
       "      <th>alpha</th>\n",
       "      <th>beta</th>\n",
       "      <th>beta_2</th>\n",
       "      <th>full_eb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1791</th>\n",
       "      <td>cobbty01</td>\n",
       "      <td>11435</td>\n",
       "      <td>4189</td>\n",
       "      <td>1916.500000</td>\n",
       "      <td>0.366331</td>\n",
       "      <td>Ty Cobb</td>\n",
       "      <td>0.362864</td>\n",
       "      <td>155.962532</td>\n",
       "      <td>390.890474</td>\n",
       "      <td>4344.962532</td>\n",
       "      <td>7636.890474</td>\n",
       "      <td>0.362629</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.385753</td>\n",
       "      <td>0.805572</td>\n",
       "      <td>-0.191325</td>\n",
       "      <td>0.385753</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.805572</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.191325</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>185.094214</td>\n",
       "      <td>430.674688</td>\n",
       "      <td>7676.674688</td>\n",
       "      <td>0.300590</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4512</th>\n",
       "      <td>hornsro01</td>\n",
       "      <td>8173</td>\n",
       "      <td>2930</td>\n",
       "      <td>1926.291667</td>\n",
       "      <td>0.358497</td>\n",
       "      <td>Rogers Hornsby</td>\n",
       "      <td>0.354065</td>\n",
       "      <td>153.158919</td>\n",
       "      <td>393.694088</td>\n",
       "      <td>3083.158919</td>\n",
       "      <td>5636.694088</td>\n",
       "      <td>0.353579</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.275787</td>\n",
       "      <td>0.920190</td>\n",
       "      <td>-0.195976</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>170.805055</td>\n",
       "      <td>444.963847</td>\n",
       "      <td>5687.963847</td>\n",
       "      <td>0.277385</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4720</th>\n",
       "      <td>jacksjo01</td>\n",
       "      <td>4981</td>\n",
       "      <td>1772</td>\n",
       "      <td>1914.071429</td>\n",
       "      <td>0.355752</td>\n",
       "      <td>Shoeless Joe Jackson</td>\n",
       "      <td>0.348884</td>\n",
       "      <td>149.024948</td>\n",
       "      <td>397.828058</td>\n",
       "      <td>1921.024948</td>\n",
       "      <td>3606.828058</td>\n",
       "      <td>0.347517</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0.414766</td>\n",
       "      <td>0.772548</td>\n",
       "      <td>-0.187314</td>\n",
       "      <td>0.414766</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.772548</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.187314</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>177.635826</td>\n",
       "      <td>438.133076</td>\n",
       "      <td>3647.133076</td>\n",
       "      <td>0.288478</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>465</th>\n",
       "      <td>barnero01</td>\n",
       "      <td>2391</td>\n",
       "      <td>860</td>\n",
       "      <td>1875.333333</td>\n",
       "      <td>0.359682</td>\n",
       "      <td>Ross Barnes</td>\n",
       "      <td>0.345865</td>\n",
       "      <td>142.898198</td>\n",
       "      <td>403.954808</td>\n",
       "      <td>1002.898198</td>\n",
       "      <td>1934.954808</td>\n",
       "      <td>0.341371</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.937880</td>\n",
       "      <td>0.086291</td>\n",
       "      <td>-0.024171</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>158.979334</td>\n",
       "      <td>456.789569</td>\n",
       "      <td>1987.789569</td>\n",
       "      <td>0.258180</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2377</th>\n",
       "      <td>delahed01</td>\n",
       "      <td>7510</td>\n",
       "      <td>2597</td>\n",
       "      <td>1895.500000</td>\n",
       "      <td>0.345806</td>\n",
       "      <td>Ed Delahanty</td>\n",
       "      <td>0.341626</td>\n",
       "      <td>152.452674</td>\n",
       "      <td>394.400332</td>\n",
       "      <td>2749.452674</td>\n",
       "      <td>5307.400332</td>\n",
       "      <td>0.341256</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0.654586</td>\n",
       "      <td>0.472539</td>\n",
       "      <td>-0.127126</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.000000</td>\n",
       "      <td>-0.0</td>\n",
       "      <td>169.855381</td>\n",
       "      <td>445.913521</td>\n",
       "      <td>5358.913521</td>\n",
       "      <td>0.275843</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       playerID     AB     H       yearID   average                  name  \\\n",
       "1791   cobbty01  11435  4189  1916.500000  0.366331               Ty Cobb   \n",
       "4512  hornsro01   8173  2930  1926.291667  0.358497        Rogers Hornsby   \n",
       "4720  jacksjo01   4981  1772  1914.071429  0.355752  Shoeless Joe Jackson   \n",
       "465   barnero01   2391   860  1875.333333  0.359682           Ross Barnes   \n",
       "2377  delahed01   7510  2597  1895.500000  0.345806          Ed Delahanty   \n",
       "\n",
       "      eb_estimate     alpha_0      beta_0      alpha_1       beta_1    new_eb  \\\n",
       "1791     0.362864  155.962532  390.890474  4344.962532  7636.890474  0.362629   \n",
       "4512     0.354065  153.158919  393.694088  3083.158919  5636.694088  0.353579   \n",
       "4720     0.348884  149.024948  397.828058  1921.024948  3606.828058  0.347517   \n",
       "465      0.345865  142.898198  403.954808  1002.898198  1934.954808  0.341371   \n",
       "2377     0.341626  152.452674  394.400332  2749.452674  5307.400332  0.341256   \n",
       "\n",
       "      bats_B  bats_L    year_1    year_2    year_3  year_1:bats_L  \\\n",
       "1791       0       1  0.385753  0.805572 -0.191325       0.385753   \n",
       "4512       0       0  0.275787  0.920190 -0.195976       0.000000   \n",
       "4720       0       1  0.414766  0.772548 -0.187314       0.414766   \n",
       "465        0       0  0.937880  0.086291 -0.024171       0.000000   \n",
       "2377       0       0  0.654586  0.472539 -0.127126       0.000000   \n",
       "\n",
       "      year_1:bats_B  year_2:bats_L  year_2:bats_B  year_3:bats_L  \\\n",
       "1791            0.0       0.805572            0.0      -0.191325   \n",
       "4512            0.0       0.000000            0.0      -0.000000   \n",
       "4720            0.0       0.772548            0.0      -0.187314   \n",
       "465             0.0       0.000000            0.0      -0.000000   \n",
       "2377            0.0       0.000000            0.0      -0.000000   \n",
       "\n",
       "      year_3:bats_B       alpha        beta       beta_2   full_eb  \n",
       "1791           -0.0  185.094214  430.674688  7676.674688  0.300590  \n",
       "4512           -0.0  170.805055  444.963847  5687.963847  0.277385  \n",
       "4720           -0.0  177.635826  438.133076  3647.133076  0.288478  \n",
       "465            -0.0  158.979334  456.789569  1987.789569  0.258180  \n",
       "2377           -0.0  169.855381  445.913521  5358.913521  0.275843  "
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_batting.head()"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}