import numpy as np


def beta_posterior(alpha_prior, beta_prior, heads, tails):
    return alpha_prior + heads, beta_prior + tails


def beta_mean(alpha, beta):
    return alpha / (alpha + beta)


def beta_credible_interval(alpha, beta, n_samples=200_000, level=0.95):
    samples = np.random.beta(alpha, beta, size=n_samples)
    lower = np.percentile(samples, 100 * (1 - level) / 2)
    upper = np.percentile(samples, 100 * (1 + level) / 2)
    return lower, upper


def sequential_update(alpha0, beta0, flips):
    """Update the posterior one flip at a time; final result must
    match the closed-form batch update, since Bayesian updating is
    order-independent for i.i.d. observations."""
    alpha, beta = alpha0, beta0
    history = [(alpha, beta)]
    for flip in flips:
        if flip == "H":
            alpha += 1
        else:
            beta += 1
        history.append((alpha, beta))
    return history


if __name__ == "__main__":
    np.random.seed(0)
    true_theta = 0.7
    flips = [
        "H" if np.random.rand() < true_theta else "T"
        for _ in range(50)
    ]
    heads, tails = flips.count("H"), flips.count("T")

    alpha0, beta0 = 1, 1  # uniform prior (Beta(1,1))
    alpha_n, beta_n = beta_posterior(alpha0, beta0, heads, tails)

    print(f"Observed: {heads} heads, {tails} tails")
    print(f"Posterior mean: {beta_mean(alpha_n, beta_n):.3f}")
    lo, hi = beta_credible_interval(alpha_n, beta_n)
    print(f"95% credible interval: [{lo:.3f}, {hi:.3f}]")

    # verify sequential updating matches the closed-form batch result
    history = sequential_update(alpha0, beta0, flips)
    assert history[-1] == (alpha_n, beta_n)
    print("Sequential update matches batch update:", history[-1])
