logistic-regression-advi.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 0 [0%]: ELBO = -190906.77\n",
      "Iteration 1000 [10%]: ELBO = -118799.35\n",
      "Iteration 2000 [20%]: ELBO = -112010.98\n",
      "Iteration 3000 [30%]: ELBO = -60375.62\n",
      "Iteration 4000 [40%]: ELBO = -53889.02\n",
      "Iteration 5000 [50%]: ELBO = -38572.64\n",
      "Iteration 6000 [60%]: ELBO = -38198.56\n",
      "Iteration 7000 [70%]: ELBO = -38449.0\n",
      "Iteration 8000 [80%]: ELBO = -29149.07\n",
      "Iteration 9000 [90%]: ELBO = -29631.96\n",
      "Finished [100%]: ELBO = -30219.02\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys, os\n",
    "sys.path.insert(0, os.path.expanduser('~/work/git/github/taku-y/pymc3'))\n",
    "\n",
    "import numpy as np\n",
    "from numpy import random\n",
    "from scipy import stats\n",
    "import pymc3 as pm\n",
    "from pymc3 import Model, Normal, Bernoulli, variational\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import theano.tensor as tt\n",
    "\n",
    "def invlogit(x):\n",
    "    return np.exp(x) / (1 + np.exp(x))\n",
    "\n",
    "def tinvlogit(x):\n",
    "    lower = 1e-6\n",
    "    upper = 1 - 1e-6\n",
    "    return lower + (upper - lower) * tt.exp(x) / (1 + tt.exp(x))\n",
    "\n",
    "# Set up basic parameters\n",
    "num_features = 20\n",
    "n = 100000\n",
    "\n",
    "# Choose random values for the actual alpha and betas\n",
    "rng = np.random.RandomState(10)\n",
    "alpha_a = rng.normal(size=1)\n",
    "\n",
    "betas_a = rng.normal(size = num_features)\n",
    "\n",
    "# Create fake predictor data\n",
    "predictors = rng.normal(size=(n, num_features))\n",
    "\n",
    "# Calculate the outcomes\n",
    "outcomes = rng.binomial(1, invlogit(alpha_a + np.sum(betas_a[None, :] * predictors, 1)))\n",
    "\n",
    "model = Model()\n",
    "\n",
    "with model:\n",
    "    alpha = Normal('alpha', mu=0, tau=2.**-2, shape=(1))\n",
    "    betas = Normal('betas', mu=0, tau=2. ** -2, shape=(1, num_features))\n",
    "\n",
    "    p = tinvlogit(alpha+tt.sum(betas * predictors, 1))\n",
    "\n",
    "    o = Bernoulli('o', p, observed=outcomes)\n",
    "\n",
    "with model:\n",
    "    advi_fit = variational.advi(n=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  },
  "widgets": {
   "state": {},
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}