Confinterval

notebooks/Module 1.5.2 - Bayesian Inference Basics.ipynb

@@ -428,7 +428,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 1,
    "metadata": {
     "attributes": {
      "classes": [],
@@ -531,22 +531,23 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "0.3999999999999999"
+       "0.75"
       ]
      },
-     "execution_count": 23,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "prior(0.8, 1, 2)  # In a beta distribution with params (1, 2), there is a 40% chance of seeing 80% heads"
+    "#In a beta distribution with a = 1, b = 2, there is a 75% chance of seeing less than 50% heads.\n",
+    "scipy.stats.beta(1, 2).cdf(.5) - scipy.stats.beta(1, 2).cdf(0)"
    ]
   },
   {
@@ -903,6 +904,13 @@
     "    ..."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "See `solutions/credible_interval.py` for solution"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -986,9 +994,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.9.6"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

notebooks/Module 2.3.2 - Confidence Intervals.ipynb

@@ -674,8 +674,8 @@
     "    \"\"\"Computer confidence interval for the given values\"\"\"\n",
     "    assert 0 < ci <= 1\n",
     "    n = len(values)\n",
-    "    lower = int(n * (1-ci)/2)\n",
-    "    upper = int(n * ci / 2)\n",
+    "    lower = int(n * (1-ci)/2) \n",
+    "    upper = int(n * (1-((1-ci)/2)))\n",
     "    assert upper > lower  # Can be lower == upper if not enough samples\n",
     "    sorted_values = np.sort(values)\n",
     "    return sorted_values[lower], sorted_values[upper]\n"
@@ -786,9 +786,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.9.6"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

notebooks/solutions/credible_interval.py

@@ -0,0 +1,54 @@
+def shortest_interval(posterior, alpha=0.05, sensitivity=.01):
+    """
+    A (1 - α) credible interval
+
+    Parameters
+    ----------
+    posterior : function
+        A posterior probability density function (normalized)
+
+    alpha : float
+        A significance level (of sorts)
+
+    Returns
+    -------
+        A tuple (lower, upper) giving the lower and upper bounds for
+        the (1 - α) credible interval
+    """
+
+    curr_prob = 0.0
+    target = 1 - alpha
+    min_length = 1
+    optimal_start = 0
+
+    # Initialize starting and ending indexes
+    start = 0
+    end = 0
+    while (end < 1):
+
+        # Increase upper bound while
+        # sum is smaller than or equal to 1-a
+        while (curr_prob <= target and end < 1):
+            end += sensitivity
+            curr_prob, err = integrate.quad(posterior, start, end)
+
+
+        # If current sum becomes greater than x
+        # start increasing lower bound, and save 
+        # smallest interval.
+        while (curr_prob > target and start < 1):
+
+            # new min length and new start
+            if min_length > (end - start):
+                min_length = end - start
+                optimal_start = start
+
+            start += sensitivity
+            curr_prob, err = integrate.quad(posterior, start, end)
+
+
+    return (optimal_start, optimal_start+min_length)
+
+beta = scipy.stats.beta(100,200)
+low, high = hpd(beta.pdf, alpha=.10, sensitivity=.005)
+beta.cdf(high) - beta.cdf(low) #this is equal to 1-a