8schools_compare.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "J = 8\n",
    "y = np.array([28.,  8., -3.,  7., -1.,  1., 18., 12.])\n",
    "sigma = np.array([15., 10., 16., 11.,  9., 11., 10., 18.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyMC3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import theano.tensor as tt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as Centered_eight:\n",
    "    mu = pm.Normal('mu', mu=0, sd=5)\n",
    "    log_tau = pm.Normal('log_tau', mu=5, sd=1)\n",
    "    theta_prime = pm.Normal('theta_prime', mu=0, sd=1, shape=J)\n",
    "    theta = mu + tt.exp(log_tau) * theta_prime\n",
    "    obs = pm.Normal('obs', mu=theta, sd=sigma, observed=y)\n",
    "    trace = pm.sample()\n",
    "\n",
    "theta3 = trace['mu'][:, np.newaxis] + np.exp(trace['log_tau'])[:, np.newaxis] * trace['theta_prime']\n",
    "theta3.mean(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyMC4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc4 as pm\n",
    "import tensorflow as tf # For Random Variable operation\n",
    "from tensorflow_probability import edward2 as ed \n",
    "from pymc4.inference.sampling.sample import sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pm.Model(num_schools=J, y=y, sigma=sigma )\n",
    "@model.define\n",
    "def process(cfg):\n",
    "    mu = ed.Normal(loc=0., scale=5., name=\"mu\")  # `mu` above\n",
    "    # Due to the lack of HalfCauchy distribution.\n",
    "    log_tau = ed.Normal(\n",
    "        loc=5., scale=1., name=\"log_tau\")  # `log(tau)` above\n",
    "    theta_prime = ed.Normal(\n",
    "        loc=tf.zeros(cfg.num_schools),\n",
    "        scale=tf.ones(cfg.num_schools),\n",
    "        name=\"theta_prime\")  # `theta_prime` above\n",
    "    theta = mu + tf.exp(\n",
    "        log_tau) * theta_prime  # `theta` above\n",
    "    y = ed.Normal(\n",
    "        loc=theta,\n",
    "        scale=np.float32(cfg.sigma),\n",
    "        name=\"y\")  # `y` above\n",
    "    \n",
    "    return y\n",
    "\n",
    "trace = pm.sample(model)\n",
    "theta4 = trace['mu'][:, np.newaxis] + np.exp(trace['log_tau'])[:, np.newaxis] * trace['theta_prime']\n",
    "\n",
    "theta4.mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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>PyMC3_mu</th>\n",
       "      <th>PyMC3_std</th>\n",
       "      <th>PyMC4_mu</th>\n",
       "      <th>PyMC4_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>12.948958</td>\n",
       "      <td>10.727330</td>\n",
       "      <td>12.726427</td>\n",
       "      <td>10.590535</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5.875544</td>\n",
       "      <td>7.512560</td>\n",
       "      <td>5.851380</td>\n",
       "      <td>7.713525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.379641</td>\n",
       "      <td>10.265549</td>\n",
       "      <td>0.752818</td>\n",
       "      <td>10.213339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5.330624</td>\n",
       "      <td>8.099070</td>\n",
       "      <td>5.279244</td>\n",
       "      <td>7.985706</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.804756</td>\n",
       "      <td>7.540049</td>\n",
       "      <td>0.577965</td>\n",
       "      <td>7.336870</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1.992809</td>\n",
       "      <td>8.381048</td>\n",
       "      <td>2.091519</td>\n",
       "      <td>8.090734</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>11.580781</td>\n",
       "      <td>8.296009</td>\n",
       "      <td>11.756276</td>\n",
       "      <td>8.154613</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6.023647</td>\n",
       "      <td>9.903032</td>\n",
       "      <td>5.516689</td>\n",
       "      <td>10.676764</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    PyMC3_mu  PyMC3_std   PyMC4_mu  PyMC4_std\n",
       "0  12.948958  10.727330  12.726427  10.590535\n",
       "1   5.875544   7.512560   5.851380   7.713525\n",
       "2   1.379641  10.265549   0.752818  10.213339\n",
       "3   5.330624   8.099070   5.279244   7.985706\n",
       "4   0.804756   7.540049   0.577965   7.336870\n",
       "5   1.992809   8.381048   2.091519   8.090734\n",
       "6  11.580781   8.296009  11.756276   8.154613\n",
       "7   6.023647   9.903032   5.516689  10.676764"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({\"PyMC3_mu\": theta3.mean(axis=0), \"PyMC4_mu\": theta4.mean(axis=0), \"PyMC3_std\": theta3.std(axis=0), \"PyMC4_std\": theta4.std(axis=0)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (pymc3)",
   "language": "python",
   "name": "pymc3"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}