<jupyter_start><jupyter_text>(notebook_name)= Chapter 3 The Beta-Binomial Bayesian Model[Chapter 3 The Beta-Binomial Bayesian Model](https://www.bayesrulesbook.com/chapter-3.html):::{post} March 24, 2022:tags: python:category: beginner:author: <NAME>:::<jupyter_code>import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pyreadr
import requests
import seaborn as sns

from scipy import integrate, stats
%config InlineBackend.figure_format = 'retina'
%load_ext watermark
RANDOM_SEED = 1301
np.random.seed(RANDOM_SEED)
az.style.use("arviz-darkgrid")<jupyter_output><empty_output><jupyter_text>[Chapter 3 The Beta-Binomial Bayesian Model](https://www.bayesrulesbook.com/chapter-3.html) [3.1.2 Tuning the Beta prior](https://www.bayesrulesbook.com/chapter-3.htmltuning-the-beta-prior)<jupyter_code>def plot_beta(alpha, beta, mean=False, mode=False):
    θ = np.linspace(0, 1, 100)
    p_theta_given_y = stats.beta.pdf(θ, alpha, beta)
    plt.plot(θ, p_theta_given_y, lw=4, color="k")
    plt.xlabel(r"$\pi$")
    plt.ylabel("$f(x)$")
# Plot the Beta(45, 55) prior
plot_beta(45, 55)
def summarize_beta(alpha, beta):
    r"""Summarize a Beta Model for \eqn{\pi}

    @param alpha,beta positive shape parameters of the Beta model

    Return Pandas Series with summary
    """

    mean = alpha / (alpha + beta)
    var = alpha * beta / ((alpha + beta) ** 2 * (alpha + beta + 1))
    sd = np.sqrt(var)

    if alpha < 1 and beta < 1:
        mode = "0 and 1"
    elif alpha <= 1 and beta > 1:
        mode = 0
    elif alpha > 1 and beta < 1:
        mode = 1
    else:
        mode = (alpha - 1) / (alpha + beta - 2)
    return pd.Series({"mean": mean, "mode": mode, "var": var, "sd": sd})
summarize_beta(45, 55)<jupyter_output><empty_output><jupyter_text>[3.3 The Beta posterior model](https://www.bayesrulesbook.com/chapter-3.htmlthe-beta-posterior-model)<jupyter_code>def plot_beta_binomial(
    alpha, beta, y=None, n=None, prior=True, likelihood=True, posterior=True
) -> None:
    """Plot a Beta-Binomial Bayesian Model

    @param alpha,beta positive shape parameters of the prior Beta model
    @param y observed number of successes
    @param n observed number of trials
    @param prior a logical value indicating whether the prior model should be plotted
    @param likelihood a logical value indicating whether the scaled likelihood should be plotted
    @param posterior a logical value indicating whether posterior model should be plotted
    """

    if y is None or n is None:
        print("Warning: to visualize the posterior specify function parameters y and n")

    θ = np.linspace(0, 1, 100)
    p_theta_given_y = stats.beta.pdf(θ, alpha, beta)
    plt.fill_between(θ, p_theta_given_y, lw=4, color="k", label="prior", alpha=0.2)
    alpha_post = alpha + y
    beta_post = beta + n - y
    p_theta_given_y_post = stats.beta.pdf(θ, alpha_post, beta_post)
    plt.fill_between(θ, p_theta_given_y_post, lw=4, color="r", label="posterior", alpha=0.2)
    likelihood = stats.binom.pmf(y, n, θ)
    scale_factor = integrate.simpson(likelihood, θ)
    plt.fill_between(
        θ, likelihood / scale_factor, lw=4, color="b", label="likelihood scaled", alpha=0.2
    )
    plt.xlabel(r"$\pi$")
    plt.ylabel("density")
    plt.legend()
# FIGURE 3.7: The prior pdf, scaled likelihood function, and posterior pdf of Michelle’s election support π.
plot_beta_binomial(alpha=45, beta=55, y=30, n=50)
def summarize_beta_binomial(alpha, beta, y=None, n=None):
    """Summarize a Beta-Binomial Bayesian model

    @param alpha,beta positive shape parameters of the prior Beta model
    @param y number of successes
    @param n number of trials

    Return: Pandas dataframe summarizing beta binomial
    """

    def beta_mean(a, b):
        return a / (a + b)

    def beta_mode(a, b):
        if a < 1 and b < 1:
            return "0 and 1"
        elif a <= 1 and b > 1:
            return 0
        elif a > 1 and b < 1:
            return 1
        else:
            return (a - 1) / (a + b - 2)

    def beta_var(a, b):
        return a * b / ((a + b) ** 2 * (a + b + 1))

    prior_mean = beta_mean(alpha, beta)
    prior_mode = beta_mode(alpha, beta)
    prior_var = beta_var(alpha, beta)
    prior_sd = np.sqrt(prior_var)
    if y is None and n is None:
        summary = pd.DataFrame(
            {
                "alpha": alpha,
                "beta": beta,
                "mean": prior_mean,
                "mode": prior_mode,
                "var": prior_var,
                "sd": prior_sd,
            },
            index=["prior"],
        )
    else:
        post_alpha = y + alpha
        post_beta = n - y + beta
        post_mean = beta_mean(post_alpha, post_beta)
        post_mode = beta_mode(post_alpha, post_beta)
        post_var = beta_var(post_alpha, post_beta)
        post_sd = np.sqrt(post_var)
        summary = pd.DataFrame(
            {
                "alpha": [alpha, post_alpha],
                "beta": [beta, post_beta],
                "mean": [prior_mean, post_mean],
                "mode": [prior_mode, post_mode],
                "var": [prior_var, post_var],
                "sd": [prior_sd, post_sd],
            },
            index=["prior", "posterior"],
        )
    return summary
summarize_beta_binomial(alpha=45, beta=55, y=30, n=50)<jupyter_output><empty_output><jupyter_text>[3.5 Simulating the Beta-Binomial](https://www.bayesrulesbook.com/chapter-3.htmlchapter-3-simulation)<jupyter_code>michelle_sim = pd.DataFrame({"pi": np.random.beta(a=45, b=55, size=10_000)})
michelle_sim["y"] = [np.random.binomial(n=50, p=pi) for pi in michelle_sim["pi"]]
michelle_sim.head()
ax = michelle_sim.plot.scatter(x="pi", y="y")
ax.axhline(y=30, color="w", linewidth=4)
ax.text(0.67, 30, "y == 30", fontsize=12, verticalalignment="center")
ax.set_title(
    "FIGURE 3.8: A scatterplot of 10,000 simulated pairs of Michelle’s support π and polling outcome y."
);
# Keep only the simulated pairs that match our data
michelle_posterior = michelle_sim.query("y == 30")
# Plot the remaining pi values
ax = michelle_posterior["pi"].plot.density()
ax.set_xlabel("pi")
ax.set_title(
    "FIGURE 3.9: A density plot of simulated π values that produced polling outcomes in which Y=30 voters supported Michelle."
);
michelle_posterior["pi"].describe()<jupyter_output><empty_output><jupyter_text>[3.6 Example: Milgram’s behavioral study of obedience](https://www.bayesrulesbook.com/chapter-3.htmlmilgram-3)<jupyter_code># Beta(1, 10) prior
plot_beta(alpha=1, beta=10);
# FIGURE 3.10: A Beta(1,10) prior model of π.
summarize_beta_binomial(alpha=1, beta=10, y=26, n=40)
# FIGURE 3.11: The Beta prior pdf, scaled Binomial likelihood function, and Beta posterior pdf for π,
# the proportion of subjects that would follow the given instructions.

plot_beta_binomial(alpha=1, beta=10, y=26, n=40)
%load_ext watermark
%watermark -n -u -v -iv -w<jupyter_output>Last updated: Sun Mar 27 2022

Python implementation: CPython
Python version       : 3.10.4
IPython version      : 8.1.1

Watermark: 2.3.0
{'fraction_non_alphanumeric': 0.1235543018335684, 'fraction_numerical': 0.026657263751763048, 'mean_word_length': 3.198342214328005, 'pattern_counts': {'":': 17, '<': 29, '<?xml version=': 0, '>': 26, 'https://': 6, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '471707', 'license': 'MIT', 'n_tokens_mistral': 2633, 'n_tokens_neox': 2384, 'n_words': 677, 'path': 'Bayes_Rules/Chapter_3_The_Beta-Binomial_Bayesian_Model.ipynb', 'repo_name': 'rserran/resources', 'size': '7101'}