diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index c1b7989..44f2fe8 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -51,3 +51,6 @@ jobs:
     - name: Normal Mean Hypothesis Testing
       run: |
         python test/hypothesis_test.py
+    - name: Delta Binomial Testing
+      run: |
+        python test/delta_binomial_test.py
diff --git a/bbai/model/__init__.py b/bbai/model/__init__.py
new file mode 100644
index 0000000..7c3e339
--- /dev/null
+++ b/bbai/model/__init__.py
@@ -0,0 +1,6 @@
+from ._delta_binomial import DeltaBinomialModel
+
+__all__ = [
+        'DeltaBinomialModel',
+]
+
diff --git a/bbai/model/_delta_binomial.py b/bbai/model/_delta_binomial.py
new file mode 100644
index 0000000..aaf927c
--- /dev/null
+++ b/bbai/model/_delta_binomial.py
@@ -0,0 +1,89 @@
+from bbai._computation._bridge import mdldb
+
+class DeltaBinomialModel:
+    """Model the difference of the probabilities of two binomial distributions.
+
+    The delta of binomial distribution probabilities was studied by Laplace in [1, p. 59], where
+    Laplace sought to determine whether London has a higher boys-to-girls birth rate than Paris.
+
+    The model provides a posterior distribution for the difference of two binomial distribution
+    probabilities. It can be fit with three different priors:
+
+       * uniform: the original uniform prior used by Laplace
+       * Jeffreys: Jeffreys rule prior
+       * Reference (recommended): the reference, as determined by Proposition 0.2 of [2].
+
+    References:
+      [1]: Laplace, P. (1778). Mémoire sur les probabilités. Translated by Richard J. Pulskamp.
+      [2]: Berger, J., J. Bernardo, and D. Sun (2022). Objective bayesian inference and its 
+           relationship to frequentism.
+    
+    Parameters
+    ----------
+    prior : default='reference'
+        The prior to use. Possible values are 'uniform', 'jeffreys', and 'reference'.
+
+    Example
+    --------
+    >>> from bbai.model import DeltaBinomialModel
+    >>> import numpy as np
+    >>> 
+    >>> a1, b1, a2, b2 = 5, 3, 2, 7
+    >>> model = DeltaBinomialModel(prior='reference')
+    >>>
+    >>> # Fit a posterior distribution with likelihood function
+    >>> #     L(theta, x) = (theta + x)^a1 * (1 - theta - x)^b1 * x^a2 (1-x)^b2 
+    >>> # where theta represents the difference of the two binomial distribution probabilities
+    >>> model.fit(a1, b1, a2, b2)
+    >>>
+    >>> # Print the probability that theta < 0.123
+    >>> print(model.cdf(0.123))
+    >>>     # Prints 0.10907436812863071
+    """
+    UNIFORM = 0
+    JEFFREYS = 1
+    REFERENCE = 2
+
+    def __init__(self, prior='reference'):
+        if prior == 'uniform':
+            self.prior_ = DeltaBinomialModel.UNIFORM
+        elif prior == 'jeffreys':
+            self.prior_ = DeltaBinomialModel.JEFFREYS
+        elif prior == 'reference':
+            self.prior_ = DeltaBinomialModel.REFERENCE
+        else:
+            raise RuntimeError("unknown prior {}".format(prior))
+        self.model_ = None
+
+    def fit(self, a1, b1, a2, b2):
+        """Fit a posterior distribution for the parameter p1 - p2 where p1 and p2
+        represent probabilities for binomial distributions with likelihood functions
+            L1(p1) = p1^a1 * (1 - p1)^b1
+            L2(p2) = p2^a2 * (1 - p2)^b2
+        """
+        self.model_ = mdldb.delta_binomial_fit(self.prior_, a1, b1, a2, b2)
+
+    def cdf(self, t):
+        """Cumulative distribution function for the posterior in p1-p2"""
+        assert -1 <= t and t <= 1
+        assert self.model_ is not None
+        return mdldb.delta_binomial_cdf(self.model_, t)
+
+    def pdf(self, t):
+        """Probability density function for the posterior in p1-p2"""
+        assert -1 <= t and t <= 1
+        assert self.model_ is not None
+        return mdldb.delta_binomial_pdf(self.model_, t)
+
+    def prior_pdf(self, *args):
+        """Probability density function for the prior in p1-p2
+            model.prior_pdf(theta)
+                 # returns the marginal prior PDF 
+            model.prior_pdf(theta, x)
+                 # returns the prior PDF 
+        """
+        assert self.model_ is not None
+        if len(args) == 1:
+            return mdldb.delta_binomial_prior_marginal_pdf(self.model_, args[0])
+        elif len(args) == 2:
+            return mdldb.delta_binomial_prior_pdf(self.model_, args[0], args[1])
diff --git a/example/21-delta-binomial-coverage.ipynb b/example/21-delta-binomial-coverage.ipynb
new file mode 100644
index 0000000..0878923
--- /dev/null
+++ b/example/21-delta-binomial-coverage.ipynb
@@ -0,0 +1,570 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Runs frequentist coverage testing for the difference of two binomial distribution probabilities using 95% posterior credible sets (see [2] for background).\n",
+    "\n",
+    "The delta of binomial distribution probabilities was studied by Laplace in [1, p. 59]. Seeking to determine whether London had a higher boys-to-girls birth rate than Paris, Laplace modeled birth rates as binomial distributions and computed a posterior for the delta of the two probabilities using a uniform prior.\n",
+    "\n",
+    "Of course, we know now that the uniform prior is arbitrary and dependent on the measurement scale. This notebook tests Laplace's original uniform prior against Jeffreys rule prior; and the modern approach of reference priors, Proposition 0.2 of [2].\n",
+    "\n",
+    "### References:\n",
+    "[1]: Laplace, P. (1778). Mémoire sur les probabilités. Translated by Richard J. Pulskamp.\n",
+    "\n",
+    "[2]: Berger, J., J. Bernardo, and D. Sun (2022). Objective bayesian inference and its \n",
+    "           relationship to frequentism."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Imports"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from bbai.model import DeltaBinomialModel\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from collections import defaultdict\n",
+    "import math\n",
+    "import pandas as pd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Set up coverage simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "alpha = 0.95\n",
+    "low = (1 - alpha) / 2.0\n",
+    "high = 1.0 - low\n",
+    "\n",
+    "def coverage(m, theta, x, n1, n2):\n",
+    "    p1 = theta + x\n",
+    "    p2 = x\n",
+    "    theta = p1 - p2\n",
+    "    res = 0.0\n",
+    "    for k1 in range(n1+1):\n",
+    "        prob1 = math.comb(n1, k1) * p1 ** k1 * (1 - p1)**(n1 - k1)\n",
+    "        for k2 in range(n2+1):\n",
+    "            prob = prob1 * math.comb(n2, k2) * p2 ** k2 * (1 - p2)**(n2 - k2)\n",
+    "            m.fit(k1, n1-k1, k2, n2-k2)\n",
+    "            cdf = m.cdf(theta)\n",
+    "            if low <= cdf and cdf <= high:\n",
+    "                res += prob\n",
+    "    return res"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def run_simulation(params, n1, n2):\n",
+    "    covs = []\n",
+    "    for prior in ['uniform', 'jeffreys', 'reference']:\n",
+    "        row = [prior]\n",
+    "        for theta, xs in params:\n",
+    "            for x in xs:\n",
+    "                m = DeltaBinomialModel(prior=prior)\n",
+    "                cov = coverage(m, theta, x, n1, n2)\n",
+    "                row.append(cov)\n",
+    "        covs.append(row)\n",
+    "    return covs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def make_coverage_table(params, covs):\n",
+    "    cols = []\n",
+    "    cols.append(('', 'prior'))\n",
+    "    for theta, xs in params:\n",
+    "        for x in xs:\n",
+    "            cols.append(('theta={}'.format(theta), 'x={}'.format(x)))\n",
+    "    df = pd.DataFrame(\n",
+    "        covs,\n",
+    "        columns = pd.MultiIndex.from_tuples(cols)\n",
+    "    )\n",
+    "    df = df.set_index(df.columns[0])\n",
+    "    df.index.name = 'prior'\n",
+    "    return df"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Run Simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "params = [\n",
+    "    (0.0, [0.0, 0.1, 0.25, 0.5]),\n",
+    "    (0.5, [0.0, 0.1, 0.25]),\n",
+    "    (0.98, [0.0, 0.005, 0.01]),\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### $n_1=3, n_2=3$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "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 tr th {\n",
+       "        text-align: left;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead tr:last-of-type th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr>\n",
+       "      <th></th>\n",
+       "      <th colspan=\"4\" halign=\"left\">theta=0.0</th>\n",
+       "      <th colspan=\"3\" halign=\"left\">theta=0.5</th>\n",
+       "      <th colspan=\"3\" halign=\"left\">theta=0.98</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th></th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.1</th>\n",
+       "      <th>x=0.25</th>\n",
+       "      <th>x=0.5</th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.1</th>\n",
+       "      <th>x=0.25</th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.005</th>\n",
+       "      <th>x=0.01</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>prior</th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>uniform</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.998542</td>\n",
+       "      <td>0.986816</td>\n",
+       "      <td>0.96875</td>\n",
+       "      <td>0.875</td>\n",
+       "      <td>0.927288</td>\n",
+       "      <td>0.949219</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.00000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>jeffreys</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.998542</td>\n",
+       "      <td>0.986816</td>\n",
+       "      <td>0.96875</td>\n",
+       "      <td>0.875</td>\n",
+       "      <td>0.927288</td>\n",
+       "      <td>0.949219</td>\n",
+       "      <td>0.941192</td>\n",
+       "      <td>0.941408</td>\n",
+       "      <td>0.94148</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>reference</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.998542</td>\n",
+       "      <td>0.986816</td>\n",
+       "      <td>0.96875</td>\n",
+       "      <td>1.000</td>\n",
+       "      <td>0.974160</td>\n",
+       "      <td>0.962402</td>\n",
+       "      <td>0.941192</td>\n",
+       "      <td>0.941408</td>\n",
+       "      <td>0.94148</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "          theta=0.0                              theta=0.5            \\\n",
+       "              x=0.0     x=0.1    x=0.25    x=0.5     x=0.0     x=0.1   \n",
+       "prior                                                                  \n",
+       "uniform         1.0  0.998542  0.986816  0.96875     0.875  0.927288   \n",
+       "jeffreys        1.0  0.998542  0.986816  0.96875     0.875  0.927288   \n",
+       "reference       1.0  0.998542  0.986816  0.96875     1.000  0.974160   \n",
+       "\n",
+       "                    theta=0.98                     \n",
+       "             x=0.25      x=0.0   x=0.005   x=0.01  \n",
+       "prior                                              \n",
+       "uniform    0.949219   0.000000  0.000000  0.00000  \n",
+       "jeffreys   0.949219   0.941192  0.941408  0.94148  \n",
+       "reference  0.962402   0.941192  0.941408  0.94148  "
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "covs = run_simulation(params, 3, 3)\n",
+    "make_coverage_table(params, covs)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### $n_1=10, n_2=3$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "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 tr th {\n",
+       "        text-align: left;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead tr:last-of-type th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr>\n",
+       "      <th></th>\n",
+       "      <th colspan=\"4\" halign=\"left\">theta=0.0</th>\n",
+       "      <th colspan=\"3\" halign=\"left\">theta=0.5</th>\n",
+       "      <th colspan=\"3\" halign=\"left\">theta=0.98</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th></th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.1</th>\n",
+       "      <th>x=0.25</th>\n",
+       "      <th>x=0.5</th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.1</th>\n",
+       "      <th>x=0.25</th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.005</th>\n",
+       "      <th>x=0.01</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>prior</th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>uniform</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.989656</td>\n",
+       "      <td>0.983693</td>\n",
+       "      <td>0.985596</td>\n",
+       "      <td>0.945312</td>\n",
+       "      <td>0.932969</td>\n",
+       "      <td>0.945239</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>jeffreys</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.989592</td>\n",
+       "      <td>0.978478</td>\n",
+       "      <td>0.956299</td>\n",
+       "      <td>0.945312</td>\n",
+       "      <td>0.960052</td>\n",
+       "      <td>0.951203</td>\n",
+       "      <td>0.817073</td>\n",
+       "      <td>0.846899</td>\n",
+       "      <td>0.877521</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>reference</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.894403</td>\n",
+       "      <td>0.928312</td>\n",
+       "      <td>0.948242</td>\n",
+       "      <td>0.944336</td>\n",
+       "      <td>0.962462</td>\n",
+       "      <td>0.951790</td>\n",
+       "      <td>0.983822</td>\n",
+       "      <td>0.975868</td>\n",
+       "      <td>0.966160</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "          theta=0.0                               theta=0.5            \\\n",
+       "              x=0.0     x=0.1    x=0.25     x=0.5     x=0.0     x=0.1   \n",
+       "prior                                                                   \n",
+       "uniform         1.0  0.989656  0.983693  0.985596  0.945312  0.932969   \n",
+       "jeffreys        1.0  0.989592  0.978478  0.956299  0.945312  0.960052   \n",
+       "reference       1.0  0.894403  0.928312  0.948242  0.944336  0.962462   \n",
+       "\n",
+       "                    theta=0.98                      \n",
+       "             x=0.25      x=0.0   x=0.005    x=0.01  \n",
+       "prior                                               \n",
+       "uniform    0.945239   0.000000  0.000000  0.000000  \n",
+       "jeffreys   0.951203   0.817073  0.846899  0.877521  \n",
+       "reference  0.951790   0.983822  0.975868  0.966160  "
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "covs = run_simulation(params, 10, 3)\n",
+    "make_coverage_table(params, covs)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### $n_1=10, n_2=10$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "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 tr th {\n",
+       "        text-align: left;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead tr:last-of-type th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr>\n",
+       "      <th></th>\n",
+       "      <th colspan=\"4\" halign=\"left\">theta=0.0</th>\n",
+       "      <th colspan=\"3\" halign=\"left\">theta=0.5</th>\n",
+       "      <th colspan=\"3\" halign=\"left\">theta=0.98</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th></th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.1</th>\n",
+       "      <th>x=0.25</th>\n",
+       "      <th>x=0.5</th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.1</th>\n",
+       "      <th>x=0.25</th>\n",
+       "      <th>x=0.0</th>\n",
+       "      <th>x=0.005</th>\n",
+       "      <th>x=0.01</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>prior</th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "      <th></th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>uniform</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.990960</td>\n",
+       "      <td>0.965160</td>\n",
+       "      <td>0.957809</td>\n",
+       "      <td>0.944336</td>\n",
+       "      <td>0.947102</td>\n",
+       "      <td>0.954700</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>jeffreys</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.949782</td>\n",
+       "      <td>0.915044</td>\n",
+       "      <td>0.947739</td>\n",
+       "      <td>0.934570</td>\n",
+       "      <td>0.947157</td>\n",
+       "      <td>0.934719</td>\n",
+       "      <td>0.817073</td>\n",
+       "      <td>0.817698</td>\n",
+       "      <td>0.817907</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>reference</th>\n",
+       "      <td>1.0</td>\n",
+       "      <td>0.949782</td>\n",
+       "      <td>0.915044</td>\n",
+       "      <td>0.947739</td>\n",
+       "      <td>0.978516</td>\n",
+       "      <td>0.950859</td>\n",
+       "      <td>0.934762</td>\n",
+       "      <td>0.983822</td>\n",
+       "      <td>0.983311</td>\n",
+       "      <td>0.983141</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "          theta=0.0                               theta=0.5            \\\n",
+       "              x=0.0     x=0.1    x=0.25     x=0.5     x=0.0     x=0.1   \n",
+       "prior                                                                   \n",
+       "uniform         1.0  0.990960  0.965160  0.957809  0.944336  0.947102   \n",
+       "jeffreys        1.0  0.949782  0.915044  0.947739  0.934570  0.947157   \n",
+       "reference       1.0  0.949782  0.915044  0.947739  0.978516  0.950859   \n",
+       "\n",
+       "                    theta=0.98                      \n",
+       "             x=0.25      x=0.0   x=0.005    x=0.01  \n",
+       "prior                                               \n",
+       "uniform    0.954700   0.000000  0.000000  0.000000  \n",
+       "jeffreys   0.934719   0.817073  0.817698  0.817907  \n",
+       "reference  0.934762   0.983822  0.983311  0.983141  "
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "covs = run_simulation(params, 10, 10)\n",
+    "make_coverage_table(params, covs)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/test/delta_binomial_test.py b/test/delta_binomial_test.py
new file mode 100644
index 0000000..66900cd
--- /dev/null
+++ b/test/delta_binomial_test.py
@@ -0,0 +1,83 @@
+import pytest
+import numpy as np
+from bbai.model import DeltaBinomialModel
+import scipy
+
+def test_delta_binomial_model():
+    for prior in ['uniform', 'jeffreys', 'reference']:
+        a1 = 2
+        b1 = 1
+        a2 = 3
+        b2 = 4
+        model = DeltaBinomialModel(prior=prior)
+        model.fit(a1, b1, a2, b2)
+
+        # cdf and pdf match up
+        t = 0.123
+        h = 1.0e-6
+        d = (model.cdf(t + h) - model.cdf(t)) / h
+        expected = model.pdf(t)
+        assert d == pytest.approx(expected, rel=1.0e-4)
+
+        d = (model.cdf(-t + h) - model.cdf(-t)) / h
+        expected = model.pdf(-t)
+        assert d == pytest.approx(expected, rel=1.0e-4)
+
+        # end point behavior
+        assert model.cdf(-1.0) == 0.0
+        assert model.cdf(1.0) == 1.0
+
+        # prior integrates to 1
+        delta = 1.0e-3
+        tot = 0.0
+        for t in np.arange(-0.9999, 0.9999, delta):
+            tot += model.prior_pdf(t)
+        assert tot * delta == pytest.approx(1.0, rel=1.0e-2)
+
+        # marginal prior matches up with pdf
+        delta = 1.0e-6
+        tot = 0.0
+        theta = 0.123
+        val = model.prior_pdf(theta)
+        eps = 1.0e-7
+        def f(x):
+          return model.prior_pdf(theta, x)
+        res = scipy.integrate.quad(f, eps, 1.0 - theta - eps)
+        assert res[0] == pytest.approx(val, rel=1.0e-3)
+
+        # symmetry
+        t = 0.123
+        model.fit(a1, b1, a2, b2)
+        val1 = model.cdf(t)
+        model.fit(b1, a1, b2, a2)
+        val2 = 1 - model.cdf(-t)
+        assert val1 == pytest.approx(val2)
+
+def test_precompute():
+    a1 = 3
+    b1 = 2
+    a2 = 1
+    b2 = 2
+
+    # uniform
+    model = DeltaBinomialModel(prior='uniform')
+    model.fit(a1, b1, a2, b2)
+    val = model.pdf(0.123)
+    assert val == pytest.approx(1.38587)
+
+    # Jeffreys
+    model = DeltaBinomialModel(prior='jeffreys')
+    model.fit(a1, b1, a2, b2)
+    val = model.pdf(0.123)
+    assert val == pytest.approx(1.22913, rel=1.0e-4)
+
+    # Reference
+    model = DeltaBinomialModel(prior='reference')
+    model.fit(a1, b1, a2, b2)
+    val = model.pdf(0.123)
+    assert val == pytest.approx(1.08699, rel=1.0e-4)
+
+
+if __name__ == "__main__":
+    raise SystemExit(pytest.main([__file__]))
+
diff --git a/wheel/bbai-1.6.0-py3-none-manylinux1_x86_64.whl b/wheel/bbai-1.7.1-py3-none-manylinux1_x86_64.whl
similarity index 67%
rename from wheel/bbai-1.6.0-py3-none-manylinux1_x86_64.whl
rename to wheel/bbai-1.7.1-py3-none-manylinux1_x86_64.whl
index d8c30ad..388a81c 100644
Binary files a/wheel/bbai-1.6.0-py3-none-manylinux1_x86_64.whl and b/wheel/bbai-1.7.1-py3-none-manylinux1_x86_64.whl differ