基于马尔可夫链蒙特卡洛(MCMC)原理的复杂分布数值模拟系统
项目介绍
本项目是一个基于MATLAB开发的综合性数值模拟平台,旨在通过马尔可夫链蒙特卡洛(MCMC)方法解决统计学与数值分析中难以直接解析的概率密度函数抽样问题。系统通过构建平稳分布为目标分布的随机过程,实现了在多维状态空间中的随机游走与样本采集。该系统广泛应用于贝叶斯推断、复杂物理系统模拟及高维积分计算,能够帮助研究人员准确获取目标分布的统计特性及演化行为。
功能特性
- 双核心算法支持:内置主流的Metropolis-Hastings算法与Gibbs采样算法,满足不同场景下的抽样需求。
- 高维空间处理:支持多维参数空间的联合分布演化,能够处理具有相关性的变量抽样。
- 完善的收敛性诊断:提供燃尽期(Burn-in)设定、参数演化轨迹分析(Trace Plot)及自相关性分析(ACF),确保样本的独立性与可靠性。
- 可视化统计分析:自动生成后验分布直方图、二维联合空间散点图以及详细的数值评估报告。
- 灵活的参数配置:用户可自定义迭代步数、步长、初始位置及目标分布函数。
系统要求- 软件环境:MATLAB R2016b 或更高版本。
- 依赖工具箱:Statistics and Machine Learning Toolbox(用于调用目标密度评价函数)。
功能实现逻辑说明系统的核心逻辑遵循 MCMC 的基本流程,通过以下步骤实现数值模拟:
- 环境初始化与参数预设:程序启动后,首先初始化状态空间维度、迭代总次数(如10000次)以及燃尽期长度(如2000次)。同时,定义目标二元高斯分布的均值向量和协方差矩阵。
- 状态转移与样本采集:
- 系统运行Metropolis-Hastings逻辑:从初始点出发,利用随机游走产生建议样本,根据目标概率密度之比计算接受率,决定是否转移到新状态。
- 系统运行Gibbs采样逻辑:利用解析的条件概率分布,交替抽取各个维度的条件样本,从而构建完整的状态向量。
- 数据预处理:自动剔除燃尽期内的初始样本,消除初值偏差对比拟合分布的影响,保留后续平稳阶段的有效样本。
- 统计量估算:对有效样本进行均值、方差及接受率(MH算法专用)的数值计算。
- 多维度可视化:将计算结果转化为六联图表,包括不同参数的时间序列轨迹、边缘分布直方图、样本自相关系数及联合采样空间分布。
关键核心算法分析
- Metropolis-Hastings 逻辑实现:
该部分通过建议分布(以当前点为中心的正态分布)进行随机搜索。核心在于接受概率 alpha 的计算,即目标 PDF 在新旧位置的比例。程序通过生成的均匀随机数与 alpha 比较,实现了对高能态区域的重点采样。
针对二元正态分布,该算法通过解析推导出的条件分布公式进行采样。程序计算了两个变量间的相关系数 rho 以及条件均值与方差,实现了在不显式计算接受概率的情况下,通过维度间的迭代更新逼近联合分布。
系统内置了手动实现的自相关函数计算逻辑。通过对有效样本序列进行滞后相关性分析,计算从0阶到100阶的滞后系数。如果自相关系数随滞后阶数迅速下降,则表明链的混合程度良好,样本独立性高。
轨迹图通过绘制迭代步数与对应取值的关系,展示了马尔可夫链的遍历性。联合分布散点图则通过设置透明度标记(MarkerFaceAlpha),直观地展示了高概率区域的样本聚集情况。
使用方法
- 启动系统:在MATLAB编辑器中打开代码并运行。
- 配置参数:在代码顶部的参数设置区域,根据模拟需求修改总迭代步数(total_samples)、燃尽期(burn_in)以及步长(step_size)。
- 定义目标分布:在目标密度函数句柄处修改对应的概率密度公式。
- 查看报告:运行结束后,系统将在命令行窗口输出“MCMC 模拟报告”,包含MH算法的接受率、估计均值与真实值的偏差等数值。
- 分析图表:系统会自动弹出可视化界面,用户可通过分析自相关图确定样本的有效性,通过轨迹图确认燃尽期设置是否合理。