Chapter 4: Introduction to Bayesian Statistics

Introduction

In the previous chapters, we learned the fundamentals of descriptive statistics and probability. In this chapter, as an application of those ideas, we will cover Bayesian Statistics, whose importance has been growing steadily in recent years. Bayesian statistics is a mathematical framework that formalizes something close to the human learning process: updating our knowledge each time new data is observed.

In Frequentist Statistics, which we have assumed so far, population parameters (for example, the probability that a coin lands heads) are treated as unknown but fixed values, and we estimate that value from observed data. Bayesian statistics, by contrast, treats the parameter itself as a probability distribution and updates that distribution every time new data is observed. This difference is not merely a matter of mathematical formalism — it represents an important shift in how uncertainty itself is understood.

1. A Deeper Look at Bayes' Theorem

In Chapter 1, we learned Bayes' theorem as a relationship between the probabilities of events (positive/negative, disease/healthy, and so on). In this chapter, we revisit it from the perspective of parameter estimation.

The basic idea behind Bayesian estimation is to update our belief about an unknown parameter $\theta$ given data $D$. Written in the context of parameter estimation, Bayes' theorem takes the following form.

$$P(\theta|D) = \frac{P(D|\theta)P(\theta)}{P(D)}$$

The name and role of each term are as follows.

• $P(\theta)$: the Prior Distribution — our belief about the parameter $\theta$ before observing any data
• $P(D|\theta)$: the Likelihood — the probability (density) of actually observing data $D$, given that the parameter is $\theta$
• $P(\theta|D)$: the Posterior Distribution — our updated belief about the parameter $\theta$ after observing data $D$
• $P(D)$: the Marginal Likelihood or Evidence — the probability (density) of observing data $D$ across all possible values of $\theta$. It acts as a normalizing constant

When $\theta$ takes continuous values, the marginal likelihood is given by the following integral.

$$P(D) = \int P(D|\theta)P(\theta)\,d\theta$$

Since this integral often cannot be solved analytically, it poses a significant challenge in practical Bayesian statistics. MCMC, which we cover later in this chapter, is one of the techniques developed to work around this challenge.

2. The Relationship Between Prior, Likelihood, and Posterior

The overall picture of Bayesian estimation can be organized as the interaction of the following three elements.

2.1 The Prior Distribution: Knowledge Before Observation

The prior distribution expresses knowledge or assumptions held before looking at the data. Past research, expert opinion, or even a state of "knowing nothing" can all be formalized as a prior distribution.

• Uninformative Prior: a prior distribution that treats every value as roughly equally likely, incorporating almost no prior knowledge
• Weakly Informative Prior: a prior distribution that only constrains a rough range of values
• Informative Prior: a prior distribution with strong constraints, based on past data or domain expertise

2.2 The Likelihood: What the Data Tell Us

The likelihood function models the process that generated the data. For example, the observation of $k$ heads out of $n$ coin flips can be modeled with a binomial distribution using the parameter $\theta$ (the probability of heads).

$$P(D|\theta) = \binom{n}{k}\theta^k(1-\theta)^{n-k}$$

2.3 The Posterior Distribution: Updated Knowledge

The posterior distribution is obtained by multiplying the prior distribution by the likelihood (and then normalizing). As the amount of data grows, the influence of the likelihood becomes relatively larger, while the influence of the prior distribution fades. Conversely, while data is still scarce, the prior distribution strongly shapes the form of the posterior distribution.

3. Using Conjugate Priors

A concept that is extremely convenient for carrying out Bayesian computations analytically is the Conjugate Prior.

Below are some representative combinations of conjugate priors.

3.1 Conjugacy Between the Beta Distribution and Bernoulli Trials

The Beta Distribution is a probability distribution for a random variable taking values in the interval $[0,1]$, making it well suited to expressing the uncertainty in a probability $\theta$ itself. Its probability density function is as follows.

$$f(\theta;\alpha,\beta) = \frac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{B(\alpha,\beta)}$$

Here $B(\alpha,\beta)$ is the beta function, which acts as a normalizing constant.

If we set the prior distribution to Beta($\alpha$, $\beta$) and observe $k$ successes out of $n$ Bernoulli Trials (trials with two outcomes, success/failure), the posterior distribution can be obtained analytically as follows.

$$\theta \sim \text{Beta}(\alpha, \beta) \quad \Rightarrow \quad \theta|D \sim \text{Beta}(\alpha+k, \beta+n-k)$$

As this update rule shows, the parameter $\alpha$ can be interpreted as a "hypothetical number of successes" and $\beta$ as a "hypothetical number of failures"; the posterior distribution is obtained simply by adding the observed number of successes and failures to each, respectively.

3.2 Python Implementation: Beta-Bernoulli Updating

Let's implement a repeated coin flip experiment and observe how the posterior distribution changes as more data is observed.

import numpy as np
from scipy import stats

# Generate coin flips using a true probability (set for simulation purposes, treated as unknown)
np.random.seed(42)
true_p = 0.7
n_flips = 50
flips = np.random.binomial(1, true_p, size=n_flips)  # 1=heads, 0=tails

# Prior distribution Beta(2, 2): leans toward fair but weakly informative
alpha_prior, beta_prior = 2, 2

# Update the posterior parameters as the number of observations increases
checkpoints = [0, 5, 20, 50]

print("=== Progression of Beta-Bernoulli Updating ===")
for n in checkpoints:
    successes = flips[:n].sum()
    failures = n - successes
    alpha_post = alpha_prior + successes
    beta_post = beta_prior + failures
    posterior_mean = alpha_post / (alpha_post + beta_post)
    print(f"After n={n:2d} observations: {successes} successes, "
          f"Beta({alpha_post}, {beta_post}), posterior mean={posterior_mean:.3f}")

Execution Result:

=== Progression of Beta-Bernoulli Updating ===
After n= 0 observations: 0 successes, Beta(2, 2), posterior mean=0.500
After n= 5 observations: 3 successes, Beta(5, 4), posterior mean=0.556
After n=20 observations: 14 successes, Beta(16, 8), posterior mean=0.667
After n=50 observations: 39 successes, Beta(41, 13), posterior mean=0.759

Let's visualize how the posterior distribution grows sharper as more observations are added to the prior.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

np.random.seed(42)
true_p = 0.7
n_flips = 50
flips = np.random.binomial(1, true_p, size=n_flips)

alpha_prior, beta_prior = 2, 2
checkpoints = [0, 5, 20, 50]
theta = np.linspace(0.001, 0.999, 200)

plt.figure(figsize=(10, 6))
for n in checkpoints:
    successes = flips[:n].sum()
    failures = n - successes
    alpha_post = alpha_prior + successes
    beta_post = beta_prior + failures
    pdf = stats.beta.pdf(theta, alpha_post, beta_post)
    plt.plot(theta, pdf, label=f'n={n} ({successes} successes)', linewidth=2)

plt.axvline(true_p, color='black', linestyle=':', label=f'True probability={true_p}')
plt.xlabel('θ (probability of heads)', fontsize=12)
plt.ylabel('Probability density', fontsize=12)
plt.title('Change in the Posterior Distribution as Observations Increase', fontsize=14)
plt.legend()
plt.grid(alpha=0.3)
plt.show()

4. Implementing Bayesian Estimation

Once we have a posterior distribution, we need to extract concrete estimates and measures of uncertainty from it.

4.1 Point Estimation: Posterior Mean and MAP Estimation

• Posterior Mean: the expected value $E[\theta|D]$ of the posterior distribution. The most commonly used point estimator
• MAP Estimation (Maximum A Posteriori estimation): the value of $\theta$ at which the posterior density is maximized, i.e., the mode of the posterior distribution

For the Beta($\alpha$, $\beta$) distribution, each can be obtained analytically using the following formulas.

$$E[\theta|D] = \frac{\alpha}{\alpha+\beta} \qquad \theta_{\text{MAP}} = \frac{\alpha-1}{\alpha+\beta-2}$$

(Note that the MAP estimation formula holds only when $\alpha > 1$ and $\beta > 1$)

4.2 Interval Estimation: The Credible Interval

A Credible Interval is an interval of the posterior distribution that contains a specified probability mass (for example, 95%). It plays a role very similar to the Confidence Interval in frequentist statistics, but its interpretation is fundamentally different (we examine this difference in detail in Section 6).

4.3 Python Implementation: Point and Interval Estimation

import numpy as np
from scipy import stats

def posterior_summary(alpha_post, beta_post, cred_mass=0.95):
    """
    Compute point estimates and a credible interval from a beta posterior distribution

    Parameters:
    -----------
    alpha_post, beta_post : float
        Parameters of the posterior distribution Beta(alpha_post, beta_post)
    cred_mass : float
        Probability mass of the credible interval (default 95%)

    Returns:
    --------
    mean, map_estimate, (lower, upper)
    """
    mean = alpha_post / (alpha_post + beta_post)
    if alpha_post > 1 and beta_post > 1:
        map_estimate = (alpha_post - 1) / (alpha_post + beta_post - 2)
    else:
        map_estimate = None
    lower = stats.beta.ppf((1 - cred_mass) / 2, alpha_post, beta_post)
    upper = stats.beta.ppf(1 - (1 - cred_mass) / 2, alpha_post, beta_post)
    return mean, map_estimate, (lower, upper)

# The coin-flip data from the previous section (39 successes out of 50, prior Beta(2,2))
alpha_post, beta_post = 41, 13
mean, map_est, (lower, upper) = posterior_summary(alpha_post, beta_post)

print(f"Posterior mean: {mean:.4f}")
print(f"MAP estimate: {map_est:.4f}")
print(f"95% credible interval: [{lower:.4f}, {upper:.4f}]")

Execution Result:

Posterior mean: 0.7593
MAP estimate: 0.7692
95% credible interval: [0.6379, 0.8624]

4.4 Practical Example: A Bayesian Comparison for A/B Testing

Let's put Bayesian estimation into practice using an A/B test comparing the click-through rates of two website designs (Design A and Design B).

import numpy as np
from scipy import stats

# Design A: 48 clicks out of 500 visitors, Design B: 63 clicks out of 480 visitors
alpha_prior, beta_prior = 1, 1  # Uninformative prior Beta(1,1) = uniform distribution

n_a, success_a = 500, 48
n_b, success_b = 480, 63

alpha_a = alpha_prior + success_a
beta_a = beta_prior + (n_a - success_a)
alpha_b = alpha_prior + success_b
beta_b = beta_prior + (n_b - success_b)

def posterior_summary(alpha_post, beta_post, cred_mass=0.95):
    mean = alpha_post / (alpha_post + beta_post)
    map_estimate = (alpha_post - 1) / (alpha_post + beta_post - 2)
    lower = stats.beta.ppf((1 - cred_mass) / 2, alpha_post, beta_post)
    upper = stats.beta.ppf(1 - (1 - cred_mass) / 2, alpha_post, beta_post)
    return mean, map_estimate, (lower, upper)

mean_a, map_a, ci_a = posterior_summary(alpha_a, beta_a)
mean_b, map_b, ci_b = posterior_summary(alpha_b, beta_b)

print(f"Design A: posterior mean={mean_a:.4f}, 95% credible interval=[{ci_a[0]:.4f}, {ci_a[1]:.4f}]")
print(f"Design B: posterior mean={mean_b:.4f}, 95% credible interval=[{ci_b[0]:.4f}, {ci_b[1]:.4f}]")

# Estimate P(Design B is better) using Monte Carlo simulation
np.random.seed(0)
n_mc = 200000
samples_a = np.random.beta(alpha_a, beta_a, n_mc)
samples_b = np.random.beta(alpha_b, beta_b, n_mc)
prob_b_better = np.mean(samples_b > samples_a)

print(f"\nP(Design B's click rate > Design A's click rate) = {prob_b_better:.4f}")

Execution Result:

Design A: posterior mean=0.0976, 95% credible interval=[0.0732, 0.1250]
Design B: posterior mean=0.1328, 95% credible interval=[0.1040, 0.1644]

P(Design B's click rate > Design A's click rate) = 0.9582

5. An Introduction to Markov Chain Monte Carlo (MCMC)

When a conjugate prior is available, as with the beta distribution and Bernoulli trials, the posterior distribution can be obtained analytically. However, in many real-world Bayesian models, the combination of likelihood and prior is not conjugate, and the posterior distribution cannot be written down in closed form. This is exactly the situation where Markov Chain Monte Carlo (MCMC) proves powerful.

MCMC refers to a family of techniques that generate a sequence of samples following the posterior distribution without ever computing the normalizing constant $P(D)$. With the resulting set of samples, quantities such as the posterior mean and credible intervals can be approximated using histograms or sample averages.

5.1 Markov Chains and the Monte Carlo Method

The Monte Carlo Method is a general term for techniques that use repeated random trials to approximately compute quantities (such as integrals or expected values) that are difficult to obtain analytically.

A Markov Chain is a stochastic process with the property (the Markov property) that the next state depends only on the current state and not on any earlier history. MCMC constructs a Markov chain whose Stationary Distribution is the target posterior distribution, and generates samples by performing a random walk along that chain.

5.2 The Metropolis Algorithm

Among MCMC methods, the Metropolis Algorithm is one of the most intuitive to understand. The procedure is as follows.

1. Set an arbitrary initial value $\theta_0$
2. From the current state $\theta_t$, generate a new candidate $\theta^*$ using a Proposal Distribution (often a symmetric distribution such as a normal distribution)
3. Compute the Acceptance Probability $a = \min\left(1, \dfrac{P(\theta^*|D)}{P(\theta_t|D)}\right)$
4. Accept the candidate with probability $a$, setting $\theta_{t+1}=\theta^*$; otherwise reject it and set $\theta_{t+1}=\theta_t$
5. Repeat steps 2–4 as many times as needed

The early portion of the generated sample sequence is strongly influenced by the initial value $\theta_0$ and may not yet have converged to the stationary distribution. This early portion is typically discarded as the Burn-in period, with only the remaining samples used as an approximation of the posterior distribution.

5.3 Python Implementation: A Hand-Written Metropolis Sampler

Because the beta-Bernoulli model has an analytically tractable posterior distribution, we can compare the results obtained via MCMC against the analytical solution. This provides a good way to check that an MCMC implementation is working correctly.

import numpy as np
from scipy import stats

def log_posterior(theta_val, successes, failures, alpha_prior, beta_prior):
    """Compute the log of the unnormalized posterior distribution (returns -inf outside 0<theta<1)"""
    if theta_val <= 0 or theta_val >= 1:
        return -np.inf
    log_lik = successes * np.log(theta_val) + failures * np.log(1 - theta_val)
    log_prior = (alpha_prior - 1) * np.log(theta_val) + (beta_prior - 1) * np.log(1 - theta_val)
    return log_lik + log_prior

def metropolis_sampler(log_post_fn, n_samples, init, proposal_std, seed=42):
    """
    An MCMC sampler based on the Metropolis algorithm

    Parameters:
    -----------
    log_post_fn : callable
        A function returning the log posterior distribution (normalization not required)
    n_samples : int
        Number of samples to generate
    init : float
        Initial value
    proposal_std : float
        Standard deviation of the proposal distribution (normal distribution)

    Returns:
    --------
    samples : ndarray
        The generated sequence of samples
    acceptance_rate : float
        The acceptance rate
    """
    rng = np.random.default_rng(seed)
    samples = np.zeros(n_samples)
    current = init
    current_log_p = log_post_fn(current)
    n_accepted = 0

    for i in range(n_samples):
        proposal = current + rng.normal(0, proposal_std)
        proposal_log_p = log_post_fn(proposal)
        log_accept_ratio = proposal_log_p - current_log_p

        if np.log(rng.uniform()) < log_accept_ratio:
            current = proposal
            current_log_p = proposal_log_p
            n_accepted += 1
        samples[i] = current

    acceptance_rate = n_accepted / n_samples
    return samples, acceptance_rate

# Use the same coin-flip data as in Section 3 (39 successes out of 50), with prior Beta(2,2)
successes_total, failures_total = 39, 11
alpha_prior, beta_prior = 2, 2
target = lambda th: log_posterior(th, successes_total, failures_total, alpha_prior, beta_prior)

n_samples = 20000
samples, acc_rate = metropolis_sampler(target, n_samples, init=0.5, proposal_std=0.1)

burn_in = 2000
post_burn = samples[burn_in:]

analytic_mean = (alpha_prior + successes_total) / (alpha_prior + beta_prior + successes_total + failures_total)

print(f"Acceptance rate: {acc_rate:.3f}")
print(f"Posterior mean from MCMC: {post_burn.mean():.4f}")
print(f"Posterior mean from the analytical solution (beta distribution): {analytic_mean:.4f}")
print(f"Posterior standard deviation from MCMC: {post_burn.std():.4f}")

Execution Result:

Acceptance rate: 0.548
Posterior mean from MCMC: 0.7601
Posterior mean from the analytical solution (beta distribution): 0.7593
Posterior standard deviation from MCMC: 0.0575

Let's visualize the behavior of the MCMC samples and compare them against the analytical solution.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def log_posterior(theta_val, successes, failures, alpha_prior, beta_prior):
    if theta_val <= 0 or theta_val >= 1:
        return -np.inf
    log_lik = successes * np.log(theta_val) + failures * np.log(1 - theta_val)
    log_prior = (alpha_prior - 1) * np.log(theta_val) + (beta_prior - 1) * np.log(1 - theta_val)
    return log_lik + log_prior

def metropolis_sampler(log_post_fn, n_samples, init, proposal_std, seed=42):
    rng = np.random.default_rng(seed)
    samples = np.zeros(n_samples)
    current = init
    current_log_p = log_post_fn(current)
    n_accepted = 0
    for i in range(n_samples):
        proposal = current + rng.normal(0, proposal_std)
        proposal_log_p = log_post_fn(proposal)
        if np.log(rng.uniform()) < proposal_log_p - current_log_p:
            current = proposal
            current_log_p = proposal_log_p
            n_accepted += 1
        samples[i] = current
    return samples, n_accepted / n_samples

successes_total, failures_total = 39, 11
alpha_prior, beta_prior = 2, 2
target = lambda th: log_posterior(th, successes_total, failures_total, alpha_prior, beta_prior)
samples, acc_rate = metropolis_sampler(target, 20000, init=0.5, proposal_std=0.1)
burn_in = 2000
post_burn = samples[burn_in:]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Trace plot: how the samples evolve over iterations
axes[0].plot(samples, linewidth=0.5, alpha=0.7)
axes[0].axvline(burn_in, color='red', linestyle='--', label='Burn-in boundary')
axes[0].set_xlabel('Iteration')
axes[0].set_ylabel('θ')
axes[0].set_title('Trace Plot')
axes[0].legend()

# Histogram vs. analytical solution
theta_grid = np.linspace(0.001, 0.999, 200)
analytic_pdf = stats.beta.pdf(theta_grid, alpha_prior + successes_total, beta_prior + failures_total)
axes[1].hist(post_burn, bins=50, density=True, alpha=0.6, label='MCMC samples')
axes[1].plot(theta_grid, analytic_pdf, 'r-', linewidth=2, label='Analytical solution Beta(41, 13)')
axes[1].set_xlabel('θ')
axes[1].set_ylabel('Probability density')
axes[1].set_title('MCMC Samples vs. Analytical Solution')
axes[1].legend()

plt.tight_layout()
plt.show()

6. Comparing Bayesian and Frequentist Statistics

Let's now organize the differences between the Bayesian statistics we've covered so far and the frequentist statistics we assumed up until this chapter. Both approaches draw inferences from the same data using different underlying philosophies, so it is important to understand the strengths and limitations of each.

6.1 Different Interpretations of Confidence Intervals and Credible Intervals

One of the most commonly confused points is the difference in interpretation between the frequentist Confidence Interval and the Bayesian Credible Interval.

• 95% Confidence Interval: a property of the overall procedure — "if the same data-collection and interval-computation procedure is repeated 100 times, on average 95 of the resulting intervals will contain the true parameter value." Whether any single interval contains the true value is either 0 or 1, and no probabilistic statement can be made about it
• 95% Credible Interval: a direct probabilistic statement about the parameter itself — "given the observed data, the probability that the parameter falls within this interval is 95%"

6.2 Differences in Small-Sample Behavior: An Implementation Comparison

The difference between the two approaches can become more pronounced when the sample size is small. Let's compare the frequentist Wald confidence interval and the Bayesian credible interval using a small dataset of 8 successes out of 10 trials.

import numpy as np
from scipy import stats

n = 10
k = 8

# Frequentist: maximum likelihood estimate and Wald confidence interval
p_hat = k / n  # Maximum Likelihood Estimate (MLE)
se = np.sqrt(p_hat * (1 - p_hat) / n)
z = 1.96
freq_ci = (p_hat - z * se, p_hat + z * se)

# Bayesian: posterior distribution using an uninformative prior Beta(1,1)
alpha_prior, beta_prior = 1, 1
alpha_post = alpha_prior + k
beta_post = beta_prior + (n - k)
bayes_mean = alpha_post / (alpha_post + beta_post)
bayes_ci = stats.beta.ppf([0.025, 0.975], alpha_post, beta_post)

print(f"Maximum likelihood estimate: {p_hat:.4f}")
print(f"Frequentist 95% confidence interval (Wald method): [{freq_ci[0]:.4f}, {freq_ci[1]:.4f}]")
print(f"\nBayesian posterior mean: {bayes_mean:.4f}")
print(f"Bayesian 95% credible interval: [{bayes_ci[0]:.4f}, {bayes_ci[1]:.4f}]")

Execution Result:

Maximum likelihood estimate: 0.8000
Frequentist 95% confidence interval (Wald method): [0.5521, 1.0479]

Bayesian posterior mean: 0.7500
Bayesian 95% credible interval: [0.4822, 0.9398]

This example is not about which method is "absolutely correct" — rather, it shows that each method operates under different assumptions and approximations. As the sample size grows and a weak prior is used, the results from frequentist and Bayesian methods generally converge to similar values. It is important to develop the judgment needed to decide which framework is appropriate for a given situation.

7. Summary and Next Steps

In this chapter, we studied the ideas behind Bayesian statistics consistently, from theory through to implementation in Python.

Reviewing the Learning Objectives

Let's revisit the learning objectives presented at the beginning of this chapter.

• ✅ Understand the way of thinking behind Bayesian statistics — organized as the relationship between the prior distribution, the likelihood, and the posterior distribution
• ✅ Explain the relationship between the prior distribution, the likelihood, and the posterior distribution — learned through equations and concrete examples in Sections 1 and 2
• ✅ Implement Bayesian estimation using conjugate priors — implemented the beta-Bernoulli model in Sections 3 and 4
• ✅ Understand the basic ideas behind MCMC — implemented a hand-written Metropolis sampler in Section 5 and checked it against the analytical solution
• ✅ Explain the differences between Bayesian statistics and frequentist statistics — confirmed the difference in interpretation between confidence intervals and credible intervals through concrete examples in Section 6

Next Steps

Building on the estimation framework covered in this chapter, the next chapter moves on to further topics in statistical inference, such as hypothesis testing. The ideas behind Bayesian statistics also form the foundation for more advanced techniques that appear frequently in machine learning, such as Bayesian optimization and Bayesian neural networks.

Practice Problems

from scipy import stats

alpha_prior, beta_prior = 2, 2
n, k = 10, 7

alpha_post = alpha_prior + k
beta_post = beta_prior + (n - k)

mean = alpha_post / (alpha_post + beta_post)
map_estimate = (alpha_post - 1) / (alpha_post + beta_post - 2)

print(f"Posterior distribution: Beta({alpha_post}, {beta_post})")
print(f"Posterior mean: {mean:.4f}")
print(f"MAP estimate: {map_estimate:.4f}")

# Output:
# Posterior distribution: Beta(9, 5)
# Posterior mean: 0.6429
# MAP estimate: 0.6667

import numpy as np
from scipy import stats

n, k = 10, 7
alpha_post, beta_post = 9, 5

# Bayesian credible interval
bayes_ci = stats.beta.ppf([0.025, 0.975], alpha_post, beta_post)

# Frequentist Wald confidence interval
p_hat = k / n
se = np.sqrt(p_hat * (1 - p_hat) / n)
freq_ci = (p_hat - 1.96 * se, p_hat + 1.96 * se)

print(f"Bayesian 95% credible interval: [{bayes_ci[0]:.4f}, {bayes_ci[1]:.4f}]")
print(f"Frequentist 95% confidence interval (Wald method): [{freq_ci[0]:.4f}, {freq_ci[1]:.4f}]")

# Output:
# Bayesian 95% credible interval: [0.3857, 0.8614]
# Frequentist 95% confidence interval (Wald method): [0.4160, 0.9840]

import numpy as np

def log_posterior(theta_val, successes, failures, alpha_prior, beta_prior):
    if theta_val <= 0 or theta_val >= 1:
        return -np.inf
    log_lik = successes * np.log(theta_val) + failures * np.log(1 - theta_val)
    log_prior = (alpha_prior - 1) * np.log(theta_val) + (beta_prior - 1) * np.log(1 - theta_val)
    return log_lik + log_prior

def metropolis_sampler(log_post_fn, n_samples, init, proposal_std, seed=42):
    rng = np.random.default_rng(seed)
    samples = np.zeros(n_samples)
    current = init
    current_log_p = log_post_fn(current)
    n_accepted = 0
    for i in range(n_samples):
        proposal = current + rng.normal(0, proposal_std)
        proposal_log_p = log_post_fn(proposal)
        if np.log(rng.uniform()) < proposal_log_p - current_log_p:
            current = proposal
            current_log_p = proposal_log_p
            n_accepted += 1
        samples[i] = current
    return samples, n_accepted / n_samples

target = lambda th: log_posterior(th, successes=7, failures=3, alpha_prior=2, beta_prior=2)

for std in [0.01, 1.0]:
    samples, acc_rate = metropolis_sampler(target, 5000, init=0.5, proposal_std=std)
    print(f"proposal_std={std}: acceptance rate={acc_rate:.3f}, sample standard deviation={samples[500:].std():.4f}")

# Output:
# proposal_std=0.01: acceptance rate=0.975, sample standard deviation=0.1299
# proposal_std=1.0: acceptance rate=0.156, sample standard deviation=0.1208