| Software | Version |
|---|
| Python | 3.6.8 64bit [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] |
| IPython | 7.2.0 |
| OS | Darwin 18.7.0 x86_64 i386 64bit |
| pandas | 0.23.4 |
| numpy | 1.15.4 |
| Thu Jan 14 23:56:45 2021 EST |
"
],
"text/latex": [
"\\begin{tabular}{|l|l|}\\hline\n",
"{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
"Python & 3.6.8 64bit [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE\\_401/final)] \\\\ \\hline\n",
"IPython & 7.2.0 \\\\ \\hline\n",
"OS & Darwin 18.7.0 x86\\_64 i386 64bit \\\\ \\hline\n",
"pandas & 0.23.4 \\\\ \\hline\n",
"numpy & 1.15.4 \\\\ \\hline\n",
"\\hline \\multicolumn{2}{|l|}{Thu Jan 14 23:56:45 2021 EST} \\\\ \\hline\n",
"\\end{tabular}\n"
],
"text/plain": [
"Software versions\n",
"Python 3.6.8 64bit [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]\n",
"IPython 7.2.0\n",
"OS Darwin 18.7.0 x86_64 i386 64bit\n",
"pandas 0.23.4\n",
"numpy 1.15.4\n",
"Thu Jan 14 23:56:45 2021 EST"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%load_ext version_information\n",
"%version_information pandas, numpy"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Alex Galea \n",
"\n",
"CPython 3.6.8\n",
"IPython 7.2.0\n",
"\n",
"pymc3 3.6\n",
"matplotlib 3.0.2\n",
"seaborn 0.9.0\n",
"pandas 0.23.4\n"
]
}
],
"source": [
"%load_ext watermark\n",
"%watermark -a \"Alex Galea\" -v -p pymc3,matplotlib,seaborn,pandas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bayesian Modeling Discussion\n",
"\n",
"We can model the probability of an outcome $y$ as $P_t(y)$ using a **Gamma distribution** where time $t$ is defined below.\n",
"\n",
"The PMF is given by\n",
"\n",
"$$\n",
"P_t(\\alpha, \\beta) = \\frac{\\beta^\\alpha t^{\\alpha-1}e^{-\\beta t}}{\\Gamma (\\alpha)}\n",
"$$\n",
"\n",
"where $t$ is the time metric and $\\alpha$ & $\\beta$ are solved for with MCMC.\n",
"\n",
"Based on a set of goalie pull observations $X$ from 2003-2019 NHL games, we'll solve for the posterior distribution $P_t(y|X)$, the probability of the outcome $y$, given the observations. This is done computationally using markov chain monte carlo and the `pymc3` library.\n",
"\n",
"The outcomes we're interested in are $y = \\big\\{\\mathrm{goal\\;for}, \\mathrm{goal\\;against}, \\mathrm{no\\;goal}\\big\\}$. \n",
"\n",
"We'll use a **uniform prior** over the domain of times (last 5mins). Note: when gathering the observations, we throw out goalie pulls greater than 5 minutes from the end of the game (due to high likelihood of false positives when parsing goalie pulls from the raw game table).\n",
"\n",
"Once we find the posteriors discussed above, we can study the risk reward of pulling a goalie. We'll compare posteriors to find the odds of scoring a goal (and the odds of getting scored on) over time $t$ where:\n",
" - **t = Time elapsed** e.g. if there's 3 minutes left, what is the chance that pulling the goalie will result in a goal for?\n",
" - **t = Time since goalie pull** e.g. after the goalie has been pulled for 1 minute, what is the chance of getting a goal?"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"HALF_SPACE = u'\\u2000'"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"import pymc3 as pm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load the training data"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"20032004_goalie_pulls_2019-03-01.pkl 20122013_goalie_pulls_2019-04-25.pkl\r\n",
"20052006_goalie_pulls_2019-03-01.pkl 20132014_goalie_pulls_2019-04-25.pkl\r\n",
"20062007_goalie_pulls_2019-03-01.pkl 20142015_goalie_pulls_2019-04-25.pkl\r\n",
"20072008_goalie_pulls_2019-04-25.pkl 20152016_goalie_pulls_2019-04-25.pkl\r\n",
"20082009_goalie_pulls_2019-04-25.pkl 20162017_goalie_pulls_2019-04-25.pkl\r\n",
"20092010_goalie_pulls_2019-04-25.pkl 20172018_goalie_pulls_2019-04-25.pkl\r\n",
"20102011_goalie_pulls_2019-04-25.pkl 20182019_goalie_pulls_2019-04-25.pkl\r\n",
"20112012_goalie_pulls_2019-04-25.pkl \u001b[1m\u001b[36mtmp\u001b[m\u001b[m/\r\n"
]
}
],
"source": [
"ls ../../data/processed/pkl/"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def load_data():\n",
" files = glob.glob('../../data/processed/pkl/*.pkl')\n",
" files = sorted(files)\n",
" print(files)\n",
" return pd.concat((pd.read_pickle(f) for f in files), sort=False, ignore_index=True)\n",
"\n",
"def clean_df(df):\n",
" _df = df.copy()\n",
" \n",
" len_0 = _df.shape[0]\n",
" print('Removing goal_for_time < 15 mins')\n",
" _df = _df[~(_df.goal_for_time < datetime.timedelta(seconds=15*60))]\n",
" print(f'Removed {len_0 - _df.shape[0]} total rows')\n",
" \n",
" if 'game_end_time' in df.columns:\n",
" len_0 = _df.shape[0]\n",
" print('Removing game_end_time < 15 mins')\n",
" _df = _df[~(_df.game_end_time < datetime.timedelta(seconds=60*15))]\n",
" print(f'Removed {len_0 - _df.shape[0]} total rows')\n",
"\n",
" return _df"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['../../data/processed/pkl/20032004_goalie_pulls_2019-03-01.pkl', '../../data/processed/pkl/20052006_goalie_pulls_2019-03-01.pkl', '../../data/processed/pkl/20062007_goalie_pulls_2019-03-01.pkl', '../../data/processed/pkl/20072008_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20082009_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20092010_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20102011_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20112012_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20122013_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20132014_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20142015_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20152016_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20162017_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20172018_goalie_pulls_2019-04-25.pkl', '../../data/processed/pkl/20182019_goalie_pulls_2019-04-25.pkl']\n",
"Removing goal_for_time < 15 mins\n",
"Removed 115 total rows\n",
"Removing game_end_time < 15 mins\n",
"Removed 0 total rows\n"
]
}
],
"source": [
"df = load_data()\n",
"df = clean_df(df)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def load_training_samples(\n",
" df,\n",
" cols,\n",
" masks=[],\n",
" dtype='timedelta64[s]'\n",
") -> np.ndarray:\n",
" '''\n",
" Return buckets of training data.\n",
" '''\n",
" if not masks:\n",
" masks = [None] * len(cols)\n",
" out = []\n",
" for col, m in zip(cols, masks):\n",
" if m is None:\n",
" d = df[col].dropna().astype(dtype).values\n",
" else:\n",
" d = df[col][m].dropna().astype(dtype).values\n",
" out.append(d)\n",
" print(f'Loaded {len(d)} samples for col {col}')\n",
"\n",
" out = np.array(out)\n",
" print(f'Training data shape = {out.shape}')\n",
" return out"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Model 1 - Time elapsed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load data"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loaded 1496 samples for col pull_time\n",
"Loaded 3736 samples for col pull_time\n",
"Loaded 5937 samples for col pull_time\n",
"Training data shape = (3,)\n"
]
}
],
"source": [
"# Load time of pull for eventual outcomes:\n",
"feature_names = ['goal_for', 'goal_against', 'no_goals']\n",
"\n",
"# Logic for loading the data\n",
"features = ['pull_time', 'pull_time', 'pull_time']\n",
"masks = [\n",
" ~(df.goal_for_time.isnull()),\n",
" ~(df.goal_against_time.isnull()),\n",
" ~(df.game_end_timedelta.isnull()),\n",
"]\n",
"training_samples = load_training_samples(df, features, masks)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We actuall want to fit the mirror version of the Gamma distributioin (about the y-axis). In order to do this with PyMC3, we will transform the data."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"MAX_TIME_SECONDS = 20*60\n",
"training_samples[0] = MAX_TIME_SECONDS - training_samples[0]\n",
"training_samples[1] = MAX_TIME_SECONDS - training_samples[1]\n",
"training_samples[2] = MAX_TIME_SECONDS - training_samples[2]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 81., 86., 91., 89., 119., 79., 100., 97., 55., 58.]),\n",
" array([16., 57., 32., 67., 67., 60., 82., 70., 52., 96.]),\n",
" array([ 2., 75., 132., 49., 81., 63., 69., 60., 8., 63.]))"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(training_samples[0][:10],\n",
"training_samples[1][:10],\n",
"training_samples[2][:10],)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['goal_for', 'goal_against', 'no_goals']"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"feature_names"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### PyMC3 Model"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Multiprocess sampling (4 chains in 4 jobs)\n",
"CompoundStep\n",
">Metropolis: [gamma_no_goal]\n",
">Metropolis: [gamma_goal_against]\n",
">Metropolis: [gamma_goal_for]\n",
">Metropolis: [alpha_no_goal]\n",
">Metropolis: [alpha_goal_against]\n",
">Metropolis: [alpha_goal_for]\n",
">Metropolis: [beta_no_goal]\n",
">Metropolis: [beta_goal_against]\n",
">Metropolis: [beta_goal_for]\n",
"Sampling 4 chains: 100%|██████████| 74000/74000 [01:04<00:00, 1147.70draws/s]\n",
"The estimated number of effective samples is smaller than 200 for some parameters.\n"
]
},
{
"data": {
"text/latex": [
"$$\n",
" \\begin{array}{rcl}\n",
" \\text{beta_goal_for} &\\sim & \\text{Uniform}(\\mathit{lower}=0.0,~\\mathit{upper}=100.0)\\\\\\text{beta_goal_against} &\\sim & \\text{Uniform}(\\mathit{lower}=0.0,~\\mathit{upper}=100.0)\\\\\\text{beta_no_goal} &\\sim & \\text{Uniform}(\\mathit{lower}=0.0,~\\mathit{upper}=100.0)\\\\\\text{alpha_goal_for} &\\sim & \\text{Uniform}(\\mathit{lower}=0.001,~\\mathit{upper}=10.0)\\\\\\text{alpha_goal_against} &\\sim & \\text{Uniform}(\\mathit{lower}=0.001,~\\mathit{upper}=10.0)\\\\\\text{alpha_no_goal} &\\sim & \\text{Uniform}(\\mathit{lower}=0.001,~\\mathit{upper}=10.0)\\\\\\text{gamma_goal_for} &\\sim & \\text{Gamma}(\\mathit{alpha}=\\text{alpha_goal_for},~\\mathit{beta}=\\text{beta_goal_for})\\\\\\text{p_goal_for} &\\sim & \\text{Deterministic}(\\text{Constant},~\\text{gamma_goal_for_log__})\\\\\\text{gamma_goal_against} &\\sim & \\text{Gamma}(\\mathit{alpha}=\\text{alpha_goal_against},~\\mathit{beta}=\\text{beta_goal_against})\\\\\\text{p_goal_against} &\\sim & \\text{Deterministic}(\\text{Constant},~\\text{gamma_goal_against_log__})\\\\\\text{gamma_no_goal} &\\sim & \\text{Gamma}(\\mathit{alpha}=\\text{alpha_no_goal},~\\mathit{beta}=\\text{beta_no_goal})\\\\\\text{p_no_goal} &\\sim & \\text{Deterministic}(\\text{Constant},~\\text{gamma_no_goal_log__})\\\\\\text{obs_goal_for} &\\sim & \\text{Gamma}(\\mathit{alpha}=\\text{alpha_goal_for},~\\mathit{beta}=\\text{beta_goal_for})\\\\\\text{obs_goal_against} &\\sim & \\text{Gamma}(\\mathit{alpha}=\\text{alpha_goal_against},~\\mathit{beta}=\\text{beta_goal_against})\\\\\\text{obs_no_goal} &\\sim & \\text{Gamma}(\\mathit{alpha}=\\text{alpha_no_goal},~\\mathit{beta}=\\text{beta_no_goal})\n",
" \\end{array}\n",
" $$"
],
"text/plain": [
"