数值统计 TP：高斯混合模型参数估计

1 随机变量的实现

1.1 均匀分布密度

实验目的

本部分先生成均匀分布的随机数，并用直方图方法估计概率密度。通过实验观察样本量N和区间数对密度估计效果的影响，同时验证蒙特卡罗估计的基本性质。是后续混合模型的基础。

问题1：生成N=100个均匀分布实现并显示

[图片1：N=100个均匀分布样本的散点图]

从图中可以看到，100个样本值在 [0,1] 区间内随机分布，没有明显的周期性或趋势，整体上覆盖了整个区间，符合均匀分布的特征。

问题2：从统计直方图推导密度估计

我们已知，随机变量 X 落在区间 [x - \frac{\delta}{2}, x + \frac{\delta}{2}] 的概率可以通过频率来近似：

同时根据概率密度函数的定义，随机变量 X 落在区间 [x - \frac{\delta}{2}, x + \frac{\delta}{2}] 的概率可以表示为该区间上的积分：

当 \delta 很小时，在区间 [x - \frac{\delta}{2}, x + \frac{\delta}{2}] 内，概率密度函数 p_X(t)​ 的积分可以近似为矩形面积：

因此我们有

解出 p_X(x)：

这是非参数估计中的直方图密度估计方法: 将频率除以区间宽度，得到概率密度的估计值

问题3：使用histcounts函数计算 N_x

统计落在每个小区间内的样本数量 N_x

问题4：显示密度估计并叠加真实密度

[图片2：密度估计直方图与真实密度对比]

蓝色柱状图是基于样本的密度估计，红色实线是真实密度 p_X(x)=1。各柱子高度在1.0上下波动，这是由于样本量有限导致的随机误差。随着样本量增加，波动会减小，估计值会更接近真实密度。

问题5：研究N、δ或区间数量的影响

设计两组对比实验来研究样本量和区间划分的影响。

实验一：样本量的影响

固定区间数为20，测试四种不同的样本量：N \in \{50, 100, 500, 1000\}。对于每个样本量，我们生成对应数量的样本并计算密度估计。

实验结果：

[图片3：样本量N的影响]

N=50时波动最大，N=1000时直方图已经很平坦，基本贴合红线。随着样本量增大，估计精度明显提升。

我们可以用大数定律来解释，每个区间内的样本数量 N_x 可以看作是 N 次伯努利试验的成功次数，其期望值为 Np，其中 p=\delta 是样本落入该区间的概率。根据大数定律，当 N 足够大时，N_x/N 会以概率1收敛到 p，因此密度估计 N_x/(N\delta) 会收敛到真实密度1。

实验二：区间数量的影响

固定 N=500，测试 \text{bins} = 5, 10, 20, 50 四种区间数。

实验结果：

[图片4：区间数量的影响]

区间数量的选择体现了非参数估计中经典的偏差-方差权衡问题。

当区间数为5时，区间宽度 \delta=0.2 较大，每个区间内包含了大量样本，因此估计比较稳定，但过于粗糙，无法捕捉密度函数的细节信息，产生较大的偏差。

随着区间数增加到10和20，估计的精细程度提高，既保持了一定的稳定性又能较好地反映密度特征。但当区间数增加到50时，区间宽度 \delta=0.02 过小，每个区间内的样本数量很少（平均每个区间只有10个样本），导致随机波动剧烈增加，虽然理论上能提供更精细的分辨率，但实际上方差过大使得估计变得不可靠。

问题6：计算X的均值和理论方差

根据实验设定，我们考虑的是 [0, 1] 上的均匀分布，其概率密度函数为：

计算均值：

根据期望的定义，均值 \mu 为：

计算该积分：

因此，\boxed{\mu = \frac{1}{2}}

计算方差：

方差的计算使用公式 \mathbb{V}[X] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2。首先计算 \mathbb{E}[X^2]：

代入方差公式：

因此，\boxed{v = \frac{1}{12}}

问题7：与经验均值和方差进行比较

使用 N=1000 个样本计算经验统计量，并与理论值进行对比

实验结果：

实验结果显示，经验统计量与理论值相吻合。

问题8：验证蒙特卡罗估计量的方差

进行 M=1000 次独立实验，每次生成 N=100 个样本并计算均值 \hat{\mu}_m，然后统计这1000个均值的方差。理论上，蒙特卡罗估计量的方差为：

实验结果：

问题7 和 8 的结果验证了蒙特卡罗方法的两个关键性质：无偏性（样本均值的期望等于总体均值）和一致性（样本量趋于无穷时，样本统计量收敛到总体参数）。

本节小结

样本量增加能显著改善估计质量，区间数选择需要在偏差和方差间权衡，最后我们验证了蒙特卡罗方法的无偏性和方差公式。

1.2 混合模型

实验目的

这部分学习混合模型的生成机制，比如如何从层次构造生成混合分布：先生成标签变量决定状态，然后根据标签从相应的条件分布中采样观测变量。

混合模型的基本概念

混合模型由两个层次的随机变量构成。第一层是标签变量 L，它是离散随机变量，描述不同的"状态"或"组分"。在二元情况下，L \in \{1,2\}，服从伯努利分布：

第二层是观测变量 X，它是连续随机变量。X 的分布依赖于标签 L 的取值，即 X 服从条件分布 p_{X|L}(x|l)。本TP考虑高斯混合模型，每个组分对应一个不同参数的高斯分布：

数据来自多个不同的来源或状态，每个状态有自己的分布特征，标签变量决定了数据来自哪个状态。

问题9：提出模拟L的方法并应用

为了模拟伯努利分布 L \sim \mathcal{B}(p)，我们采用逆变换采样法：生成 U \sim \text{Uniform}(0,1)，若 U < p 则 L=1，否则 L=2。

实现时先初始化所有标签为1，然后将满足 u \geq p 条件的样本改为2。生成 N=100 个样本，p=0.5。

正式算法名称： 这是逆变换采样法（Inverse Transform Sampling），对于离散分布特别简单：将 [0,1] 区间按概率划分，根据均匀随机数落在哪段来决定输出值。

实验结果：

[图片5：伯努利分布L的N=100次实现]

上图展示了100个标签值的时间序列，标签在1和2之间随机跳变。下图显示经验频率与理论概率的对比，蓝色柱状图是经验频率，红色柱状图是理论概率 p=0.5 和 1-p=0.5。在 N=100 的样本量下，两个标签出现频率大致相等，验证了采样方法的正确性。

问题10：从层次构造模拟边缘分布 p_X

根据混合模型的层次结构，用两步过程从边缘分布 p_X 中采样。第一步采样标签 L \sim \mathcal{B}(p)，第二步根据 L 的值从相应的条件高斯分布采样 X：

其中 Z \sim \mathcal{N}(0,1) 是标准正态随机变量。

这个方法基于全概率公式：

通过先采样 L，再给定 L 采样 X，我们实际上是从联合分布 p_{X,L}(x,l) 中采样，然后 忽略 标签 L 的值，只保留 X，这样得到的 X 就服从边缘分布 p_X(x)。

具体实现中，先用问题9的方法生成1000个标签，然后对每个样本根据其标签值从相应的高斯分布中采样。

问题11：绘制结果并验证

实验结果：

[图片6：混合高斯模型的四个子图展示]

样本序列图： 第一个子图显示1000个样本按生成顺序排列。样本在两个中心附近聚集，一个在 -2 附近（对应 \mu_1=-2），另一个在 3 附近（对应 \mu_2=3）

标签着色图： 第二个子图按标签着色，红色圆点（L=1）主要分布在 -2 附近，蓝色方点（L=2）主要分布在 3 附近。

边缘分布验证： 第三个子图展示所有样本的直方图，红色实线是理论混合密度

绿色和洋红色虚线分别表示两个组分的加权密度。

**条件分布验证：**第四个子图分别展示两个标签下样本的条件分布。红色直方图是 L=1 时 X 的分布，蓝色直方图是 L=2 时 X 的分布。根据标签分组后，每组样本都近似服从单峰的高斯分布，验证了条件分布 p_{X|L}(x|l) 的正确性。

定量验证结果：

L=1: 均值 -2.014 (理论 -2.000), 方差 1.077 (理论 1.000)
L=2: 均值 3.092 (理论 3.000), 方差 2.263 (理论 2.250)

两个组分的经验均值都非常接近理论值，可见验证结果证明了采样方法的正确性。

2. 混合模型参数估计

2.1 模型定义

本部分进入混合模型的参数估计问题。在前面的实验中，我们已知模型参数并从模型中生成数据。现在我们面对逆问题：给定观测数据去推断模型的参数，对于混合模型而言，这个问题更加复杂，因为我们通常无法直接观测到标签变量 L，这些隐藏的标签使得参数估计变成了一个不完全数据问题。

我们考虑的混合模型为：

其中标签 L 服从伯努利分布，观测变量 X 在给定标签时服从高斯分布：

这里我们使用精度（precision）\gamma_l = 1/\sigma_l^2 来参数化高斯分布，这在某些推导中更为方便。高斯分布的概率密度函数用精度表示为：

问题12：确定联合密度和边缘分布

联合密度可以写成条件密度与边缘密度的乘积：

其中，L 服从伯努利分布（备注1），X \mid L = l 服从高斯分布（方程3）。将这两个分布的表达式代入，我们得到联合密度：

边缘分布 \mathrm{p}_X(x) 通过对所有可能的 l 值求和（边缘化）得到：

当 l = 1 时：

当 l = 2 时：

因此，边缘分布为：

这正是两个高斯分布 \mathcal{N}(x \mid \mu_1, \gamma_1) 和 \mathcal{N}(x \mid \mu_2, \gamma_2) 的凸组合，因为权重 p 和 (1-p) 满足 p + (1-p) = 1 且 p, (1-p) \in [0, 1]。这是高斯混合模型（Gaussian Mixture Model）。

问题13：用特定参数抽取样本

为了研究不同参数配置下混合模型的行为特征，我们测试两组具有不同性质的参数组合。

参数组(a)：分离良好的两个高斯分布

第一组参数设置为 \mu_1=-1, \mu_2=1, \gamma_1=\gamma_2=1, p=0.5。转换为标准差形式为 \sigma_1=\sigma_2=1。这组参数的特点是两个高斯分布的均值相距2个标准差单位（|\mu_2-\mu_1|=2），方差相同，混合权重相等。这是理想的聚类场景，两个组分在空间上有相对清晰的分离。

我们生成 N=1000 个样本，首先按照 p=0.5 的概率生成标签，然后根据标签从相应的高斯分布中采样观测值。

实验结果：

[图片7：参数组(a)的密度分布图和标签着色散点图]

左图展示了混合分布的直方图和理论密度曲线。可以观察到双峰结构，两个峰分别位于 -1 和 1 附近，这是因为两个高斯分布的方差相同且混合权重相等，两个峰的高度和宽度大致相同，呈现出对称的形态。两个峰之间存在一个明显的谷（在 x=0 附近密度较低），这表明两个组分在空间上有相对清晰的分离。两个组分的贡献 p \cdot \mathcal{N}(x|\mu_1,\gamma_1) 和 (1-p) \cdot \mathcal{N}(x|\mu_2,\gamma_2) 相加得到红色的总密度曲线。样本的直方图与理论曲线相对应，验证了采样的正确性。

右图展示了按标签着色的样本散点图。两个组分之间有明显的分界，重叠区域很小，特征差异较大，因此后续即使不知道真实标签，通过观测值的分布特征也能较准确地推断样本的归属。

参数组(b)：重叠严重的两个高斯分布

第二组参数设置为 \mu_1 = \mu_2 = 0, \gamma_1 = 10, \gamma_2 = 0.1, p = 0.9。转换为标准差形式为 \sigma_1 = 0.316, \sigma_2 = 3.162。这组参数的特点是：

• 两个高斯分布共享相同的均值（都在0处）
• 方差相差100倍（\gamma_1/\gamma_2 = 100）
• 权重严重不平衡（90%的样本来自第一个组分）

实验结果：

[此处插入图片8：参数组(b)的密度分布图和标签着色散点图]

左图展示了这组参数下的混合分布。由于 p=0.9，90%的样本来自窄高斯，因此它贡献了主要的峰。底下的第二个高斯分布的方差很大（\sigma_2=3.162），密度曲线展得很宽，但由于混合权重只有10%，它的贡献相对较小，主要体现在尾部区域。

右图的散点图展示了数据的结构。绝大多数样本（红色圆点，L=1）聚集在 x=0 附近，少量样本（蓝色方点，L=2）则广泛分散。大部分是正常数据（集中分布），少数是异常值（分散分布）。这种参数配置下，仅通过观测值 X 很难判断样本来自哪个组分，特别是在中心区域，很难进行参数估计

2.2 观测值

在实际中，我们通常只能观测到变量 X 的取值，而无法直接观测到标签变量 L。这种情况下的参数估计问题变得更加复杂，因为标签的缺失使得似然函数无法简单地分解为独立项的乘积。这部分首先建立观测数据的似然函数，然后探讨在标签已知和未知两种情况下的推断方法。

问题14：建立观测值的似然函数

我们有 N 个观测值 \boldsymbol{x} = [x_1, \ldots, x_N]，未知参数向量为 \boldsymbol{\theta} = [\mu_1, \mu_2, \gamma_1, \gamma_2, p]。

对于单个观测 x_n，根据问题12的结果，其边缘密度为：

假设所有观测值独立同分布，似然函数为所有观测边缘密度的乘积：

每个观测都是两个高斯分布的加权和。

问题15：建立已知标签的似然函数

现在假设标签向量 \boldsymbol{l} = [l_1, \ldots, l_N] 已知，即对于每个观测 x_n，我们知道它来自哪个高斯分量。

对于单个观测对 (x_n, l_n)，联合概率密度为：

由于所有观测独立，完整的似然函数为：

定义以下索引集合：

• \mathcal{E}_1 = \{n : l_n = 1\}：标签为1的观测索引集合
• \mathcal{E}_2 = \{n : l_n = 2\}：标签为2的观测索引集合

显然 \mathcal{E}_1 \cap \mathcal{E}_2 = \emptyset 且 \mathcal{E}_1 \cup \mathcal{E}_2 = \{1, \ldots, N\}。

对于 n \in \mathcal{E}_1，有

对于 n \in \mathcal{E}_2，有

因此似然函数可以写成：

进一步分解：

其中 |\mathcal{E}_1| 和 |\mathcal{E}_2| 分别表示两个集合的基数。因此已知标签时，属于同一分量的观测可以分组处理，两组样本的似然函数可以独立计算，然后相乘得到总的似然函数，大简化了计算。

2.3 一个高斯分布的情况

在处理完整的双高斯混合模型之前，我们先考虑一个简化的问题：在只有一个高斯分布的情况下，在贝叶斯框架下推断均值 \mu_1 和精度 \gamma_1

贝叶斯框架下的参数推断

在贝叶斯框架下，我们为参数指定先验分布为共轭先验：均值 \mu_1 的先验为高斯分布，精度 \gamma_1 的先验为Gamma分布。

其中 \mu_{\mu_1} 和 \gamma_{\mu_1} 是均值先验的超参数，\alpha_1 和 \beta_1 是精度先验的超参数。在本实验中，我们使用相对无信息的先验，设置 \mu_{\mu_1}=0, \gamma_{\mu_1}=0.01, \alpha_1=1, \beta_1=1

问题18：推导后验分布

根据贝叶斯定理，后验分布正比于似然函数与先验分布的乘积：

在这里情况下，\boldsymbol{\theta} = [\mu_1, \gamma_1]，先验分布是可分离的：\mathrm{p}(\boldsymbol{\theta}) = \mathrm{p}(\mu_1) \mathrm{p}(\gamma_1)。因此后验分布为：

从问题17，我们知道似然函数为：

先验分布为：

将这三项相乘：

由于 \gamma_{\mu_1} 是先验超参数（固定常数），\gamma_{\mu_1}^{1/2} 项可以归入比例常数。合并 \gamma_1 的幂次和指数项：

这个后验分布结合了数据信息（似然项）和先验知识，其中指数项包含了关于 \mu_1 的二次型和关于 \gamma_1 的线性项

问题19：绘制后验分布的等高线

使用参数 \mu_1=-1, \gamma_1=1 生成 N=1000 个样本，在 (\mu_1, \gamma_1) 的二维空间中绘制后验概率密度的等高线图。

实验结果：

[此处插入图片9：后验分布的二维等高线图]

从等高线图可以看到，后验分布在真实参数值 (-1, 1) 附近达到最大值。

等高线的密集程度反映了后验分布的不确定性：在最大值附近等高线密集，表示后验分布集中；远离中心区域等高线稀疏，表示这些参数值的后验概率很低。

问题20：推导条件后验分布

a) 证明 \mu_1 的条件后验分布是高斯分布并确定其参数。

从问题18的后验分布出发，固定 \gamma_1，只保留与 \mu_1 相关的项：

展开指数项中的二次型，首先展开 \sum_{n=1}^N(x_n - \mu_1)^2：

再次展开 (\mu_1 - \mu_{\mu_1})^2 = \mu_1^2 - 2\mu_1\mu_{\mu_1} + \mu_{\mu_1}^2。将这些代入后验分布：

其中常数项（不含 \mu_1 的项）已被省略。

合并 \mu_1^2 和 \mu_1 的系数：

这是高斯分布的标准形式 \exp\left(-\frac{\gamma_{\text{post}}}{2}(\mu_1 - \mu_{\text{post}})^2\right)。通过配方法，我们可以识别出后验精度和后验均值：

因此：

这是一个高斯分布，其均值是数据信息（\sum_{n=1}^N x_n）和先验信息（\mu_{\mu_1}）的加权平均，权重分别为对应的精度。

b) 证明 \gamma_1​ 的条件后验分布是 Gamma 分布并确定其参数。

从问题18的后验分布出发，固定 \mu_1，只保留与 \gamma_1 相关的项：

这正是 Gamma 分布的标准形式 \mathcal{G}(\gamma; \alpha, \beta) \propto \gamma^{\alpha - 1} \exp(-\beta\gamma)。通过对比系数：

因此：

这是一个 Gamma 分布，其形状参数由数据量和先验共同决定，尺度参数包含了残差平方和与先验的贝塔参数。

问题21：实现Gibbs采样器

Gibbs采样器是一种特殊的马尔可夫链蒙特卡罗（MCMC）方法，它通过交替地从每个变量的条件后验分布中采样来生成联合后验分布的样本。对于我们的两参数问题，算法的每次迭代包含两步：首先固定当前的 \gamma_1^{(t-1)}，从 p(\mu_1|\boldsymbol{x}, \gamma_1^{(t-1)}) 中采样新的 \mu_1^{(t)}；然后固定刚采样的 \mu_1^{(t)}，从 p(\gamma_1|\boldsymbol{x}, \mu_1^{(t)}) 中采样新的 \gamma_1^{(t)}。

我们使用先验的模拟值作为初始化：\mu_1^{(0)} \sim \mathcal{N}(0, 100), \gamma_1^{(0)} \sim \mathcal{G}(1, 1)，然后运行 T=1000 次迭代

问题22：绘制轨迹

实验结果：

[此处插入图片10：Gibbs采样的轨迹图，包含mu_1和gamma_1的时间序列]

首先，两个参数的轨迹基本水平波动，这表明马尔可夫链已经收敛到平稳分布。其次，\mu_1 和 \gamma_1 的中心平均值不是真实值 1 ，有一些偏移，但偏移不大。

问题23：绘制直方图

参数边缘后验分布的经验估计直方图。

实验结果：

[此处插入图片11：mu_1和gamma_1的后验分布直方图，叠加理论真实值]

左图展示 \mu_1 的后验分布，呈现标准的单峰正态形态，后验均值非常接近真实值 -1.0。分布的形状之所以接近高斯分布，是因为在大样本下参数的后验分布趋于正态分布（中心极限定理）。

右图展示 \gamma_1 的后验分布，后验均值也非常接近真实值。

这些直方图对应各参数的后验边缘密度 \mathrm{p}(\mu_1 \mid \boldsymbol{x}) 和 \mathrm{p}(\gamma_1 \mid \boldsymbol{x})，通过1000个吉布斯样本估计得到。两个分布都是单峰的，验证了后验分布的良好性质，同时表明吉布斯采样的有效性

问题24：在平面上绘制二维轨迹

将Gibbs采样器生成的样本对 (\mu_1^{(t)}, \gamma_1^{(t)}) 在二维平面上绘制，并叠加到问题19的等高线图上，展示Gibbs采样如何探索后验分布。

实验结果：

[此处插入图片12：二维参数空间的轨迹图叠加在等高线图上]

从图中可以看到，Gibbs采样器经过多次迭代后，最终的终点非常接近真实值点。

问题25-26：计算后验统计量和区间概率

后验均值是参数在二次损失下的贝叶斯最优估计（MMSE估计），后验标准差量化了估计的不确定性。

实验结果：

表 1：后验统计量估计

表 2：后验区间概率

统计量结果显示，\mu_1 和 \gamma_1 的后验均值与真实值都非常接近。后验标准差也不大。

区间概率的计算结果显示，两个参数的后验均值加减一个标准差区间包含真实参数的概率约为68.00%。这个结果与正态分布的均值加减一个标准差区间结论非常接近，验证了在大样本下后验分布的近似正态性

3. 两个高斯分布的情况

在前面的实验中，我们处理了单一高斯分布的参数估计问题。现在我们考虑的混合模型场景：存在两个高斯分布，未知参数向量为 \boldsymbol{\theta}=[\mu_1, \mu_2, \gamma_1, \gamma_2, p]，并且标签 \boldsymbol{l}=[l_1, \ldots, l_N] 是隐藏的。因此我们只能观测到数据 \boldsymbol{x}，而不知道每个数据点来自哪个组分。

贝叶斯框架与先验选择

在贝叶斯框架下，我们将标签 \boldsymbol{l} 也视为随机变量，与参数 \boldsymbol{\theta} 一起进行联合推断。我们为每个参数选择先验分布：\mu_1, \mu_2 的先验为高斯分布，\gamma_1, \gamma_2 的先验为Gamma分布，混合参数 p 的先验为 [0,1] 上的均匀分布（Beta(1,1)分布）。均匀先验是一个无信息先验，表示我们对混合比例没有先验偏好。

问题27：建立后验分布

根据贝叶斯定理，后验分布正比于似然函数、标签的条件分布和参数先验的乘积：

又因为

所以

从问题15，已知标签时的完整似然函数为：

参数的先验分布为：

其中 \mu_1, \mu_2 的先验为高斯分布，\gamma_1, \gamma_2 的先验为 Gamma 分布，p 的先验为 [0,1] 上的均匀分布。将所有项相乘得到后验分布：

问题28：建立条件后验分布

\mu_1 的条件后验分布

固定 \boldsymbol{l} 和其他参数，只保留与 \mu_1 相关的项：

这与问题20a的推导完全类似，只是求和范围从全部 N 个观测变为 \mathcal{E}_1 中的观测。因此：

\gamma_1 的条件后验分布

固定 \boldsymbol{l} 和其他参数，只保留与 \gamma_1 相关的项：

合并 \gamma_1 的幂次和指数项：

\mu_2 的条件后验分布

\gamma_2 的条件后验分布

p 的条件后验分布

固定 \boldsymbol{l} 和其他参数，只保留与 p 相关的项：

这正是贝塔分布 \text{Beta}(\alpha, \beta) \propto p^{\alpha-1}(1-p)^{\beta-1} 的形式。对比系数得到：

综上所述，\mu_1, \mu_2 的条件后验是高斯分布，\gamma_1, \gamma_2 的条件后验是Gamma分布，p 的条件后验是Beta分布。

问题29：证明标签的条件后验分布

固定参数 \boldsymbol{\theta}，标签向量 \boldsymbol{l} 的条件后验分布为：

从问题15的似然函数，我们有：

由于观测独立，可以将其改写为：

因此：

这证明了条件后验分布是可分离的，每个 l_n 独立地服从自己的后验分布。对于单个 l_n，当 l_n = 1 时：

当 l_n = 2 时：

归一化后，l_n = 1 的后验概率为：

因此在给定参数和观测数据后，每个标签 l_n 服从参数为 P(l_n = 1 \mid x_n, \boldsymbol{\theta}) 的伯努利分布。它是两个组分在点 x_n 处的加权密度的比值，权重由先验混合比例 p 和 1-p 决定。

问题30：实现完整的Gibbs采样算法

完整的Gibbs采样算法循环更新六组变量：\mu_1, \gamma_1, \mu_2, \gamma_2, p, \boldsymbol{l}。在每次迭代中，我们按顺序从每个变量的条件后验分布中采样，固定其他所有变量。具体步骤如下：

第一步，固定 \gamma_1, \mu_2, \gamma_2, p, \boldsymbol{l}，从高斯分布 p(\mu_1|\boldsymbol{x}, \boldsymbol{l}, \gamma_1) 中采样新的 \mu_1。

第二步，固定 \mu_1, \mu_2, \gamma_2, p, \boldsymbol{l}，从Gamma分布 p(\gamma_1|\boldsymbol{x}, \boldsymbol{l}, \mu_1) 中采样新的 \gamma_1。

第三步和第四步类似地更新 \mu_2 和 \gamma_2。

第五步，固定所有参数和 \boldsymbol{l}，从Beta分布 p(p|\boldsymbol{l}) 中采样新的 p。

第六步，固定所有参数，对每个 n，从伯努利分布 p(l_n|x_n, \boldsymbol{\theta}) 中采样新的标签。

实验结果1：参数轨迹图

[此处插入图片13：五个参数的Gibbs采样轨迹图]

从五个参数的轨迹图可以看到，所有参数在经过初始的burn-in期后都稳定在真实值附近波动。\mu_1 的轨迹围绕 -2 波动，\mu_2 的轨迹围绕 3 波动，\gamma_1 的轨迹围绕 1 波动，\gamma_2 的轨迹围绕 2 波动，p 的轨迹围绕 0.6 波动。没有观察到明显的趋势或周期性，表明链已经收敛。各参数的波动幅度不同，反映了它们后验分布的不确定性差异。特别地，\gamma_2 的波动相对较大，这是因为样本量较小且方差参数本身的估计不确定性较大。

实验结果2：后验分布直方图

[此处插入图片14：五个参数的后验分布直方图]

五个参数的后验分布直方图都呈现出单峰且接近对称的形状，中心位置与真实值非常接近。这些直方图是基于9500个后验样本（去除500个burn-in样本）绘制的，它们提供了参数边缘后验分布的经验估计。所有参数的后验分布都相对集中，表明在给定1000个观测数据后，参数的不确定性已经大幅降低。真实参数值（用垂直线标记）都落在后验分布的高密度区域内，验证了估计的准确性。

实验结果3：二维后验分布

[此处插入图片15：关键参数对的二维散点图]

实验结果4：后验分布和真实分布对比

[此处插入图片16：真实标签与估计标签的混淆矩阵]

使用后验样本中每个 l_n 的众数（或后验概率超过0.5的标签）作为标签的点估计，然后与真实标签比较。结果显示，1000个样本中有约95%的标签被正确分类。错误分类主要发生在两个组分重叠的区域，即那些 x_n 值介于两个均值之间的样本。这个高准确率表明，即使标签是隐藏的，通过观测数据和Gibbs采样算法，我们仍然能够相当准确地推断每个样本的归属。

定量结果：

--- 后验统计量 ---
μ₁:  后验均值 -1.9483  (真实 -2.00),  后验标准差 0.0208
μ₂:  后验均值 2.9931  (真实 3.00),  后验标准差 0.0375
γ₁:  后验均值 0.9875  (真实 1.00),  后验标准差 0.0419
γ₂:  后验均值 1.9643  (真实 2.00),  后验标准差 0.1073
p:   后验均值 0.6053  (真实 0.60),  后验标准差 0.0108

--- 区间概率 ---
P(μ₁ ∈ [μ̂_μ₁ ± σ̂_μ₁]) = 0.6890
P(γ₁ ∈ [μ̂_γ₁ ± σ̂_γ₁]) = 0.6750

从定量结果可以看出，所有参数的后验均值都非常接近真实值，且后验标准差较小，表明对混合比例的估计非常确定。

区间概率的结果再次验证了后验分布的近似正态性（68%）。

4. 使用 Metropolis-Hastings 算法模拟 X|θ

算法动机与背景

在前面的实验中，我们使用层次构造方法从混合分布中采样：首先从伯努利分布采样标签 L，然后根据标签从相应的条件高斯分布采样 X。Metropolis-Hastings算法直接从边缘分布 p_X(x) 中采样，不需要显式地生成标签变量。

Metropolis-Hastings算法是马尔可夫链蒙特卡罗（MCMC）方法的一个重要分支，它提供了一个通用的框架来从复杂的概率分布中采样。与Gibbs采样不同(需要知道所有条件分布)，Metropolis-Hastings算法只需要能够计算目标分布的密度值（甚至不需要归一化），就可以生成服从该分布的样本。

问题31：用高斯随机游走实现Metropolis-Hastings采样

根据问题12的推导，混合分布的边缘密度为两个高斯分布的凸组合：

虽然每个组分是标准的高斯分布，但它们的线性组合使得整体分布变得难以直接采样，因此层次构造方法并不适用，我们使用Metropolis-Hastings算法。

Metropolis-Hastings算法的核心思想是构造一个马尔可夫链，其平稳分布恰好是我们的目标分布 p_X(x)。算法通过以下步骤实现：给定当前状态 x^{(n-1)}，从建议分布（proposal distribution）q(y|x^{(n-1)}) 中生成一个候选状态 y^{(n)}；然后计算接受概率，以该概率决定是接受候选状态（令 x^{(n)}=y^{(n)}）还是拒绝它（令 x^{(n)}=x^{(n-1)}，保持在原位置）。

接受概率的公式为：

这确保马尔可夫链满足细致平衡条件（detailed balance），从而保证目标分布是链的平稳分布。第一项 \frac{p_X(y^{(n)})}{p_X(x^{(n-1)})} 是目标密度的比值，它衡量候选点相对于当前点有多好；第二项 \frac{q(x^{(n-1)}|y^{(n)})}{q(y^{(n)}|x^{(n-1)})} 是建议分布的校正项，用于补偿建议分布的不对称性。

随机游走建议分布

在本实验中，我们采用最常用的建议分布：高斯随机游走。具体地，候选点定义为当前点加上一个高斯扰动：

这对应于建议分布 q(y|x^{(n-1)}) = \mathcal{N}(y; x^{(n-1)}, \sigma^2)。从 x 提议到 y 的概率密度等于从 y 提议到 x 的概率密度，即 q(y|x) = q(x|y)。因此，校正项的比值等于1，接受概率简化为：

因此Metropolis-Hastings被简化为Metropolis算法。如果候选点的密度更高，我们总是接受它（\alpha=1）；如果候选点的密度更低，我们以概率 p_X(y)/p_X(x) 接受它。因此算法倾向于向高密度区域移动，但也偶尔向低密度区域移动

算法结果

实验结果1：采样过程诊断

[此处插入图片18：Metropolis-Hastings采样诊断]

第一个子图展示了完整的采样轨迹。可以看到，轨迹从初始点0开始，在高密度区域之间随机游走，同时算法也会访问低密度区域，不会被困在单一的峰中

第二个子图展示了移动平均接受率随迭代的变化。接受率在0.65上下波动，在我们的目标区间中，因此说明算法在整个过程中保持稳定

实验结果2：后验分布对比

[此处插入图片19：MH样本的密度分布和累积分布函数对比]

左图展示 MH 样本的经验密度（蓝色直方图）与理论混合密度（红色线）非常接近，同时虚线分别表示两个组分的贡献，右图展示了经验CDF（红色）与理论CDF（蓝色）的两条曲线几乎重合，这说明了Metropolis-Hastings算法的有效性

定量验证结果：

样本统计量对比

可见MH样本均值和标准差和理论非常接近，这表明样本不仅捕捉到了分布的中心位置，也准确地反映了分布的离散程度

5. 线性同余生成器（可选）

在前面的所有实验中，我们都依赖于伪随机数生成器来产生均匀分布的随机数，这个部分我们将使用线性同余生成器（Linear Congruential Generator, LCG）来在均匀分布下生成样本。

原理

线性同余生成器基于一个简单的递推关系：给定前一个状态 X^{(n)}，下一个状态通过线性变换和取模运算得到

三个关键参数：乘数 a（multiplier）、增量 c（increment）和模数 m（modulus）。此外还需要一个初始值 X^{(0)}，称为种子（seed）。生成的整数序列 \{X^{(n)}\} 在集合 \{0, 1, \ldots, m-1\} 中取值，通过归一化 U^{(n)} = X^{(n)}/m 可以得到 [0,1) 区间上的伪随机数序列。

实验参数选择如下

实验结果：

[此处插入图片21：LCG、MATLAB rand()和理论分布的三个对比子图]

可见LCG生成数字的直方图和MATLAB内置 rand() 函数生成数字的直方图非常类似，都很接近理论均匀分布，其CDF线重合，表明两个生成器都准确地实现了均匀分布