Probabilistic+Programming+Workbook+-+Chapters+1+and+2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "\n",
    "%matplotlib inline\n",
    "from IPython.core.pylabtools import figsize\n",
    "from matplotlib import pyplot as plt\n",
    "figsize(11, 9)\n",
    "\n",
    "import scipy.stats as stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recreating count data example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "count_data = np.loadtxt('./Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/Chapter1_Introduction/data/txtdata.csv') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_texts = count_data.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    idx = np.arange(len(count_data))\n",
    "    tau = pm.DiscreteUniform(\"tau\", lower=0, upper=len(count_data) - 1)\n",
    "    lambda1 = pm.Exponential(\"lambda1\", 1./avg_texts)\n",
    "    lambda2 = pm.Exponential(\"lambda2\", 1./avg_texts)\n",
    "    lambda_ = pm.math.switch(tau >= idx, lambda1, lambda2)\n",
    "    count = pm.Poisson(\"count\", lambda_, observed=count_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(20000, step=step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "burnt_trace = trace[10000:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(311)\n",
    "plt.hist(burnt_trace['lambda1'], histtype='stepfilled', bins=35, normed=True)\n",
    "plt.xlim([15,30])\n",
    "plt.title(r\"Posterior of $\\lambda_1,\\lambda_2$ and $\\tau$\")\n",
    "\n",
    "ax = plt.subplot(312)\n",
    "plt.hist(burnt_trace['lambda2'], histtype='stepfilled', bins=35, normed=True)\n",
    "plt.xlim([15,30])\n",
    "\n",
    "ax = plt.subplot(313)\n",
    "plt.hist(burnt_trace['tau'], histtype='stepfilled', bins=35, normed=True)\n",
    "plt.xlabel('Switching day')\n",
    "plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recreating Bayesian A/B testing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_conversion_rate_A = 0.05\n",
    "true_conversion_rate_B = 0.03\n",
    "N_A = 1500\n",
    "N_B = 750\n",
    "\n",
    "observations_A = stats.bernoulli.rvs(true_conversion_rate_A, size=N_A)\n",
    "observations_B = stats.bernoulli.rvs(true_conversion_rate_B, size=N_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as ab_testing:\n",
    "    conversion_rate_A = pm.Uniform(\"conversion_rate_A\", lower=0, upper=1)\n",
    "    conversion_rate_B = pm.Uniform(\"conversion_rate_B\", lower=0, upper=1)\n",
    "    \n",
    "    delta = pm.Deterministic(\"delta\", conversion_rate_A - conversion_rate_B)\n",
    "    \n",
    "    conversions_A = pm.Bernoulli(\"conversions_A\", conversion_rate_A, observed=observations_A)\n",
    "    conversions_B = pm.Bernoulli(\"conversions_B\", conversion_rate_B, observed=observations_B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with ab_testing:\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(20000, step=step)\n",
    "    burnt_trace = trace[10000:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ax = plt.subplot(311)\n",
    "plt.title('Posterior distribution for conversion rate A')\n",
    "plt.hist(burnt_trace['conversion_rate_A'], histtype='stepfilled', normed=True)\n",
    "plt.axvline(0.05, color='black')\n",
    "\n",
    "ax = plt.subplot(312)\n",
    "plt.title('Posterior distribution for conversion rate A')\n",
    "plt.hist(burnt_trace['conversion_rate_B'], histtype='stepfilled', normed=True)\n",
    "plt.axvline(0.03, color='black')\n",
    "plt.tight_layout()\n",
    "ax = plt.subplot(313)\n",
    "plt.title('Posterior distribution for conversion rate $\\Delta$')\n",
    "plt.hist(burnt_trace['delta'], histtype='stepfilled', normed=True)\n",
    "plt.axvline(0.02, color='black')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = sum(burnt_trace['conversion_rate_A'] > burnt_trace['conversion_rate_B']) / len(burnt_trace['conversion_rate_A'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Probability of A being superior to B is {:.2%}\".format(p))"
   ]
  },
  {
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}