1. 频率派和贝叶斯派的本质区别

可以这么理解：贝叶斯学派发展的模型（如概率图模型）核心是计算后验分布，这通常涉及积分问题。当解析解不可得时，会采用数值积分方法，例如蒙特卡洛采样（如MCMC）来近似积分。而频率派发展的模型（通常称为统计机器学习）本质上是优化算法，步骤包括：第一步设计模型，第二步定义损失函数，第三步通过优化算法（例如梯度下降或牛顿法）最小化损失函数，从而将问题转化为求解最优解。

【参考】：b站-BV1cW411C7RS

这是大神shuhuai008在讲解贝叶斯派和频率派区别时，我认为的精彩之处。相应的，我有几个问题需要解决：
Q1：贝叶斯派经典模型有哪些？概率图模型是一个统称还是个例？
Q2：简单了解，蒙特卡洛方法是啥，干啥的？MCMC？
Q3：频率派经典模型有哪些？

这些问题稍后解决，不重要的问题解答，我都会放到文末，首先看主要内容。

---

对概率的诠释有两大学派，一种是概率派，另一种是贝叶斯派。我们假设有这样一组数据（观测集），共有N个样本，每个样本p个维度：
$$\left(\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1p}\\ x_{21}& x_{22}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{np}\end{array}\right)$$
假设这是一个概率模型，每个观测都是由$p(x|\theta )$生成的，即$x\stackrel{\text{i.i.d.}}{\sim }p(x|\theta )$。

2. 频率派的观点

频率派认为，模型的参数$\theta$是一个常量，虽然是未知的，数据$x$是随机变量。频率派统计的目标是基于随机样本对固定参数$\theta$进行推断，常用的方法就是极大似然估计，为了简化计算加上$log$运算：
$$\begin{array}{rl}{\theta }_{MLE}& =\underset{\theta }{\mathrm{argmax}}\log p(X|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}\log \prod _{i=1}^{N}p(x_i|\theta )\\ & =\underset{\theta }{\mathrm{argmax}}\sum _{i=1}^N\log p(x_i|\theta )\end{array}$$

学过极大似然估计，不过忘得一干二净，我想我需要简单了解一下。

Q4：简单了解下极大似然估计。
解：在抽样之后，我们已经得到了一组确定的观测值，这组数据肯定服从某个分布，具有某些确定的参数，但是，我们仍然把它当做一个更宏大的随机过程的一个结果来看待，所以$X$是随机变量。

极大似然估计，听起来高大上，拆开来看：
似然：在某个参数$\theta$下，观测到当前这组数据$X$的可能性有多大。
极大：最大化。
估计：找出这个参数$\theta$的值。
整理一下就是：找到一个参数$\theta$的值，使得在当前参数下，观察到这组数据的可能性达到最大。

一个简单的比喻：你有个朋友要打靶，但是你不知道他是高手还是菜鸟，并且你看不到靶子，只能根据他的成绩推断他是高手还是菜鸟。成绩出来了，他连续10枪命中10环。

那么，在“他是高手”和“他是菜鸟”两种参数设定下，连续10枪命中10环这件事发生的可能性哪个更大？你肯定会估计他是个高手，因为只有他是高手才可能打出这样的好成绩。当然，极大似然估计也是这么推断的。

下面分解公式：
${\theta }_{MLE}$就是最终要求的参数，即极大似然估计（Maximum Likelihood Estimation, MLE）出的参数值。
$\underset{\theta }{\mathrm{argmax}}$的意思是，让后面这个函数$logp(X|\theta )$取得最大值时，所对应的$\theta$的取值。
$p(X|\theta )$是似然函数，意思是：在参数为$\theta$时，观察到整个数据集$X$的概率密度。加上$log$为了简化计算。

我对自己的要求是，既要学会推导数学公式，也要用代码复现。即使在初学阶段，我想也要整个代码加深理解。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize # 优化算法库，用于寻找函数的最小值
from scipy.stats import norm    # 统计函数库，这里导入了正态分布
import seaborn as sns

# 在导入 matplotlib 之后，绘图之前添加以下代码
"""
解决绘图不能正常显示中文的问题
"""
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False    # 正确显示负号

# 设置随机种子以确保结果可重现
np.random.seed(42)

# 1. 生成模拟数据
def generate_data(n_samples=1000, true_mu=5, true_sigma=2):
    """
    生成服从正态分布的模拟数据
    """
    data = np.random.normal(true_mu, true_sigma, n_samples)
    """
    true_mu是分布的均值，true_sigma是标准差，n_samples是样本数量
    """
    return data

# 2. 定义似然函数和对数似然函数
def log_likelihood_normal(params, data):
    """
    正态分布的对数似然函数
    params: [mu, sigma]
    data: 观测数据
    """
    mu, sigma = params
    if sigma <= 0:  # 标准差必须为正，否则返回负无穷
        return -np.inf
    n = len(data)
    # 正态分布的对数似然函数的公式
    log_likelihood = -n/2 * np.log(2*np.pi) - n*np.log(sigma) - 1/(2*sigma**2)*np.sum((data - mu)**2)
    return log_likelihood

def negative_log_likelihood(params, data):
    """
    负对数似然函数（用于最小化）
    """
    return -log_likelihood_normal(params, data)
    

这部分代码很简单，主要是生成了1000个服从正态分布的模拟数据，然后定义了似然函数。有2个点需要明确说明一下：

第一，第36行代码，在定义log_likelihood_normal函数的这部分：

    if sigma <= 0:  # 标准差必须为正，否则返回负无穷
        return -np.inf

Q5：为什么要规定，标准差为负时，返回负无穷？
解：在概率分布中，标准差$\sigma$衡量的是数据的离散程度，显然$\sigma >0$。当然，$\sigma =0$意味着所有数据点的值完全没有波动，完全相同，这个分布退化成了一个点，虽然失去了分布的意义，但也是可以接受的。然而，$\sigma <0$没有任何数学意义了，距离和离散程度不可能是负数。所以，在程序中对$\sigma <0$的情况加以处理是自然而然的，那么问题进而转化为：为什么要返回负无穷？

因为此程序中，我们借助minimize函数来寻找最小值，但我们的目标是最大化似然函数，所以定义了负对数似然函数。最小化负对数似然函数就等价于最大化对数似然函数。当出现$\sigma <0$的情况后，返回负无穷，后续再设置阈值会排除掉这种极端情况。

为什么多此一举转化为寻找最小值，不直接寻找最大值？我查了一下资料，无非是更有利，于此不需要了解。

第二，第40行代码，正态分布的对数似然函数的公式是如何推导的：

    # 正态分布的对数似然函数的公式
    log_likelihood = -n/2 * np.log(2*np.pi) - n*np.log(sigma) - 1/(2*sigma**2)*np.sum((data - mu)**2)

回忆一下正态分布的概率密度函数（PDF）：
$$p(x|\mu ,\sigma )=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)$$
对于所有数据点独立同分布（$i.i.d.$）的情况，似然函数是每个数据点概率的乘积：
$$L(\mu ,\sigma |X)=\prod _{i=1}^{n}p(x_i|\mu ,\sigma )$$
取对数后，乘积就变成了求和：
$$\log L(\mu ,\sigma |X)=\sum _{i=1}^n\log p(x_i|\mu ,\sigma )$$
将正态分布的PDF带入上式，经过简单推导（我手推了一下，很简单，看着复杂而已），就可以得到程序中计算对数似然函数的公式：
$$\log L=-\frac{n}{2}\log (2\pi )-n\log (\sigma )-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(x_i-\mu {)}^2$$

---

继续展示代码，下面是核心内容：使用优化算法进行估计。

# 3. 使用优化算法进行极大似然估计
def mle_normal(data, initial_params=None):
    """
    使用MLE估计正态分布的参数
    """
    if initial_params is None:
        initial_params = [np.mean(data), np.std(data)]
    
    # 最小化负对数似然函数
    result = minimize(
        negative_log_likelihood,    # 需要最小化的函数
        initial_params,     # 初始化参数值
        args=(data,),
        bounds=[(None, None), (1e-6, None)],  # mu无限制，sigma>0
        method='L-BFGS-B'   # 优化方法
    )
    
    if result.success:
        return result.x
    else:
        raise ValueError("优化失败: " + result.message)

其中，bounds参数规定取值范围：$\mu$可以是任意实数，$\sigma$必须大于一个很小的正数（1e-6），也可以是1e-8或1e-10，但1e-6是一个安全和常用的起点，这个值已经能接近“物理意义”上的0。不能直接取0，似然函数在$\sigma =0$处未定义。
L-BFGS-B是众多优化算法中的一种，适用于有边界约束（比如$\sigma >0$）的优化问题，作用是找到负对数似然函数最小的参数值。

将真实参数和估计参数进行可视化对比：

# 4. 可视化函数
def plot_results(data, true_params, estimated_params):
    """
    可视化真实分布、估计分布和数据直方图
    """
    true_mu, true_sigma = true_params
    est_mu, est_sigma = estimated_params
    
    # 创建图形
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
    
    # 左图：数据直方图和分布曲线
    sns.histplot(data, bins=30, stat='density', alpha=0.7, ax=ax1, label='数据分布')
    
    # 生成x值范围
    x_min, x_max = min(data), max(data)
    x = np.linspace(x_min, x_max, 1000)
    
    # 绘制真实分布
    true_pdf = norm.pdf(x, true_mu, true_sigma)
    ax1.plot(x, true_pdf, 'r-', linewidth=2, label=f'真实分布: N({true_mu:.2f}, {true_sigma:.2f}$^2$)')
    
    # 绘制估计分布
    est_pdf = norm.pdf(x, est_mu, est_sigma)
    ax1.plot(x, est_pdf, 'g--', linewidth=2, label=f'MLE估计: N({est_mu:.2f}, {est_sigma:.2f}$^2$)')
    
    ax1.set_title('数据分布与参数估计')
    ax1.set_xlabel('值')
    ax1.set_ylabel('密度')
    ax1.legend()
    
    # 右图：似然函数曲面（简化版，固定一个参数）
    sigma_vals = np.linspace(0.5, 3.5, 100)
    log_likelihoods = [log_likelihood_normal([est_mu, s], data) for s in sigma_vals]
    
    ax2.plot(sigma_vals, log_likelihoods, 'b-')
    ax2.axvline(x=est_sigma, color='g', linestyle='--', label=f'MLE估计 σ={est_sigma:.2f}')
    ax2.axvline(x=true_sigma, color='r', linestyle='--', label=f'真实 σ={true_sigma:.2f}')
    ax2.set_title('对数似然函数随σ的变化 (μ固定)')
    ax2.set_xlabel('σ')
    ax2.set_ylabel('对数似然')
    ax2.legend()
    
    plt.tight_layout()
    plt.show()

图1

简单看一下，左图展示了数据的分布情况，一眼正态分布，毕竟我们是生成服从正态分布的模拟数据。真实分布和MLE估计的分布是非常接近的。右图展示了当固定$\mu$时（样本均值），对数似然函数随$\sigma$变化的情况：

• 负对数似然函数是一个凸函数，有明确的最大值点。正态分布的对数似然函数是凹的，相比于凹函数，许多优化算法在求解凸函数时效率最高，能保证找到全局最优解，而且在计算机计算时，可以防止数据溢出问题。这在另一种角度上，解释了Q5中，为什么当$\sigma <0$时，需要返回负无穷。
• MLE估计的$\sigma$值（1.9575）正好位于曲线的最大值点，与真实$\sigma$值（2.0）非常接近。

最后一部分代码，设置主程序，计算误差。

# 5. 主程序
def main():
    # 设置真实参数
    true_mu = 5.0
    true_sigma = 2.0
    n_samples = 1000
    
    print("开始生成模拟数据...")
    print(f"真实参数: μ = {true_mu}, σ = {true_sigma}")
    print(f"样本数量: {n_samples}")
    
    # 生成数据
    data = generate_data(n_samples, true_mu, true_sigma)
    
    print("\n计算样本统计量...")
    sample_mean = np.mean(data)
    sample_std = np.std(data)
    print(f"样本均值: {sample_mean:.4f}")
    print(f"样本标准差: {sample_std:.4f}")
    
    print("\n使用MLE估计参数...")
    # 使用样本统计量作为初始值
    initial_params = [sample_mean, sample_std]
    estimated_params = mle_normal(data, initial_params)
    
    print(f"MLE估计结果: μ = {estimated_params[0]:.4f}, σ = {estimated_params[1]:.4f}")
    print(f"与真实参数的差异: Δμ = {abs(estimated_params[0]-true_mu):.4f}, Δσ = {abs(estimated_params[1]-true_sigma):.4f}")
    
    # 可视化结果
    print("\n生成可视化图表...")
    plot_results(data, (true_mu, true_sigma), estimated_params)
    
    # 计算置信区间（使用Fisher信息矩阵近似）
    n = len(data)
    mu_se = estimated_params[1] / np.sqrt(n)  # 均值标准误
    sigma_se = estimated_params[1] / np.sqrt(2 * n)  # 标准差标准误
    
    mu_ci = (estimated_params[0] - 1.96 * mu_se, estimated_params[0] + 1.96 * mu_se)
    sigma_ci = (estimated_params[1] - 1.96 * sigma_se, estimated_params[1] + 1.96 * sigma_se)
    
    print(f"\n95%置信区间:")
    print(f"μ: [{mu_ci[0]:.4f}, {mu_ci[1]:.4f}]")
    print(f"σ: [{sigma_ci[0]:.4f}, {sigma_ci[1]:.4f}]")
    
    # 检查真实值是否在置信区间内
    mu_in_ci = mu_ci[0] <= true_mu <= mu_ci[1]
    sigma_in_ci = sigma_ci[0] <= true_sigma <= sigma_ci[1]
    
    print(f"真实μ在置信区间内: {mu_in_ci}")
    print(f"真实σ在置信区间内: {sigma_in_ci}")

if __name__ == "__main__":
    main()

打印结果：

开始生成模拟数据...
真实参数: μ = 5.0, σ = 2.0
样本数量: 1000

计算样本统计量...
样本均值: 5.0387
样本标准差: 1.9575

使用MLE估计参数...
MLE估计结果: μ = 5.0387, σ = 1.9575
与真实参数的差异: Δμ = 0.0387, Δσ = 0.0425

生成可视化图表...
libpng warning: iCCP: known incorrect sRGB profile
libpng warning: iCCP: known incorrect sRGB profile
libpng warning: iCCP: known incorrect sRGB profile
libpng warning: iCCP: known incorrect sRGB profile
libpng warning: iCCP: known incorrect sRGB profile

95%置信区间:
μ: [4.9173, 5.1600]
σ: [1.8717, 2.0432]
真实μ在置信区间内: True
真实σ在置信区间内: True

这个有个值得怀疑的问题，那就是：
Q6：为什么？对于正态分布，用样本均值$\bar{x}=5.0387$和样本标准差$s=1.9575$作为初始值，MLE优化后得到的估计值${\hat{\mu }}_{MLE}$和${\hat{\sigma }}_{MLE}$竟然还是$\bar{x}$和$s$。这是有什么错误吗？
解：这不是错误或巧合，而是必然。推导一下：
上文已求得正态分布的对数似然函数：
$$\log L(\mu ,\sigma |X)=-\frac{n}{2}\log (2\pi )-n\log (\sigma )-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(x_i-\mu {)}^2$$

为了找到极大值点，分别对$\mu$和$\sigma$求偏导并令其为零。

• 对$\mu$求偏导：
  $$\frac{\partial \log L}{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{i=1}^n(x_i-\mu )=0$$
  解得：
  $$\sum _{i=1}^n(x_i-\mu )=0⟹\sum _{i=1}^n{x}_i=n\mu ⟹{\hat{\mu }}_{MLE}=\frac{1}{n}\sum _{i=1}^n{x}_i=\bar{x}$$
• 对$\sigma$求偏导：
  $$\frac{\partial \log L}{\partial \sigma }=-\frac{n}{\sigma }+\frac{1}{{\sigma }^3}\sum _{i=1}^n(x_i-\mu {)}^2=0$$
  解得：
  $$-n{\sigma }^2+\sum _{i=1}^n(x_i-\mu {)}^2=0⟹{\hat{\sigma }}_{MLE}^2=\frac{1}{n}\sum _{i=1}^n(x_i-{\hat{\mu }}_{MLE}{)}^2$$

显然，经过数学推导的结果与程序结果一致。值得注意的是，这里MLE推导结果是有偏估计。
Q7：为什么是有偏估计呢？

有意思，接下来我想知道，如果我给优化算法一个很垃圾的初始值，比如假设均值是0，标准差是1，MLE优化后得到的估计值${\hat{\mu }}_{MLE}$和${\hat{\sigma }}_{MLE}$还是不是$\bar{x}$和$s$？
Q8：给定一个垃圾初始值，可视化一下算法估计参数的路径。

---

3. 贝叶斯派的观点

贝叶斯派认为$p(x|\theta )$中的$\theta$不是一个常量，而是一个随机变量，满足一个预设的先验的概率分布$\theta$ ~ $p(\theta )$。于是，根据贝叶斯定理，把参数$\theta$的后验$p(\theta |X)$和先验$p(\theta )$通过似然$p(X|\theta )$联系起来：
$$p(\theta |X)=\frac{p(X|\theta )\cdot p(\theta )}{p(X)}=\frac{p(X|\theta )\cdot p(\theta )}{\underset{\begin{array}{c}\theta \end{array}}{\int }p(X|\theta )\cdot p(\theta )d\theta }\text{\hspace{1em}(3-1)}$$
其中，$\underset{\begin{array}{c}\theta \end{array}}{\int }p(X|\theta )\cdot p(\theta )d\theta$是全概率公式的连续形式。一般来说，这个积分很难计算：

• 高维积分：当 $\theta$ 包含多个参数时，这是一个多重积分。
• 函数繁杂：似然函数和先验分布可能是很复杂的函数。
• 解析解难求：只有在"共轭先验"这种特殊情况下才能得到解析解。

因此，一般需要采用其他方法解决问题，比如：

• 数值积分：在低维情况下使用。
• 变分推断：寻找一个近似分布。
• MCMC方法：通过采样来近似这个积分。

这些方法我还不懂，先做个了解算了。

为了求$\sigma$的值，我们先引入一个参数估计方法——最大后验（Maximum A Posteriori, MAP）。这个方法和MLE是非常类似的：
$${\theta }_{MAP}=\underset{\theta }{\mathrm{argmax}}p(\theta |X)=\underset{\theta }{\mathrm{argmax}}p(X|\theta )\cdot p(\theta )\text{\hspace{1em}(3-2)}$$

其中，第二个等号是由于分母和$\theta$没有关系。可是$p(X)=\underset{\begin{array}{c}\theta \end{array}}{\int }p(X|\theta )\cdot p(\theta )d\theta$，我打眼一看，谁说没有关系？

首先要明确的是，在优化问题中，什么是“与$\theta$无关”？那就是：当我们在优化过程中改变θ的值时，分母的值不会随之变化。

若$X$是离散型变量，则$p(X)=\sum _{\theta }p(X|\theta )\cdot p(\theta )$，若X是连续型变量，则$p(X)=\underset{\begin{array}{c}\theta \end{array}}{\int }p(X|\theta )\cdot p(\theta )d\theta$。其实不论是离散还是连续，在求和、求积分的过程中，都要取遍$\theta$所有可能的值，最终$p(X)$得到的结果是一个常数$c$。而我们的优化变量$\theta$只是这所有可能的值中的、能够使得函数取最大值的，那一个特例。

所以，虽然$p(X)$看上去与$\theta$有关，但最后的结果总是一个常数$c$，优化变量$\theta$不论取谁，$p(X)$的值不变，因为会遍历所有的$\theta$值来计算，优化变量$\theta$爱怎么取怎么取。

另外，从数学证明的角度看，对于任意两个$\theta$值，例如比较${\theta }_1$和${\theta }_2$的后验概率：
$$\frac{p({\theta }_1|X)}{p({\theta }_2|X)}=\frac{\frac{p(X|{\theta }_1)p({\theta }_1)}{p(X)}}{\frac{p(X|{\theta }_2)p({\theta }_2)}{p(X)}}=\frac{p(X|{\theta }_1)p({\theta }_1)}{p(X|{\theta }_2)p({\theta }_2)}$$
发现了吗？分母$p(X)$被约掉了。所以，比较两个$\theta$值的后验概率的相对大小时，根本不需要知道$p(X)$的具体值。因此，我们得出：
$${\theta }_{MAP}=\underset{\theta }{\mathrm{argmax}}p(\theta |X)=\underset{\theta }{\mathrm{argmax}}p(X|\theta )\cdot p(\theta )$$

求解出这个$\theta$值后，代入$\frac{p(X|\theta )\cdot p(\theta )}{\underset{\begin{array}{c}\theta \end{array}}{\int }p(X|\theta )\cdot p(\theta )d\theta }$，就得到了参数的后验分布$p(\theta |X)$。

得到了参数的后验分布后，我们可以将这个分布用于贝叶斯预测：
$$p(x_{new}|X)=\underset{\theta }{\int }p(x_{new},\theta |X)d\theta =\underset{\theta }{\int }p(x_{new}|\theta )\cdot p(\theta |X)d\theta \text{\hspace{1em}(3-3)}$$

这里的思想是，想用老数据$X$预测新数据$x_{new}$，中间让$\theta$作为一个桥梁。我们以“从联合分布到边际分布”的视角来理解这个公式。$p(x_{new},\theta |X)$ 是新数据和参数的联合后验分布。如果我们想只知道新数据的分布，而不管参数的具体值，就需要把参数"积分掉"。

例如：

• 我们知道联合概率 $p($身高, 体重$)$。
• 如果只想知道身高的分布 $p($身高$)$，就需要对体重积分：$p($身高$)=\int p($身高, 体重$)d$体重。

这么一理解，贝叶斯预测就简单许多了。接下来，上代码！

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from matplotlib.widgets import Slider, Button   # 交互式滑块和按钮

# 设置中文字体和样式
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.figsize'] = [12, 8]

# 1. 生成模拟数据
np.random.seed(42)
n_samples = 100  # 增加样本数量以获得更好的可视化效果
true_slope = 2.5    # 真实斜率
true_intercept = 1.0    # 真实截距
noise_std = 1.5     # 噪声标准差

# 生成x和y数据
x = np.linspace(-3, 3, n_samples)   # 在[-3, 3]区间生成100个均匀分布的点
y = true_slope * x + true_intercept + np.random.normal(0, noise_std, n_samples) # 生成带噪声的线性数据

我以确定性的线性关系，然后加上服从正态分布的随机噪声（均值为0，标准差为1.5），来生成模拟数据。噪声模拟现实世界中观测值的不确定性。这个线性回归模型可以表示为：
$$y=w_0+w_{1}x+ϵ\text{\hspace{1em}(3-4)}$$
其中，$y$是因变量，$x$是自变量，$w_0$是截距,$w_1$是斜率，$ϵ$是噪声项，服从正态分布：
$$ϵ\sim \mathcal{N}(0,{\sigma }^2)$$

# 2. 设置先验分布p(θ)
prior_mean = np.array([0.0, 0.0])  # 先验均值μ₀，假设斜率和截距都接近0
prior_cov = np.array([[2.0, 0.0],  # 先验协方差Σ₀矩阵
                      [0.0, 2.0]])  # 假设斜率和截距独立，方差都为2

根据公式(3-1)，下一步我们要设计先验$p(\theta )$。因为我们有两个优化参数截距$w_0$和斜率$w_1$，所以需要设计二维高斯分布。在一维情况下，高斯分布由两个参数定义：$\mu$和$\sigma$。在多维高斯分布中就变为了：均值向量$\mu$和协方差矩阵$\Sigma$，均值向量描述每个维度的均值，协方差矩阵描述变量间的相关性。

请注意，我并不是因为噪声属于高斯分布，在此设置先验分布时就选择高斯分布，实际上两者之间没什么关系。贝叶斯方法的精髓在于：所有未知量都应该用概率分布来表示其不确定性。我为截距$w_0$和斜率$w_1$设置高斯先验，是因为高斯分布的一些优势，我选择使用高斯分布来估计截距$w_0$和斜率$w_1$：

1. 共轭先验（先验分布和后验分布属于同一分布族）。对于线性回归，似然函数是高斯分布，如果选择高斯先验，那么后验分布也是高斯分布。这允许我们得到解析解，而不需要复杂的数值积分。
2. 最大熵原理。在给定均值和方差的约束下，高斯分布是最不确定的分布（具有最大熵）。也就是说，我们只表达了均值和方差的信息，没有引入额外的假设。
3. 计算效率。只需要存储均值和协方差。

首先，我们有两个要估计的参数，一般来说默认为0，所以:
$$\theta =\left[\begin{array}{c}{w}_0\\ w_1\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

假设斜率和截距之间没有线性相关，且方差设定为2.0：
$$\Sigma =\left[\begin{array}{cc}Var(w_0)& Cov(w_0,w_1)\\ Cov(w_1,w_0)& Var(w_1)\end{array}\right]=\left[\begin{array}{cc}2.0& 0.0\\ 0.0& 2.0\end{array}\right]$$

之所以设置方差为2.0，是因为我们并不强烈相信均值为0，所以设置了2.0的方差。

数学公式可表达为：
$$p(\theta )=\mathcal{N}(\theta |{\mu }_0,{\Sigma }_0)$$

我们的观测集$X$，也就是我们的模拟数据，是根据线性回归+高斯噪声生成的，这造成了我在选择预测模型时不得不选择贝叶斯线性回归，而估计参数也不得不确定为截距$w_0$和斜率$w_1$。而实际情况中并不是这样，实际中，我们需要对观测集$X$进行探索性分析，如果数据的线性关系比较明显、数据量较小、有先验知识以及需要量化不确定性时，贝叶斯线性回归模型会适合使用。而且，实际应用往往会选择出多个模型，最终预测模型需要根据验证集的效果来确定。

因此，选择了贝叶斯线性回归作为预测模型，那么似然函数 $p(X|\theta )$ 就相应确定为线性模型加上高斯噪声的形式，也就是公式(3-4)。

在公式(3-1)中，$X$代表完整的观测数据集，包括输入特征和对应的输出值。在代码中，这对应于$(x,y)$对。这可以分解为 $p(X|\theta )=p(y|x,\theta )p(x|\theta )=p(y|x,\theta )$，代码中，输入特征$x$是固定的与$\theta$无关。于是，此处的似然函数公式：
$$p(y|X,\theta )=\mathcal{N}(y|Xw,{\sigma }^{2}I)$$
其中，$x$表示单个数据量，为了使用多个数据点，我们设计数据矩阵$X$，包括所有数据特征：
$$X=\left[\begin{array}{cc}1& x_1\\ 1& x_2\\ ⋮& ⋮\\ 1& x_n\end{array}\right]$$

接下来，就是计算公式(3-2)。最大化这个参数后验MAP。

# 3. 定义贝叶斯更新函数
def bayesian_update(x, y, prior_mean, prior_cov, noise_std):
    """
    贝叶斯线性回归的解析解更新
    返回后验均值和协方差
    """
    # 设计矩阵
    X_design = np.column_stack([np.ones_like(x), x])    # 对应 X
    
    # 计算后验分布参数
    # 计算先验精度矩阵（协方差矩阵的逆）
    prior_precision = np.linalg.inv(prior_cov)  # 对应 Σ₀⁻¹
    # 计算噪声精度（方差的倒数）
    noise_precision = 1.0 / (noise_std ** 2)    # 对应 1/σ²
    # 计算后验协方差矩阵
    posterior_cov = np.linalg.inv(prior_precision + noise_precision * X_design.T @ X_design) # 对应 Σ_N = (Σ₀⁻¹ + (1/σ²)XᵀX)⁻¹
    # 计算后验均值
    posterior_mean = posterior_cov @ (prior_precision @ prior_mean + noise_precision * X_design.T @ y) # 对应 μ_N = Σ_N(Σ₀⁻¹μ₀ + (1/σ²)Xᵀy)
    
    return posterior_mean, posterior_cov

在贝叶斯线性回归的设定中，我们假设先验分布和似然函数都是高斯分布，因此后验分布也是高斯分布。高斯分布有一个很好的性质：它的众数（即最大值点，也就是MAP估计）恰好等于均值（也等于中位数）。因此，我们不需要使用梯度下降或其他优化算法来寻找最大值，而是可以直接通过解析公式计算出后验均值，这个后验均值就是MAP估计。

说人话就是，对于高斯后验分布，均值就是MAP估计，对应代码中的posterior_mean。

对应数学公式：
$$p(\theta |X,y)=\mathcal{N}(\theta |{\mu }_N,{\Sigma }_N)$$
其中：
$${\Sigma }_N={\left({\Sigma }_0^{-1}+\frac{1}{{\sigma }^2}{X}^{T}X\right)}^{-1}$$
$${\theta }_{MAP}={\mu }_N={\Sigma }_N\left({\Sigma }_0^{-1}{\mu }_0+\frac{1}{{\sigma }^2}{X}^{T}y\right)$$

下面一股脑给出剩余的代码：

# 4. 计算不同数据量下的后验分布
sample_sizes = [0, 5, 10, 20, 50, 100]  # 要显示的数据点数量
posterior_means = []    # 储存不同样本量下的后验均值
posterior_covs = []     # 储存不同样本量下的后验协方差

for n in sample_sizes:
    if n == 0:
        # 先验分布（没有数据时）
        posterior_means.append(prior_mean)
        posterior_covs.append(prior_cov)
    else:
        # 使用前n个数据点更新后验分布
        mean, cov = bayesian_update(x[:n], y[:n], prior_mean, prior_cov, noise_std)
        posterior_means.append(mean)
        posterior_covs.append(cov)

# 转换为numpy数组方便后续处理
posterior_means = np.array(posterior_means)
posterior_covs = np.array(posterior_covs)

# 5. 创建交互式可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
plt.subplots_adjust(bottom=0.25)  # 为滑块留出空间

# 左上图：数据空间和回归线
scatter = ax1.scatter([], [], alpha=0.7, label='观测数据')  # 创建空散点图
regression_lines = []   # 储存多条回归线
for _ in range(20):
    line, = ax1.plot([], [], 'r-', alpha=0.2)   # 创建20条透明的红色回归线
    regression_lines.append(line)
mean_line, = ax1.plot([], [], 'b-', linewidth=2, label='平均预测')
true_line, = ax1.plot([], [], 'g--', linewidth=2, label='真实关系')

ax1.set_xlim(-3.5, 3.5)
ax1.set_ylim(-8, 12)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title('数据空间和回归线')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 右上图：预测分布
pred_line, = ax2.plot([], [], 'purple', linewidth=2, label='预测分布')
pred_fill = ax2.fill_between([], [], [], alpha=0.3, color='purple') # 预测分布填充区域
pred_mean_line = ax2.axvline(x=0, color='purple', linestyle='--', label='预测均值')
true_value_line = ax2.axvline(x=0, color='green', linestyle='--', label='真实值')

ax2.set_xlabel('y')
ax2.set_ylabel('概率密度')
ax2.set_title('预测分布 (x=0)')
ax2.set_ylim(0, 0.6)
ax2.legend()
ax2.grid(True, alpha=0.3)

# 添加滑块控制数据点数量
ax_slider = plt.axes([0.2, 0.1, 0.6, 0.03]) # 定义滑块位置
n_slider = Slider(
    ax=ax_slider,
    label='数据点数量',
    valmin=0,
    valmax=len(sample_sizes)-1,
    valinit=0,
    valstep=1,
)

# 添加测试点位置滑块
ax_test_slider = plt.axes([0.2, 0.05, 0.6, 0.03])
test_slider = Slider(
    ax=ax_test_slider,
    label='测试点 x 位置',
    valmin=-3.0,
    valmax=3.0,
    valinit=0.0,
    valstep=0.1,
)

# 更新函数
def update(val):
    # 获取当前滑块值
    idx = int(n_slider.val) # 数据点数量索引
    n = sample_sizes[idx]   # 实际数据点数量
    test_x = test_slider.val    # 测试点x位置
    
    # 获取当前后验分布
    current_mean = posterior_means[idx] # 当前后验均值
    current_cov = posterior_covs[idx]   # 当前后验协方差
    
    # 更新数据空间图
    scatter.set_offsets(np.column_stack([x[:n], y[:n]]))    # 更新散点图数据
    
    # 更新回归线
    x_test = np.linspace(-3.5, 3.5, 100)    # 生成测试x值
    for line in regression_lines:
        # 从后验分布中采样参数
        w_sample = np.random.multivariate_normal(current_mean, current_cov)
        # 计算预测y值
        y_test = w_sample[0] + w_sample[1] * x_test
        # 更新回归线
        line.set_data(x_test, y_test)
    
    # 更新平均预测线
    y_mean = current_mean[0] + current_mean[1] * x_test # 使用后验均值计算预测
    mean_line.set_data(x_test, y_mean)
    
    # 更新真实关系线
    y_true = true_intercept + true_slope * x_test   # 计算真实关系
    true_line.set_data(x_test, y_true)
    
    # 更新预测分布图
    X_design_star = np.array([1, test_x])   # 测试点的设计矩阵
    pred_mean = current_mean @ X_design_star    # 预测均值
    pred_var = X_design_star @ current_cov @ X_design_star.T + noise_std**2 # 预测方差
    
    # 生成预测分布曲线
    x_pred = np.linspace(pred_mean - 3*np.sqrt(pred_var), pred_mean + 3*np.sqrt(pred_var), 100)
    y_pred = stats.norm.pdf(x_pred, pred_mean, np.sqrt(pred_var))   # 正态分布概率密度
    
    pred_line.set_data(x_pred, y_pred)  # 更新预测分布曲线
    
    # 更新填充区域
    ax2.collections.clear() # 清除之前的填充
    ax2.fill_between(x_pred, 0, y_pred, alpha=0.3, color='purple')  # 重新填充
    
    # 更新预测均值线
    pred_mean_line.set_xdata([pred_mean, pred_mean])
    
    # 更新真实值线
    true_value = true_intercept + true_slope * test_x   # 计算真实值
    true_value_line.set_xdata([true_value, true_value])
    
    # 更新标题
    ax1.set_title(f'数据空间和回归线 (n={n})')
    ax2.set_title(f'预测分布 (x={test_x:.1f})')
    
    fig.canvas.draw_idle()  # 重绘图形

# 注册更新函数
n_slider.on_changed(update)
test_slider.on_changed(update)

# 添加重置按钮
reset_ax = plt.axes([0.8, 0.025, 0.1, 0.04])    # 定义按钮位置
reset_button = Button(reset_ax, '重置', hovercolor='0.975') # 创建重置按钮

def reset(event):
    n_slider.reset()    # 重置数据点数量滑块
    test_slider.reset() # 重置测试点位置滑块
reset_button.on_clicked(reset)  # 重置注册函数

# 初始化
update(0)   # 初始更新

plt.show()

# 6. 额外可视化：参数收敛情况
plt.figure(figsize=(10, 6))

# 计算所有可能数据量下的后验均值
all_means = []  # 储存所有后验均值
for n in range(0, n_samples+1):     # 从0到100个数据点
    if n == 0:
        all_means.append(prior_mean)    # 先验均值
    else:
        mean, _ = bayesian_update(x[:n], y[:n], prior_mean, prior_cov, noise_std)   # 更新后验
        all_means.append(mean)

all_means = np.array(all_means) # 转换为numpy按钮

# 绘制参数收敛图
plt.plot(range(n_samples+1), all_means[:, 0], 'b-', label='截距 (w0)')
plt.plot(range(n_samples+1), all_means[:, 1], 'r-', label='斜率 (w1)')
plt.axhline(y=true_intercept, color='b', linestyle='--', alpha=0.7, label='真实截距')
plt.axhline(y=true_slope, color='r', linestyle='--', alpha=0.7, label='真实斜率')

plt.xlabel('数据点数量')
plt.ylabel('参数值')
plt.title('参数收敛路径')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 7. 打印最终结果
final_mean = all_means[-1]  # 最终后验均值
final_cov = bayesian_update(x, y, prior_mean, prior_cov, noise_std)[1]  # 最终后验协方差

print("贝叶斯线性回归结果")
print("=" * 50)
print(f"真实参数: 截距 = {true_intercept}, 斜率 = {true_slope}, 噪声标准差 = {noise_std}")
print(f"先验参数: 均值 = {prior_mean}, 协方差矩阵 = \n{prior_cov}")
print(f"样本数量: {n_samples}")

print("\n最终后验分布:")
print(f"后验均值: {final_mean}")
print(f"后验协方差矩阵: \n{final_cov}")

# 计算参数估计误差
w0_error = abs(final_mean[0] - true_intercept)  # 截距误差
w1_error = abs(final_mean[1] - true_slope)  # 斜率误差

print(f"\n参数估计误差:")
print(f"截距误差: {w0_error:.4f}")
print(f"斜率误差: {w1_error:.4f}")

这部分代码比较繁琐，主要是设计可视化结果，占据了大量的篇幅，可视化是辅助理解的东西，这部分代码没必要看，我只说明两个重要的点。

按照流程，得到了参数的后验分布后，我们可以将这个分布用于贝叶斯预测。

    # 更新预测分布图
    X_design_star = np.array([1, test_x])   # 测试点的设计矩阵
    pred_mean = current_mean @ X_design_star    # 预测均值
    pred_var = X_design_star @ current_cov @ X_design_star.T + noise_std**2 # 预测方差
    
    # 生成预测分布曲线
    x_pred = np.linspace(pred_mean - 3*np.sqrt(pred_var), pred_mean + 3*np.sqrt(pred_var), 100)
    y_pred = stats.norm.pdf(x_pred, pred_mean, np.sqrt(pred_var))   # 正态分布概率密度

这几行代码是贝叶斯预测的核心。在贝叶斯线性回归中，如果选择高斯似然和高斯先验，那么后验分布和预测分布都是高斯分布。经过所有计算后，我们得到了预测分布的两个关键参数：预测分布的均值${\mu }_{pred}$和预测分布的方差${\sigma }_{pred}^2$。一旦我们有了这些参数，对于任何新的输入 $x_{new}$，预测就变得非常简单：

1. 构建设计向量。$x_{new}^{\ast }=[1,x_{new}{]}^T$

    X_design_star = np.array([1, test_x])   # 测试点的设计矩阵

2. 计算预测均值。${\mu }_{pred}=(x_{new}^{\ast }{)}^T{\mu }_N$

    pred_mean = current_mean @ X_design_star    # 预测均值

3. 计算预测方差。${\sigma }_{pred}^2=(x_{new}^{\ast }{)}^T{\Sigma }_N{x}_{new}^{\ast }+{\sigma }^2$

    pred_var = X_design_star @ current_cov @ X_design_star.T + noise_std**2 # 预测方差

4. 生成预测分布范围。生成覆盖预测分布 $\pm 3\sigma$ 范围的$x$值，这包含了约99.7%的概率质量（根据正态分布的3σ原则）

    x_pred = np.linspace(pred_mean - 3*np.sqrt(pred_var), pred_mean + 3*np.sqrt(pred_var), 100)

5. 计算概率密度。$p(y_{new}|x_{new},X)=\frac{1}{\sqrt{2\pi {\sigma }_{pred}^2}}\exp \left(-\frac{(y_{new}-{\mu }_{pred}{)}^2}{2{\sigma }_{pred}^2}\right)$

    y_pred = stats.norm.pdf(x_pred, pred_mean, np.sqrt(pred_var))   # 正态分布概率密度

这部分整体的数学表达式为：
$$p(y_{new}|x_{new},X,y)=\mathcal{N}(y_{new}|{\mu }_{pred},{\sigma }_{pred}^2)\text{\hspace{1em}(3-5)}$$

这部分代码最终得到的是预测分布的概率密度函数。如果想得到点预测，直接使用 ${\mu }_{pred}$ 作为 $y_{new}$ 的预测值，也就是pred_mean。

有个疑问啊，看公式(3-3)是积分形式，为什么这个案例中没有涉及到积分呢？
Q9：怎么没有涉及到积分呢？
解：对于线性回归和高斯先验的特殊情况，这个积分是通过解析解而不是数值积分来计算的，这个积分有闭合形式的解，就是(3-5)式。这个解析解已经包含了积分的效果，数学上可以证明：
$$\int p(y_{new}|x_{new},\theta )p(\theta |X,y)d\theta =\mathcal{N}(y_{new}|x_{new}^T{\mu }_N,x_{new}^T{\Sigma }_N{x}_{new}+{\sigma }^2)=\mathcal{N}(y_{new}|{\mu }_{pred},{\sigma }_{pred}^2)$$
所以代码中直接计算的就是积分的结果。

---

另一个比较有意思的点，是这部分代码：

# 4. 计算不同数据量下的后验分布
sample_sizes = [0, 5, 10, 20, 50, 100]  # 要显示的数据点数量
posterior_means = []    # 储存不同样本量下的后验均值
posterior_covs = []     # 储存不同样本量下的后验协方差

for n in sample_sizes:
    if n == 0:
        # 先验分布（没有数据时）
        posterior_means.append(prior_mean)
        posterior_covs.append(prior_cov)
    else:
        # 使用前n个数据点更新后验分布
        mean, cov = bayesian_update(x[:n], y[:n], prior_mean, prior_cov, noise_std)
        posterior_means.append(mean)
        posterior_covs.append(cov)

这部分可视化代码，通过不同数据量下的后验分布变化来演示贝叶斯学习的关键特性。体现出来贝叶斯统计的核心思想：随着数据量的增加，似然函数的影响逐渐超过先验分布。

• n=0：只有先验分布，没有数据影响。
• n=5：少量数据，开始影响后验。
• n=100：大量数据，后验高度集中于真实参数附近。
  

图3 n=0

图4 n=20

图5 n=100

图6 参数收敛路径

图片中展示出来的信息直观明了，看看就好。到了该结束的时候了，我想，回头再看看频率派和贝叶斯派的本质区别，会更加容易理解。

频率派和贝叶斯派的本质区别

可以这么理解：贝叶斯学派发展的模型（如概率图模型）核心是计算后验分布，这通常涉及积分问题。当解析解不可得实，会采用数值积分方法，例如蒙特卡洛采样（如MCMC）来近似积分。而频率派发展的模型（通常称为统计机器学习）本质上是优化算法，步骤包括：第一步设计模型，第二步定义损失函数，第三步通过优化算法（例如梯度下降或牛顿法）最小化损失函数，从而将问题转化为求解最优解。

---

问题解答

Q1：贝叶斯派经典模型有哪些？概率图模型是一个统称还是个例？
解：经典的贝叶斯模型有：贝叶斯线性回归、朴素贝叶斯分类器、潜在狄利克雷分配、高斯过程和贝叶斯神经网络。

概率图模型是一个统称，用图这种数据结构来表达变量之间的复杂概率依赖关系。许多贝叶斯模型，如潜在狄利克雷分配，使用概率图模型来直观地定义其结构并进行推断。同时，概率图模型也可以用于表示频率派模型，如马尔科夫随机场中的参数通常由频率派方法学习。

也就是说，概率图模型是一种思想，是一种框架。此处应该是想表达：从贝叶斯角度出发，常常利用概率图模型这一强大工具来构建和推断复杂的模型。

Q2：简单了解，蒙特卡洛方法是啥，干啥的？MCMC？
解：蒙特卡洛方法的核心思想是：与其绞尽脑汁去推导一个复杂的公式，不如用计算机生成大量随机数来模拟它，最终得到一个足够精确的近似解。一个经典例子是计算圆周率$\pi$，在一个边长为2的正方形中画一个半径为1的圆，随机往里面扔豆子，豆子落在圆内的概率近似为圆周率。

回想一下贝叶斯分析，我们的目标是计算后验分布，但这个分布往往没有一个标准的形态，像一座崎岖不平的山，这就形成了简单蒙特卡洛困境：我想往“后验分布”这座山上扔豆子，通过豆子的密集程度来描绘山的形状，但通常分布函数很复杂，没办法直接扔出符合它形状的豆子。

MCMC全程马尔科夫链蒙特卡洛，就是解决这个问题的。就了解到这儿吧，了解深了也没啥价值。

Q3：频率派经典模型有哪些？
解：线性模型及其扩展：线性回归、逻辑回归。基于距离和间隔的模型：K-最近邻、支持向量机。基于树和集成学习的模型：决策树、随机森林、梯度提升机（如XGBoost）。无监督学习模型：K-均值聚类。

Q7：为什么是有偏估计呢？
解：因为用数据估计参数会消耗信息。在估计$\sigma$之前，它已经用了同一个样本估计均值$\mu$，消耗了一个自由度。用样本均值来估计总体均值时：
$$\hat{\mu }=\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$$

这个过程实际上是在利用数据满足一个约束条件：所有数据点的平均值等于$\bar{x}$。这意味着，一旦固定了$\bar{x}和$$n-1$个数据值的点，第$n$个数据点的值就不能自由变化了，它必须满足：
$$x_n=n\bar{x}-\sum _{i=1}^{n-1}{x}_i$$
是不是有点绕，简单来说，你确定了$n-2$个值，还剩下两个值，只要满足约束，你仍然有很多组合来选择这两个值，但确定了$n-1$个值时，最后一个值在约束下已经被确定了，你没得选择！

这就是消耗了一个自由度，你只剩下$n-1$个自由的数据点来估计剩余的$\sigma$了。体现在$\sigma$的估计上就是：因为你不知道真实的$\mu$值，只能用先估计出来的$\hat{\mu }=\bar{x}$来代替：
$${\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(x_i-{\bar{x}}^2)$$

这就会产生一个问题，因为一个数学性质：样本均值$\bar{x}$，是使得平方偏差之和最小的那个值（不证明，当一个结论记忆算了，在此篇文章中不重要），所以有：
$$\sum _{i=1}^n(x_i-\bar{x}{)}^2\le \sum _{i=1}^n(x_i-\mu {)}^2$$
对于任何$\mu \ne \bar{x}$都成立。综上，用$\bar{x}$代替$\mu$会系统性地低估真实的平方偏差之和，这就是有偏的含义。

那么为什么当分母是$n-1$，即$\frac{1}{n-1}{\sum }_{i=1}^n(x_i-\bar{x}{)}^2$就是无偏估计了呢？简单证明一下：虽然$\sum (x_i-\bar{x}{)}^2$仍然会系统性地偏小，但是它的数学期望满足：
$$E\left[\sum _{i=1}^n(x_i-\bar{x}{)}^2\right]=(n-1){\sigma }^2$$
因此，如果我们使用：
$$s^2=\frac{1}{n-1}\sum _{i=1}^n(x_i-\bar{x}{)}^2$$
那么这个估计量的期望值就是：
$$E[s^2]=E\left[\frac{1}{n-1}\sum _{i=1}^n(x_i-\bar{x}{)}^2\right]=\frac{1}{n-1}\times (n-1){\sigma }^2={\sigma }^2$$
这就成为了一个无偏估计。解释到此为止吧，因为我也是半瓶子醋哈哈，况且再深入就跑题了。

Q8：给定一个垃圾初始值，可视化一下算法估计参数的路径。
解：这部分直接上代码，这部分只是起到一个补充说明的作用，额外增加了一些可视化内容，所以在我看来，没必要仔细分析代码，不过，我仍然会在代码中详细注释。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.stats import norm
import seaborn as sns
from matplotlib.gridspec import GridSpec

# 设置中文字体和数学符号支持
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False    # 正确显示负号
plt.rcParams['mathtext.fontset'] = 'stix'     # 使用 STIX 字体

# 设置随机种子以确保结果可重现
np.random.seed(42)

# 1. 生成模拟数据
def generate_data(n_samples=1000, true_mu=5, true_sigma=2):
    """
    生成服从正态分布的模拟数据
    true_mu: 分布的均值
    true_sigma: 标准差
    n_samples: 样本数量
    """
    data = np.random.normal(true_mu, true_sigma, n_samples)
    return data

# 2. 定义似然函数和对数似然函数
def log_likelihood_normal(params, data):
    """
    正态分布的对数似然函数
    params: [mu, sigma]
    data: 观测数据
    """
    mu, sigma = params
    if sigma <= 0:  # 标准差必须为正，否则返回负无穷
        return -np.inf
    n = len(data)
    # 正态分布的对数似然函数的公式
    log_likelihood = -n/2 * np.log(2*np.pi) - n*np.log(sigma) - 1/(2*sigma**2)*np.sum((data - mu)**2)
    return log_likelihood

def negative_log_likelihood(params, data):
    """
    负对数似然函数（用于最小化）
    """
    return -log_likelihood_normal(params, data)

# 3. 使用优化算法进行极大似然估计（增强版，记录优化路径）
def mle_normal_with_path(data, initial_params=None):
    """
    使用MLE估计正态分布的参数，并记录优化路径
    """
    if initial_params is None:
        initial_params = [np.mean(data), np.std(data)]
    
    # 存储优化路径
    path = {
        'mu': [initial_params[0]],
        'sigma': [initial_params[1]],
        'nll': [negative_log_likelihood(initial_params, data)]
    }
    
    # 定义回调函数记录路径
    def callback(params):
        path['mu'].append(params[0])
        path['sigma'].append(params[1])
        path['nll'].append(negative_log_likelihood(params, data))
        return False
    
    # 最小化负对数似然函数
    result = minimize(
        negative_log_likelihood,
        initial_params,
        args=(data,),
        bounds=[(None, None), (1e-6, None)],
        method='L-BFGS-B',
        callback=callback
    )
    
    if result.success:
        # 添加最终结果到路径
        path['mu'].append(result.x[0])
        path['sigma'].append(result.x[1])
        path['nll'].append(negative_log_likelihood(result.x, data))
        return result.x, path
    else:
        raise ValueError("优化失败: " + result.message)

# 4. 增强版可视化函数
def plot_enhanced_results(data, true_params, estimated_params, path, initial_params):
    """
    增强版可视化：展示真实分布、估计分布、数据直方图和优化过程
    """
    true_mu, true_sigma = true_params
    est_mu, est_sigma = estimated_params
    init_mu, init_sigma = initial_params
    
    # 创建图形和网格布局
    fig = plt.figure(figsize=(18, 12))
    gs = GridSpec(2, 2, figure=fig, hspace=0.3, wspace=0.25)
    
    # 左上图：数据分布与参数估计
    ax1 = fig.add_subplot(gs[0, 0])
    sns.histplot(data, bins=30, stat='density', alpha=0.7, ax=ax1, label='数据分布')
    
    # 生成x值范围
    x_min, x_max = min(data), max(data)
    x = np.linspace(x_min, x_max, 1000)
    
    # 绘制真实分布
    true_pdf = norm.pdf(x, true_mu, true_sigma)
    ax1.plot(x, true_pdf, 'r-', linewidth=2, label=f'真实分布: N({true_mu:.2f}, {true_sigma:.2f}$^2$)')
    
    # 绘制估计分布
    est_pdf = norm.pdf(x, est_mu, est_sigma)
    ax1.plot(x, est_pdf, 'g--', linewidth=2, label=f'MLE估计: N({est_mu:.2f}, {est_sigma:.2f}$^2$)')
    
    ax1.set_title('数据分布与参数估计')
    ax1.set_xlabel('值')
    ax1.set_ylabel('密度')
    ax1.legend()
    
    # 右上图：参数收敛路径
    ax2 = fig.add_subplot(gs[0, 1])
    iterations = range(len(path['mu']))
    ax2.plot(iterations, path['mu'], 'b-', label='$\\mu$', linewidth=2)
    ax2.plot(iterations, path['sigma'], 'g-', label='$\\sigma$', linewidth=2)
    ax2.axhline(y=true_mu, color='r', linestyle='--', alpha=0.7, label='真实 $\\mu$')
    ax2.axhline(y=true_sigma, color='m', linestyle='--', alpha=0.7, label='真实 $\\sigma$')
    ax2.set_xlabel('迭代次数')
    ax2.set_ylabel('参数值')
    ax2.set_title('参数收敛路径')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 左下图：负对数似然值收敛
    ax3 = fig.add_subplot(gs[1, 0])
    ax3.plot(iterations, path['nll'], 'purple', linewidth=2)
    ax3.set_xlabel('迭代次数')
    ax3.set_ylabel('负对数似然值')
    ax3.set_title('负对数似然值收敛过程')
    ax3.grid(True, alpha=0.3)
    
    # 右下图：参数空间中的优化路径
    ax4 = fig.add_subplot(gs[1, 1])
    
    # 创建参数空间的网格
    mu_min = min(min(path['mu']), true_mu) - 1
    mu_max = max(max(path['mu']), true_mu) + 1
    sigma_min = max(0.1, min(min(path['sigma']), true_sigma) - 1)
    sigma_max = max(max(path['sigma']), true_sigma) + 1
    
    mu_vals = np.linspace(mu_min, mu_max, 50)
    sigma_vals = np.linspace(sigma_min, sigma_max, 50)
    M, S = np.meshgrid(mu_vals, sigma_vals)
    
    # 计算网格上每点的负对数似然值
    Z = np.zeros_like(M)
    for i in range(M.shape[0]):
        for j in range(M.shape[1]):
            Z[i, j] = negative_log_likelihood([M[i, j], S[i, j]], data)
    
    # 绘制等高线图
    contour = ax4.contour(M, S, Z, 20, colors='gray', alpha=0.5)
    ax4.clabel(contour, inline=True, fontsize=8)
    
    # 绘制优化路径
    ax4.plot(path['mu'], path['sigma'], 'ro-', linewidth=2, markersize=4, label='优化路径')
    ax4.plot(init_mu, init_sigma, 'bo', markersize=8, label='初始值')
    ax4.plot(true_mu, true_sigma, 'y*', markersize=12, label='真实参数')
    ax4.plot(est_mu, est_sigma, 'kx', markersize=10, label='最终估计')
    
    ax4.set_xlabel('$\\mu$')
    ax4.set_ylabel('$\\sigma$')
    ax4.set_title('参数空间中的优化路径')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()

# 5. 计算统计信息
def calculate_statistics(data, true_params, estimated_params):
    """
    计算并输出各种统计信息
    """
    true_mu, true_sigma = true_params
    est_mu, est_sigma = estimated_params
    
    # 计算基本统计量
    n = len(data)
    sample_mean = np.mean(data)
    sample_std = np.std(data)
    
    # 计算偏差
    mu_bias = est_mu - true_mu
    sigma_bias = est_sigma - true_sigma
    
    # 计算相对误差
    mu_rel_error = abs(mu_bias / true_mu) * 100
    sigma_rel_error = abs(sigma_bias / true_sigma) * 100
    
    # 计算置信区间
    mu_se = est_sigma / np.sqrt(n)
    sigma_se = est_sigma / np.sqrt(2 * n)
    
    mu_ci = (est_mu - 1.96 * mu_se, est_mu + 1.96 * mu_se)
    sigma_ci = (est_sigma - 1.96 * sigma_se, est_sigma + 1.96 * sigma_se)
    
    # 检查真实值是否在置信区间内
    mu_in_ci = mu_ci[0] <= true_mu <= mu_ci[1]
    sigma_in_ci = sigma_ci[0] <= true_sigma <= sigma_ci[1]
    
    # 输出结果
    print("="*50)
    print("极大似然估计结果分析")
    print("="*50)
    print(f"真实参数: μ = {true_mu:.4f}, σ = {true_sigma:.4f}")
    print(f"样本统计量: μ̄ = {sample_mean:.4f}, s = {sample_std:.4f}")
    print(f"MLE估计: μ = {est_mu:.4f}, σ = {est_sigma:.4f}")
    print(f"偏差: Δμ = {mu_bias:.4f}, Δσ = {sigma_bias:.4f}")
    print(f"相对误差: μ误差 = {mu_rel_error:.2f}%, σ误差 = {sigma_rel_error:.2f}%")
    print(f"95%置信区间:")
    print(f"  μ: [{mu_ci[0]:.4f}, {mu_ci[1]:.4f}] - 真实值{'在' if mu_in_ci else '不在'}区间内")
    print(f"  σ: [{sigma_ci[0]:.4f}, {sigma_ci[1]:.4f}] - 真实值{'在' if sigma_in_ci else '不在'}区间内")
    print("="*50)
    
    return {
        'mu_bias': mu_bias,
        'sigma_bias': sigma_bias,
        'mu_rel_error': mu_rel_error,
        'sigma_rel_error': sigma_rel_error,
        'mu_ci': mu_ci,
        'sigma_ci': sigma_ci,
        'mu_in_ci': mu_in_ci,
        'sigma_in_ci': sigma_in_ci
    }

# 6. 主程序
def main():
    # 设置真实参数
    true_mu = 5.0
    true_sigma = 2.0
    n_samples = 1000
    
    print("开始生成模拟数据...")
    print(f"真实参数: μ = {true_mu}, σ = {true_sigma}")
    print(f"样本数量: {n_samples}")
    
    # 生成数据
    data = generate_data(n_samples, true_mu, true_sigma)
    
    print("\n计算样本统计量...")
    sample_mean = np.mean(data)
    sample_std = np.std(data)
    print(f"样本均值: {sample_mean:.4f}")
    print(f"样本标准差: {sample_std:.4f}")
    
    print("\n使用MLE估计参数...")
    # 使用较差的初始值，以便更好地展示优化过程
    initial_params = [0.0, 10.0]  # 差的初始值：mu=0, sigma=10
    print(f"初始参数: μ = {initial_params[0]}, σ = {initial_params[1]}")
    
    # 使用增强版的MLE函数，记录优化路径
    estimated_params, path = mle_normal_with_path(data, initial_params)
    
    print(f"MLE估计结果: μ = {estimated_params[0]:.4f}, σ = {estimated_params[1]:.4f}")
    print(f"迭代次数: {len(path['mu']) - 1}")
    
    # 计算并显示统计信息
    stats = calculate_statistics(data, (true_mu, true_sigma), estimated_params)
    
    # 可视化结果
    print("\n生成可视化图表...")
    plot_enhanced_results(data, (true_mu, true_sigma), estimated_params, path, initial_params)

if __name__ == "__main__":
    main()

图2 参数收敛路径

右上图详细展示了参数的收敛路径，尽管我们参数的初始值给定的很离谱，但仍然能得到正确结果，没毛病。其余图做个参考，看个乐子就够了。

4. 结尾

本文的程序和图片资源：
通过网盘分享的文件：1_频率派VS贝叶斯派
链接: https://pan.baidu.com/s/1ki62WQBY9DC7g_YL12PJEA 提取码: 4wy8