含再采样机制的RJMCMC算法实现与验证系统

资 源 简 介

本项目旨在MATLAB环境下开发一套增强型的可逆跳跃马尔可夫链蒙特卡洛（RJMCMC）算法，核心在于引入再采样（Resampling）步骤以解决传统跨维采样算法在复杂多模态分布中混合效率低和容易陷入局部最优的问题。项目详细功能包括：1. 基础RJMCMC框架构建，实现参数空间维度的动态变化，包括“出生”、“死亡”、“分裂”和“合并”等可逆移动操作，并精确计算Metropolis-Hastings-Green接受率及涉及的雅可比行列式；2. 再采样机制集成，借鉴序列蒙特卡洛（SMC）思想，在迭代过程中根据粒子（或链的状态）的权重进行重采样，优胜劣汰，有效防止样本退化并促进对高概率区域的探索；3. 典型应用算例设计，采用变阶多项式曲线拟合或未知分量数的高斯混合模型（GMM）作为测试对象，生成包含噪声的仿真数据；4. 性能评估与可视化，通过对比标准RJMCMC与带有再采样步骤的RJMCMC算法，输出模型阶数的后验直方图、参数估计的收敛轨迹、数据拟合效果图以及算法的接受率，从而量化证明改进算法在模型识别准确度、收敛速度和稳健性方面的性能优势。

详 情 说 明

带有再采样机制的RJMCMC算法实现与性能验证系统

项目简介

本项目是一个基于MATLAB开发的高级统计推断系统，旨在解决变维模型选择和参数估计问题。通过结合可逆跳跃马尔可夫链蒙特卡洛（RJMCMC）与序列蒙特卡洛（SMC）中的再采样机制，本系统实现了一个针对多项式曲线拟合的增强型贝叶斯推断框架。

该算法能够在参数空间维度未知的情况下（即多项式阶数未知），同时估计最佳模型阶数、回归系数以及噪声方差。引入的再采样步骤有效地改善了传统MCMC算法在多模态分布中容易陷入局部最优和混合效率低下的问题。

主要功能特性

• 变维贝叶斯推断：支持在不同维度的参数空间之间进行自动跳转，能够自动识别数据的潜在模型复杂度（多项式阶数）。
• 混合采样策略：

* Metropolis-Hastings：用于回归系数（Beta）的随机游走更新。 * Gibbs采样：利用共轭先验特性，对噪声方差（Sigma2）进行直接采样，提高了收敛效率。

• RJMCMC移动操作：实现了三种核心移动类型：

* Update（参数更新）：保持模型阶数不变，更新模型参数。 * Birth（维度增加）：增加一个模型阶数，并提议新的系数。 * Death（维度减少）：减少一个模型阶数，移除最高阶系数。

• 粒子再采样增强：借鉴粒子滤波思想，定期根据粒子的似然权重进行系统重采样（Systematic Resampling），优胜劣汰，防止低概率样本的无效计算，提升对高概率区域的探索能力。
• 自动化可视化分析：包含完整的后验分析模块，自动绘制阶数直方图、MCMC轨迹及拟合结果。

系统要求及使用方法

系统要求

• MATLAB R2016a 或更高版本。
• 不需要额外的工具箱（代码仅使用基础数学和统计函数）。

使用方法

直接运行主脚本即可启动整个流程。系统将自动执行以下步骤：

1. 清理环境并固定随机种子以保证实验可复现。
2. 生成基于3阶多项式的仿真数据（包含高斯噪声）。
3. 初始化粒子群并开始迭代采样。
4. 在控制台实时输出迭代进度和当前平均模型阶数。
5. 迭代结束后，自动弹出图形窗口展示分析结果。

核心算法实现细节

本项目的实现严格遵循代码逻辑，主要包含以下几个关键部分：

1. 仿真数据生成

系统首先构建一个真实的非线性环境。代码在[-1, 1]区间内生成均匀分布的随机变量X，并通过一个设定的3阶多项式（$y = 1 + 2x - 1.5x^2 + 3x^3$）计算基准值，最后叠加标准差为0.5的高斯白噪声生成观测数据Y。

2. 初始化与先验设定

• 粒子群初始化：创建包含20个粒子的群体，每个粒子独立随机初始化0-2之间的阶数，并将系数初始化为零。
• 先验分布：

* 模型阶数 ($k$)：服从泊松分布，用于惩罚过高的模型复杂度。 * 回归系数 ($beta$)：服从独立的正态分布。 * 噪声方差 ($sigma^2$)：服从逆伽马分布（Inverse-Gamma），作为正态似然的共轭先验。

3. 主迭代循环 (混合RJMCMC)

算法运行5000次迭代，每次迭代遍历所有粒子并执行以下逻辑：

#### 移动类型选择 系统根据设定的概率向量 [0.4, 0.3, 0.3] 随机选择操作类型，并针对边界条件（如阶数为0时不能降维，达到最大阶数时不能升维）进行了概率归一化处理。

#### 具体移动实现

• Update (维度不变)：

* 对系数向量 $beta$ 的每个元素进行小步长的随机游走，计算Metropolis接受率。 * 对噪声方差 $sigma^2$ 使用Gibbs采样，通过残差平方和计算后验参数，直接从逆伽马分布中抽取新值。这确保了方差估计的高效性。

• Birth (升维 $k to k+1$)：

* 提议新的阶数，并从正态分布 $N(0, 0.5)$ 中抽取一个新的高阶系数。 * 构建更复杂的模型，计算似然比、先验比和提议比（Proposal Ratio）。 * 根据计算出的Metropolis-Hastings-Green接受率决定是否接受新维度。

• Death (降维 $k to k-1$)：

* 提议移除当前的最高阶系数。 * 计算相应的接受率，如果接受则模型退化为低阶多项式。

4. 再采样机制 (Resampling)

为了克服传统MCMC链容易陷入局部极值的问题，代码集成了再采样步骤：

• 频率：每50次迭代执行一次。
• 权重计算：基于当前粒子的对数似然值计算非归一化权重，并使用 exp(log_w - max_log_w) 技巧防止数值溢出。
• 系统重采样：根据累积权重分布（CDF），保留高似然粒子，淘汰低似然粒子。这使得多个马尔可夫链能够交互信息，集中算力探索高概率模型空间。

5. 结果处理与可视化

process_results 函数负责将仿真结果转化为直观的图表：

• 模型阶数后验分布：统计预烧期（Burn-in）后的采样数据，绘制直方图，并标记真实阶数（True k），用于验证算法是否正确识别了模型结构。
• 采样轨迹图 (Trace Plot)：绘制第一个粒子的模型阶数随迭代次数的变化，用于评估马尔可夫链的混合性能和收敛状态。
• 多项式拟合：根据后验样本计算统计特性，展示拟合效果（基于代码逻辑推断）。

辅助函数说明

• calculate_log_likelihood：核心评估函数。计算给定模型参数下的高斯对数似然值。如果方差非正，则返回负无穷。
• calculate_log_prior：计算当前粒子状态的联合对数先验概率，包含阶数（泊松）、系数（高斯）和方差（逆伽马）三部分。
• select_move：根据当前阶数状态和预设概率，智能选择下一步的移动类型，处理了0阶和最大阶数的边界逻辑。
• get_prob_move：计算从当前状态选择特定移动的先验概率，这是计算可逆跳跃接受率中Proposal Ratio的关键组成部分，确保了详细平衡条件。