prototype for bayesian curve regression

Author: ericmjl

notebooks/09-bayesian-curve-regression.ipynb

@@ -0,0 +1,113 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pymc3 as pm\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from utils import ECDF\n",
+    "from data import load_decay\n",
+    "import pandas as pd\n",
+    "import theano.tensor as tt\n",
+    "\n",
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "%matplotlib inline\n",
+    "%config InlineBackend.figure_format = 'retina'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.random.seed(42)\n",
+    "\n",
+    "df = load_decay()\n",
+    "df.head(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ax = df['activity'].plot()\n",
+    "# ax.set_ylim(30, 40)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with pm.Model() as model:\n",
+    "    A = pm.HalfNormal('A', sd=100)\n",
+    "    tau = pm.Exponential('tau', lam=1)\n",
+    "    C = pm.Normal('C', sd=100)\n",
+    "    \n",
+    "    sd = pm.HalfCauchy('sd', beta=1)\n",
+    "    \n",
+    "    link = A * np.exp(-df['t'].values / tau) + C\n",
+    "    \n",
+    "    like = pm.Normal('activity', mu=link, sd=sd, observed=df['activity'].values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "with model:\n",
+    "    trace = pm.sample(2000, tune=2000)\n",
+    "    # Note: Sampler may pause for a while after finishing"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "traces = pm.traceplot(trace)\n",
+    "# despine_trace(traces)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "bayesian-modelling-tutorial",
+   "language": "python",
+   "name": "bayesian-modelling-tutorial"
+  },
+  "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.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

notebooks/data.py

@@ -8,6 +8,7 @@ def load_finches_2012():
     path = '../data/finch_beaks_2012.csv'
     return load_finches(path)
 
+
 def load_finches_1975():
     path = '../data/finch_beaks_1975.csv'
     return load_finches(path)
@@ -50,3 +51,20 @@ def load_kruschke():
     df = pd.read_csv('../data/iq.csv', index_col=0)  # comment out the path to the file for students.
     df = jn.DataFrame(df).label_encode('treatment')
     return df
+
+
+# Constants for load_decay
+tau = 71.9 # indium decay half life
+A = 60  # starting magnitude
+C = 2 # measurement error
+noise_scale = 1
+
+
+def load_decay():
+    t = np.arange(0, 1000)
+    def decay_func(ts, noise):
+        return A  * np.exp(-t/tau) + C + np.random.normal(0, noise, size=(len(t)))
+
+    data = {'t': t, 'activity': decay_func(t, noise_scale)}
+    df = pd.DataFrame(data)
+    return df