Quickstart_Guide_no_errors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imports and Set up\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import theano.tensor as T\n",
    "from numpy import random, sum as nsum, concatenate\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import theano"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up for a potential error\n",
    "from __future__ import absolute_import\n",
    "from __future__ import print_function\n",
    "from __future__ import unicode_literals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', 100)\n",
    "pd.set_option('display.max_rows', 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (15, 10)\n",
    "\n",
    "plt.rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear Regression\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up basic parameters\n",
    "num_features = 10\n",
    "num_observed = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Choose random values for the actual alpha and betas\n",
    "alpha_a = random.normal(size=1)\n",
    "\n",
    "betas_a = random.normal(size = num_features)\n",
    "\n",
    "# Create fake predictor data\n",
    "X_train = random.normal(size=(num_observed, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Calculate the actual data, but put a bit of noise in\n",
    "y_a = alpha_a + nsum(betas_a[None,:] * X_train, 1) + random.normal(size=(num_observed))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Strings vs Unicode:\n",
    "\n",
    "### Tip #1: Use the appropriate string type for your Python version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up the PyMC model\n",
    "lin_reg_model = pm.Model()\n",
    "with lin_reg_model:\n",
    "    \n",
    "    # CHANGE: All names are now cast as strings.\n",
    "    alpha = pm.Normal(b'alpha', mu=0, tau=10.**-2, shape=(1))\n",
    "    betas = pm.Normal(b'betas', mu=0, tau=10. ** -2, shape=(1, num_features))\n",
    "    \n",
    "    s = pm.HalfNormal(b's', tau=1)\n",
    "    \n",
    "    temp = alpha + T.dot(betas, X_train.T)\n",
    "\n",
    "    y = pm.Normal(b'y', mu=temp , tau=s ** -2, observed=y_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Sample from the model\n",
    "with lin_reg_model:\n",
    "    step = pm.NUTS()\n",
    "    nuts_trace = pm.sample(2e3, step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the trace\n",
    "# Note: You can specify a burn-in period by indexing appropriately\n",
    "pm.traceplot(nuts_trace[1000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Look at a summary of the trace\n",
    "pm.summary(nuts_trace[1000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "betas_a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear Regression: Now with Data!\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_boston\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn.metrics import r2_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "boston = load_boston()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "boston"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "boston_features = pd.DataFrame(boston.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "boston_features.columns = boston.feature_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "boston_features.CHAS = boston_features.CHAS.astype('bool')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "boston_features.RAD = boston_features.RAD.astype('int')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "boston_target = pd.DataFrame(boston.target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_train, X_test, Y_train, Y_test = train_test_split(\n",
    "    boston_features, boston_target, test_size=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lr = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lr.fit(X_train, Y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lr.score(X_test, Y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DataFrames don't work!\n",
    "\n",
    "### Tip #2: Turn Dataframes into numpy Arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Booleans don't work!\n",
    "\n",
    "### Tip #3: Cast Booleans as ints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# CHANGE: Boolean column is cast as an integer\n",
    "X_train.CHAS = X_train.CHAS.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_input = theano.shared(np.array(X_train))\n",
    "# CHANGE: Y_train is cast as a numpy array\n",
    "model_output = theano.shared(np.array(Y_train))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up the PyMC model\n",
    "lin_reg_model_w_data = pm.Model()\n",
    "with lin_reg_model_w_data:\n",
    "    \n",
    "    alpha = pm.Normal(b'alpha', mu=0, tau=10.**-2, shape=(1))\n",
    "    betas = pm.Normal(b'betas', mu=0, tau=10. ** -2, shape=(1, len(X_test.columns)))\n",
    "\n",
    "    s = pm.HalfNormal(b's', tau=1)\n",
    "    \n",
    "    temp = alpha + T.dot(model_input, betas.T)\n",
    "\n",
    "    y = pm.Normal(b'y', temp , tau=s ** -2, observed=model_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Note: If you don't specify a sampler, PyMC3 will choose one for you\n",
    "with lin_reg_model_w_data:\n",
    "    nuts_trace = pm.sample(8e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.traceplot(nuts_trace[5000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.summary(nuts_trace[5000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lr.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lr.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can score our model\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_test.CHAS = X_test.CHAS.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Replace shared variables with testing set\n",
    "model_input.set_value(np.array(X_test))\n",
    "model_output.set_value(np.array(Y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create posterior predictive samples\n",
    "ppc = pm.sample_ppc(nuts_trace[1000:], model=lin_reg_model_w_data, samples=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred = ppc['y'].mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r2_score(Y_test, pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Hiearchical Linear Regression\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up basic parameters\n",
    "num_cats = 4\n",
    "num_per_cat = 50\n",
    "num_observed = num_per_cat * num_cats\n",
    "num_features = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up categories\n",
    "cat = concatenate([[i] * num_per_cat for i in range(num_cats)])\n",
    "\n",
    "# Simulate the features\n",
    "features = random.normal(size=(num_observed, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Choose random values to represent the actual alphas and betas\n",
    "alpha_a = random.normal(size=(num_cats))\n",
    "beta_a = random.normal(size=(num_cats, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the actual data, but put a bit of noise in\n",
    "y = alpha_a[cat] + nsum(beta_a[cat] * features, 1) + random.normal(size=(num_observed))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set up the PyMC model\n",
    "hlm = pm.Model()\n",
    "with hlm:\n",
    "    \n",
    "    # Both alpha and beta are drawn from similar distributions\n",
    "    mu_alpha = pm.Normal(b\"mu_alpha\", mu=0, sd=10, shape=(1))\n",
    "    sigma_alpha = pm.Uniform(b\"sigma_alpha\", .0, 10, testval=2.)\n",
    "    \n",
    "    mu_beta = pm.Normal(b\"mu_beta\", mu=0, sd=10, shape=(1))\n",
    "    sigma_beta = pm.Uniform(b\"sigma_beta\", .0, 10, testval=2.)\n",
    "    \n",
    "    # Simulate the alphas\n",
    "    alpha = pm.Normal(b\"alpha\", mu_alpha, sigma_alpha, shape=(num_cats))\n",
    "    \n",
    "    # Simulate the betas\n",
    "    beta = pm.Normal(b'beta', mu_beta, sigma_beta, shape=(num_cats, num_features))\n",
    "    \n",
    "    c = T.constant(cat)\n",
    "    \n",
    "    # Simulate the noise\n",
    "    s = pm.Uniform(b\"s\", .01, 10, shape=num_cats)\n",
    "\n",
    "    yd = pm.Normal(b'y', mu=alpha[c] + T.sum(beta[c]*features, 1), tau=s[c] ** -2, observed=y)\n",
    "\n",
    "    step = pm.NUTS()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Actually sample the model\n",
    "# Note: you can do this either via the \"with\" statement or by specifying the model\n",
    "tr = pm.sample(3e4, step, model=hlm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plot the variables\n",
    "pm.traceplot(tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# See a summary\n",
    "pm.summary(tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Logistic Regression\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def numpy_invlogit(x):\n",
    "    return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up basic parameters\n",
    "num_features = 10\n",
    "num_observed = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Choose random values for the actual alpha and betas\n",
    "alpha_a = random.normal(size=1)\n",
    "\n",
    "betas_a = random.normal(size = num_features)\n",
    "\n",
    "# Create fake predictor data\n",
    "X_train = random.normal(size=(num_observed, num_features))\n",
    "X_test =  random.normal(size=(num_observed, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the outcomes\n",
    "Y_train = random.binomial(1, numpy_invlogit(alpha_a + np.sum(betas_a[None, :] * X_train, 1)))\n",
    "Y_test = random.binomial(1, numpy_invlogit(alpha_a + np.sum(betas_a[None, :] * X_test, 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model_input = theano.shared(X_train)\n",
    "model_output = theano.shared(Y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "log_reg_model = pm.Model()\n",
    "\n",
    "with log_reg_model:\n",
    "    alpha = pm.Normal(b'alpha', mu=0, tau=10.**-2, shape=(1))\n",
    "    betas = pm.Normal(b'betas', mu=0, tau=10. ** -2, shape=(1, num_features))\n",
    "    \n",
    "    p = pm.invlogit(alpha + T.sum(betas*model_input, 1))\n",
    "\n",
    "    o = pm.Bernoulli(b'o', p, observed=model_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tip #4: Make sure there are no NaNs or infinities in your data or you'll get weird results with ADVI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use ADVI\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with log_reg_model:\n",
    "    v_params = pm.variational.advi(n=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(v_params.elbo_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with log_reg_model:\n",
    "    advi_trace = pm.variational.sample_vp(v_params, draws=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with log_reg_model:\n",
    "    \n",
    "    step = pm.NUTS(scaling=v_params.stds)\n",
    "\n",
    "    nuts_trace = pm.sample(8e3, step, start=v_params.means)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.traceplot(nuts_trace[1000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.summary(nuts_trace[1000:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Replace shared variables with testing set\n",
    "model_input.set_value(X_test)\n",
    "model_output.set_value(Y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create posterior predictive samples from ADVI\n",
    "ppc = pm.sample_ppc(advi_trace, model=log_reg_model, samples=1000)\n",
    "\n",
    "# Use probability of > 0.5 to assume prediction of class 1\n",
    "pred = ppc['o'].mean(axis=0) > 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('ADVI Accuracy = {}%'.format((Y_test == pred).mean() * 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create posterior predictive samples from NUTS\n",
    "ppc = pm.sample_ppc(nuts_trace[1000:], model=log_reg_model, samples=1000)\n",
    "\n",
    "# Use probability of > 0.5 to assume prediction of class 1\n",
    "pred = ppc['o'].mean(axis=0) > 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('NUTS Accuracy = {}%'.format((Y_test == pred).mean() * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "API-ify your model\n",
    "==="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note: you can save your trace and reload it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import pickle\n",
    "# fileObject = open(\"advi_trace.pickle\",'w')  \n",
    "#pickle.dump(advi_trace, fileObject)\n",
    "#fileObject.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# advi_trace = pickle.load(open(\"pickle_jar/advi_trace.pickle\",'r')  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First, create some test data\n",
    "API_test =  random.normal(size=(1, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "API_Y_test = random.binomial(1, numpy_invlogit(alpha_a + np.sum(betas_a[None, :] * API_test, 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "API_Y_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_input = theano.shared(API_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tip #5: if you train your model with observed data, you need to put in fake observed data to sample from it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Put in some fake data\n",
    "API_fake_Y = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "API_model = pm.Model()\n",
    "\n",
    "with API_model:\n",
    "    alpha = pm.Normal(b'alpha', mu=0, tau=2.**-2, shape=(1))\n",
    "    betas = pm.Normal(b'betas', mu=0, tau=2. ** -2, shape=(1, num_features))\n",
    "    \n",
    "    p = pm.invlogit(alpha + T.sum(betas*model_input, 1))\n",
    "\n",
    "    # CHANGE: The outcome is now set with an observed value of fake data\n",
    "    o = pm.Bernoulli(b'o', p, observed=API_fake_Y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Create posterior predictive samples\n",
    "ppc = pm.sample_ppc(advi_trace, model=API_model, samples=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Use probability of > 0.5 to assume prediction of class 1\n",
    "pred = ppc['o'].mean(axis=0) > 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hiearchical Logistic Regression\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up basic parameters\n",
    "num_cats = 4\n",
    "\n",
    "# Need lots of data to converge\n",
    "num_per_cat = 15000\n",
    "num_observed = num_per_cat * num_cats\n",
    "num_features = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up categories\n",
    "cat = concatenate([[i] * num_per_cat for i in range(num_cats)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Simulate the features\n",
    "features = np.random.normal(size=(num_observed, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "alpha_a = np.random.normal(size=(num_cats))\n",
    "beta_a = np.random.normal(size=(num_cats, num_features))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the actual data\n",
    "p = alpha_a[cat] + nsum(beta_a[cat] * features, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "outcomes = np.random.binomial(1, numpy_invlogit(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hlr = pm.Model()\n",
    "\n",
    "with hlr:\n",
    "    # Both alpha and beta are drawn for the same distributions\n",
    "    mu_alpha = pm.Normal(b\"mu_alpha\", mu=0, sd=10, shape=(1))\n",
    "    sigma_alpha = pm.Uniform(b\"sigma_alpha\", .0, 10, testval=2.)\n",
    "    \n",
    "    mu_beta = pm.Normal(b\"mu_beta\", mu=0, sd=10, shape=(1))\n",
    "    sigma_beta = pm.Uniform(b\"sigma_beta\", 0, 10, testval=2.)\n",
    "    \n",
    "    alpha = pm.Normal(b'alpha', mu=mu_alpha, tau=sigma_alpha, shape=(num_cats))\n",
    "    beta = pm.Normal(b'beta', mu=mu_beta, tau=sigma_beta, shape=(num_cats, num_features))\n",
    "    \n",
    "    c = T.constant(cat)\n",
    "\n",
    "    p = pm.invlogit(alpha[c] + T.sum(beta[c]*features, 1))\n",
    "    \n",
    "    o = pm.Bernoulli(b'o', p, observed=outcomes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with hlr:\n",
    "    v_params = pm.variational.advi(n=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(v_params.elbo_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with hlr:\n",
    "    \n",
    "    step = pm.NUTS(scaling=v_params.stds)\n",
    "    trace = pm.sample(20000, step, start=v_params.means)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tip #6: You can stop your trace part way through and still use it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.traceplot(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.summary(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.forestplot(trace, varnames=['mu_alpha', 'mu_beta', 'alpha', 'sigma_alpha', 'sigma_beta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}