diff --git a/docs/source/notebooks/best-api.ipynb b/docs/source/notebooks/best-api.ipynb
new file mode 100644
index 0000000000..0b89100086
--- /dev/null
+++ b/docs/source/notebooks/best-api.ipynb
@@ -0,0 +1,398 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Introduction\n",
+ "\n",
+ "This notebook is intended as an introduction to the BEST model available as a \"default model\" in PyMC3.\n",
+ "\n",
+ "John Kruschke introduced the BEST model in 2013 as a way to show how Bayesian estimation can provide a superior way of summarizing a two-sample comparison, i.e. treatment vs. control. Here, I have done the necessary software engineering to extend it to a multi-sample comparison, i.e. multiple treatments vs. control."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import pymc3 as pm\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "import matplotlib.pyplot as plt\n",
+ "import seaborn as sns\n",
+ "\n",
+ "sns.set_style('white')\n",
+ "\n",
+ "from pymc3.models.best import BEST\n",
+ "\n",
+ "%load_ext autoreload\n",
+ "%autoreload 2"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We will first start by generating some synthetic data. Here, we are assuming that there are:\n",
+ "\n",
+ "- 10 treatments (including the control)\n",
+ "- 4 replicate measurements each"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "n_treatments = 10\n",
+ "n_replicates = 4\n",
+ "means = np.random.randint(low=10, high=50, size=n_treatments)\n",
+ "sigmas = np.random.random(size=n_treatments) * 3\n",
+ "\n",
+ "data = []\n",
+ "indices = []\n",
+ "for idx, tmt in enumerate(range(n_treatments)):\n",
+ " for rep in range(n_replicates):\n",
+ " data.append(np.random.normal(loc=means[idx], scale=sigmas[idx]))\n",
+ " indices.append(idx)\n",
+ "\n",
+ "# Cast as a numpy array for convenience.\n",
+ "data = np.array(data)\n",
+ "indices = np.array(indices)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([ 21.51804962, 19.17406042, 16.51748568, 21.1254247 ,\n",
+ " 15.98845044, 16.01062808, 16.03388728, 15.98152944,\n",
+ " 11.55567225, 14.61522226, 13.82675371, 16.480401 ,\n",
+ " 38.81548662, 36.42769763, 38.78238235, 40.44465029,\n",
+ " 15.08502181, 14.0348035 , 10.01033688, 15.71514256,\n",
+ " 30.59898852, 29.99291836, 27.70848322, 27.48400447,\n",
+ " 24.09975511, 23.82401982, 21.77444056, 21.37802053,\n",
+ " 29.6104827 , 30.32411378, 35.00388415, 30.21269507,\n",
+ " 33.38952711, 35.06357773, 32.66744472, 32.57797696,\n",
+ " 39.77090017, 35.94262292, 38.0648428 , 37.21283409])"
+ ]
+ },
+ "execution_count": 3,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "data"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5,\n",
+ " 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9])"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "indices"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The BEST model accepts a `pandas` DataFrame object, and assumes that the data are structured as such:\n",
+ "\n",
+ "- Each row is one measurement\n",
+ "- Replicates are indicated by an identifying \"name\" (e.g. \"control\", \"treatment_1\", \"my_treatment_name\") etc.\n",
+ "- The replicate measurements do not necessarily have to be in order.\n",
+ "\n",
+ "While in this example we only use 10 different treatments, in principle it is possible to go much higher, even on the order of ~103 to ~104, as might be the case with replicate high throughput measurements."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " treatments | \n",
+ " measurements | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 0 | \n",
+ " 21.518050 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 0 | \n",
+ " 19.174060 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " 0 | \n",
+ " 16.517486 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " 0 | \n",
+ " 21.125425 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " 1 | \n",
+ " 15.988450 | \n",
+ "
\n",
+ " \n",
+ " | 5 | \n",
+ " 1 | \n",
+ " 16.010628 | \n",
+ "
\n",
+ " \n",
+ " | 6 | \n",
+ " 1 | \n",
+ " 16.033887 | \n",
+ "
\n",
+ " \n",
+ " | 7 | \n",
+ " 1 | \n",
+ " 15.981529 | \n",
+ "
\n",
+ " \n",
+ " | 8 | \n",
+ " 2 | \n",
+ " 11.555672 | \n",
+ "
\n",
+ " \n",
+ " | 9 | \n",
+ " 2 | \n",
+ " 14.615222 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "