贝叶斯分析助你成为优秀的调参侠：先验和似然对后验的影响

from scipy.stats import beta    
import numpy as np
import matplotlib.pyplot as plt
# 需要安装 SciencePlots 库：
# pip install SciencePlots
plt.style.use(['science', 'ieee'])

def plot_beta_dist(axis, a=1, b=1):
    '''plot the Beta distribution function
    :a: parameter $\alpha$ in Beta distribution
    :b: parameter $\beta$ in Beta distribution'''
    x = np.linspace(0, 1, 10000)
    y = beta.pdf(x, a, b)
    axis.plot(x, y)
    axis.set_ylabel(r"$P(\theta)$")
    axis.set_xlabel(r"$\theta$")
    axis.set_title(r"$\alpha=%.1f, \beta=%.1f$"%(a, b))
# Let's investigate how does Beta distribution change with alpha and beta

# parameters in 2 rows, 3 columns
params = [[(1, 1), (0.6, 0.6), (10, 10)],
          [(0.3, 1.0), (10, 1), (1, 10)]]

fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(9, 6))
for i in range(2):
    for j in range(3):
        plot_beta_dist(ax[i][j], *params[i][j])

plt.tight_layout()

from ipywidgets import interact
import ipywidgets as widgets
print_to_pdf = False

def plot_posterior(m, n, a, b):
    theta = np.linspace(0.0001, 1, 100, endpoint=False)
    #poster = theta**(a + m - 1) * (1 - theta)**(b + n - m - 1)
    ln_poster = (a +m - 1) * np.log(theta) + (b + n - m - 1) * np.log(1 - theta)
    poster = np.exp(ln_poster)
    prior =  beta.pdf(theta, a, b)
    plt.axvline(m/n, color='c', label='frequency')
    plt.plot(theta, poster/poster.max(), label="posterior")
    plt.plot(theta, prior/prior.max(), 'k--', label="prior")
    plt.legend(loc='best', fontsize="small")
    plt.xlabel(r"$\theta$")

if print_to_pdf:
    plot_posterior(m=35, n=162, a=0.6, b=10)
else:
    interact(plot_posterior, 
         m=(10, 100), n=(20,200), 
         a = (0, 10.0, 0.2), b = (0, 10.0, 0.2))

def draw_posterior(m, n, a, b, nsamples=1000000):
    def poster(theta, m=m, n=n, a=a, b=b):
        ln_poster = (a +m - 1) * np.log(theta) \
                  + (b + n - m - 1) * np.log(1 - theta)
        poster = np.exp(ln_poster)
        return poster


    samples = mcmc(poster, np.array([[0.0001, 0.9999]]),
                   n=nsamples, warmup=0.1, thin=10)

    # compute the distribution of these samples using histogram
    hist, edges = np.histogram(samples.flatten(), bins=100)

    x = 0.5 * (edges[1:] + edges[:-1])
    theta = np.linspace(0.0001, 0.9999, 100)
    y = poster(theta)

    plt.plot(x, hist/hist.max(), '-', label="MCMC")
    plt.plot(theta, y/y.max(), '-', label = r"$P(\theta | \rm m\ out\ of\ n)$")

    plt.legend(loc='best', fontsize='small')
    plt.xlabel(r'$\theta$')
    plt.xlim(x[0], x[-1])

draw_posterior(m=100, n=150, a=10, b=10, nsamples=2000000)