diff --git a/notebooks/Module 1.5.2 - Bayesian Inference Basics.ipynb b/notebooks/Module 1.5.2 - Bayesian Inference Basics.ipynb index 43c43b4..070f0ac 100644 --- a/notebooks/Module 1.5.2 - Bayesian Inference Basics.ipynb +++ b/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 } diff --git a/notebooks/Module 2.3.2 - Confidence Intervals.ipynb b/notebooks/Module 2.3.2 - Confidence Intervals.ipynb index b511261..7f21bcf 100644 --- a/notebooks/Module 2.3.2 - Confidence Intervals.ipynb +++ b/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 } diff --git a/notebooks/solutions/credible_interval.py b/notebooks/solutions/credible_interval.py new file mode 100644 index 0000000..f1d199e --- /dev/null +++ b/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 \ No newline at end of file