简介：这是一套基于BEAST算法的贝叶斯集成变点检测与时间序列分解实现，面向Matlab用户，尤其适合需要完成课程设计、期末大作业或毕业设计的电子信息、计算机、数学等专业学生。代码采用参数化编程，关键参数可方便调整，注释清晰，便于二次开发与理解算法原理。资源共266个文件，压缩包约6.14MB，以C语言核心源文件（.c/.h）、Matlab脚本（.m）、Python辅助脚本（.py）和MAT数据文件（.mat）为主，另有结果图（.png）与少量说明文档，整体结构覆盖算法实现、调用示例与运行结果。目前已吸引437人学习浏览，内容包含可直接运行的案例数据，在Matlab 2014/2019a/2021a下均可运行，下载后可按需替换数据或调节参数，快速复现变点检测与趋势分解效果，是理解BEAST算法和开展相关实验的高效参考资料。
1. 变点检测不是找“拐点”，而是给时序做贝叶斯分解
接触 BEAST 之前，我处理变点检测最常用的工具是findchangepts，它能在给定最大分段数的前提下把序列劈成几段。但真实业务里的时间序列往往同时带着趋势漂移、季节波动和多个突变，单纯二分一段并不能回答“季节分量是否也在变”这类问题。BEAST 的做法是把它变成一个贝叶斯集成问题：把序列拆成趋势项加季节项加余项，再用 MCMC 采样大量“变点数量与位置都不同”的模型，最后把所有模型对每个时间点的变点概率做加权平均。也就是说，它输出的不是一条唯一的分割线，而是一条概率曲线。
这次压缩包里的 Matlab 代码保留了完整的 C 核心，包括beastv2_COREV4.c、多线程版beastv2_COREV4_mthrd.c，以及abc_vec_avx2.c、abc_vec_avx512.c这类向量化辅助源文件，还带了可直接运行的案例数据。Matlab 2014、2019a、2021a 下都能跑。适合做信号故障定位、传感器数据分段、水文气象序列突变分析，也适合拿来当课程设计和期末大作业的算法底座。
2. 贝叶斯集成与 C 核心：BEAST 的模型平均化如何落到代码上
BEAST 的英文全称是 Bayesian Estimator of Abrupt change, Seasonal change, and Trend。它把时间序列建模成趋势分量、季节分量和余项，趋势与季节分量内部都允许出现突变。与传统“先选模型再估计参数”的两阶段思路不同，BEAST 直接对“模型集合”做后验采样，因此对变点数量不确定、位置不确定、季节谐波阶数不确定这些问题，都能给出概率意义上的答案。
2.1 集成：把“某一个模型”换成“一堆模型的后验加权”
传统变点检测里，常见流程是枚举若干个分段数，对每个分段数估计参数，再用 BIC 或 AIC 挑一个“最优”模型。这个做法在变点间隔较短、噪声较大时很不稳定，因为两个相邻候选模型的后验差异可能非常小，选任何一个都会把不确定性丢掉。
BEAST 的集成逻辑可以近似写成下面这段概念性 Matlab 代码：
% 概念演示：MCMC 每步都会得到一个“模型模式” % model_pattern(m, :) 记录第 m 个采样模型在每个时刻是否为变点 % model_weight(m) 是该模型对后验的贡献权重 for t = 1:T tcp_prob(t) = sum( model_weight(:) .* model_pattern(:, t) ); end这里的关键是：tcp_prob(t)不是某个最优模型给出的点估计，而是大量 MCMC 样本在 t 时刻是变点的平均支持度。每个采样模型有自己的变点数量、变点位置、趋势斜率、季节谐波阶数，最后把所有可能性按照后验权重合并。这样，当数据只支持“这里可能有变点，但位置模糊”时，概率曲线会摊开成一段宽峰，而不是给出一个虚假的精确位置。
这套逻辑在 C 代码里并不好读，因为需要对每个候选变点位置计算条件后验并做累加，但算法思路和上面的累加式是一致的。实际使用中，我一般不会去改 C 主体，而是通过 Matlab 参数控制 MCMC 采样规模。
2.2 源码文件：哪些是算法，哪些是数据
解压包后最先看到的是几个同名但后缀不同的 C 文件，很容易让人迷惑。我的建议是先按下面这张表把文件分类：
| 文件 | 角色 |
|---|---|
beastv2_COREV4.c | 单线程核心计算 |
beastv2_COREV4_mthrd.c | 多线程版，长序列优先用它 |
beastv2_COREV4_gui.c | 界面调用入口相关 |
abc_vec_avx2.c | AVX2 向量化辅助计算 |
abc_vec_avx512.c | AVX512 向量化辅助计算 |
plotbeast.asv | Matlab 自动保存文件，可忽略 |
mthrd后缀意味着内部用了多线程，适合在核数较多的机器上跑长序列。avx2和avx512是 CPU 指令集层面的优化，AVX512 吞吐更高，但老 CPU 不一定支持。常见的做法是先跑mex -setup确认编译器，再针对当前机器编译对应版本：
# 先确认 C 编译器 mex -setup C # 编多线程版，替代默认单线程 .mex 文件 mex -v beastv2_COREV4_mthrd.c这一步不是必须的。如果压缩包里已经带了编译好的mex文件，并且你的 Matlab 版本能直接加载，上面的命令可以跳过。plotbeast.asv是编辑器自动保存的临时文件，加路径时不用管它。
3. Matlab 参数化调用 BEAST：把包挂进路径，用一组参数跑完整分解
拿到这类源码包，第一步不是急着改算法，而是先把目录挂进 Matlab 工作区，确认函数能被识别。这套代码的调用方式是参数化的，核心参数都可以在入口处调整，不用深入 C 内部。
3.1 把 BEAST 挂进 Matlab 路径
我习惯把解压后的目录放到一个不含中文的路径下，然后统一用genpath添加所有子目录：
% 放到 D:/work/BEAST 后执行 addpath(genpath('D:/work/BEAST')); savepath; % 验证是否识别到核心函数 which beast which plotbeastaddpath(genpath(...))会把根目录和所有子目录加入搜索路径，因为有些版本会把beast和plotbeast放在src、mfiles等子目录里。savepath是把当前路径保存到 Matlab 的 pathdef 文件中，否则下次重启还要重新 addpath。最后用which确认函数存在，如果返回空，说明路径没加对或者目录名被截断。
3.2 核心参数表
BEAST 的参数风格类似属性键值对，不需要把全部参数传进去，系统会对没给的部分使用默认值。我最常改的是下面这几个：
| 参数 | 作用 |
|---|---|
y | 待分解的观测序列，建议用列向量 |
period | 季节周期长度，比如月数据为 12，日数据为 365 |
start | 起始时间标签，用于绘图横轴 |
tcp.minmax | 趋势变点数的最小值和最大值 |
scp.minmax | 季节变点数的最小值和最大值 |
sorder.minmax | 季节谐波阶数范围 |
mcmc.samples | MCMC 采样量，越大越稳定但越慢 |
mcmc.burnin | 预烧期长度，用于丢弃初值噪声 |
tcp.minmax和scp.minmax是这套算法最有意思的参数。它并不是让你指定精确变点数，而是给先验范围。范围越宽，模型空间越大，计算越慢。如果你只关心趋势突变，可以把scp.minmax设成[0, 0]，强制季节项没有变点。
3.3 跑一个最小例子
压缩包里案例数据的格式通常是一列日期、一列观测值。读取时要注意 Matlab 2014 没有readmatrix，老版本可以用csvread或xlsread。我这里按新版本写，老版本替换即可：
data = readmatrix('case_data.csv'); % 第一列日期，第二列观测值 y = data(:, 2); % 取观测序列 x = data(:, 1); % 取日期数值 out = beast(y, ... 'period', 12, ... 'start', x(1), ... 'tcp.minmax', [0, 6], ... 'scp.minmax', [0, 2], ... 'sorder.minmax', [1, 4]); plotbeast(out);这段代码的核心是把y作为第一参数，后面全部用键值对传参。start只在绘图时影响横轴刻度，不影响分解结果。tcp.minmax设成[0, 6]表示趋势项最多允许出现 6 个突变点；如果你不确定，可以把上限放大到 10，但计算时间会明显增加。
跑完之后，out里会带趋势后验均值、季节后验均值、余项、趋势变点概率、季节变点概率等字段。不同 fork 的字段名不完全一样，稳妥做法是先执行fieldnames(out)看一下再取数。
4. 用带突变与季节项的合成数据验证 BEAST 的变点概率
给别人讲 BEAST 的时候，我最喜欢用的验证方式不是直接上真实数据，而是构造一条“答案已知”的序列：在某个位置放一个突变，再加季节项和高斯噪声。这样就能直接看出变点概率峰是否出现在真实突变点附近。
4.1 构造一个有明确答案的测试序列
构造逻辑是三年日尺度数据，前 249 个点是均值为 10、斜率为 0.02 的线性趋势，第 250 个点开始整体抬升 5 个单位，并增加一个斜率变化，再叠加两个谐波的季节项和噪声：
rng(42); % 固定随机种子，保证结果可复现 n = 365 * 3; t = (1:n)'; trend = 10 + 0.02 * t; trend(250:end) = trend(249) + 5 + 0.02 * (t(250:end) - t(249)); season = 3 * sin(2 * pi * t / 365) + 1.5 * sin(4 * pi * t / 365 + 1); y = trend + season + randn(n, 1);这里rng(42)很关键，没有固定随机种子，每次跑出来的噪声不同，验证结论会漂移。trend(250:end)这一句在trend(249)的基础上加 5，相当于人为制造一个阶跃突变，同时保留原斜率，使问题干净一些。季节项用了两个谐波，对应 BEAST 里sorder.minmax能覆盖的情况。
4.2 跑 BEAST 并检查三个输出
合成数据的周期是 365，所以period设为 365。为了看得更清楚，MCMC 采样量设成 5000：
out = beast(y, ... 'period', 365, ... 'tcp.minmax', [0, 3], ... 'scp.minmax', [0, 2], ... 'mcmc.samples', 5000); % 趋势变点概率的字段在不同版本里可能叫 out.tcp 或 out.tcp.prob if isfield(out, 'tcp') && isstruct(out.tcp) tcp_prob = out.tcp.prob; else tcp_prob = out.tcp; end figure('Position', [100 100 800 700]); subplot(3, 1, 1); plot(t, y, 'Color', [0.7 0.7 0.7]); hold on; plot(t, out.trend, 'r', 'LineWidth', 1.5); title('原始序列与趋势后验均值'); subplot(3, 1, 2); plot(t, out.season, 'b'); title('季节分量'); subplot(3, 1, 3); plot(t, tcp_prob, 'r'); xlabel('时间'); ylabel('变点概率'); title('趋势变点概率');out.trend是趋势分量的后验均值，out.season是季节分量的后验均值，两者都应该是平滑的估计结果。第三幅图里，概率峰应该集中在 t=250 附近。如果曲线在 250 处出现一个明显尖峰，说明算法把人工突变位置成功找回来了。
4.3 概率阈值：不要只认 0.5
很多第一次用的人会把变点判成“概率大于 0.5 的位置”，这个习惯在信噪比高的时候没问题，但在真实业务里经常漏报。BEAST 输出的概率是“当前数据支持这里有变点的强度”，并不是一个严格的检验 p 值。小样本或突变幅度不大时，真实变点位置的概率可能只有 0.2 到 0.4。
我一般会先看概率曲线的整体形状，再找局部峰值，而不是硬切固定阈值。如果序列的突变幅度足够大，峰会很尖；如果模糊，峰会很宽。峰宽本身就是不确定性信息，这一步比单纯阈值判断更有价值。
5. 长序列优化：MEX 多线程、谐波阶数与常见参数坑
BEAST 的计算量主要集中在 MCMC 对整个模型空间进行采样。序列长度到几万点、变点数量上限放宽之后，纯 Matlab 主循环基本跑不动。这套代码把核心计算放在 C 层，目的就是让长序列优化有了实施空间。
5.1 用mthrd还是avx512
如果你的机器是 8 核以上，优先用beastv2_COREV4_mthrd.c。如果 CPU 比较新且支持 AVX512，可以把向量化辅助代码一起编进去。Linux 下常见做法是加 OpenMP 编译选项：
% 在 Matlab 命令行里执行，不要写在脚本里反复编译 mex -setup C mex -v CFLAGS="$CFLAGS -fopenmp -O2 -mavx2" beastv2_COREV4_mthrd.cWindows 下的写法略有不同，得先装 MinGW-w64 或 Microsoft C/C++ 编译器，再把-fopenmp换成对应编译器支持的选项。没有把握时不要强行编译，直接用包内已编译好的.mex文件更稳妥。
参数与速度的取舍同样重要。sorder.minmax从[1, 4]改成[1, 8]，意味着季节模型空间明显扩大，MCMC 需要更多样本来覆盖。如果数据只是年周期，阶数给 4 通常够了。mcmc.samples也不要盲目调到五万，先跑三千看稳定性，再逐步加大。
5.2 常见参数坑：NaN、周期和日期列
第一类坑是序列里带了 NaN。BEAST 的 MCMC 链路要求输入是连续完整的序列，有空洞的位置概率计算会出错。常见做法是提前用线性插值补齐：
if any(isnan(y)) y = fillmissing(y, 'linear'); end第二类坑是period给错。有人拿日数据直接填 12，结果季节项变成“每 12 天一个周期”，分解结果会非常奇怪。判断标准很简单：看最小季节回环的步长。日数据有年周期就用 365，周周期就用 7；月数据用 12；周数据用 52。
第三类坑是读入数据时把日期列当成观测列。readmatrix读出的矩阵如果有三列以上，要确认y到底取的是哪一列。我经常在调用前先画一条原始序列，肉眼确认后再跑算法，能省掉很多定位时间。
6. 把 BEAST 的变点概率自动翻译成断点列表：一个可复用函数
BEAST 给的是概率曲线，但业务系统通常要的是“哪些时间点是断点”。直接用一个固定阈值去挑，在小信号场景下会漏掉真实变点。我会写一个简单的局部峰值提取函数，只保留时间间隔足够远的峰，再把概率超过某个分位数的点当成候选断点。
function cp = find_beast_changepoints(tcp_prob, min_gap, pct) if nargin < 2 || isempty(min_gap) min_gap = 5; end if nargin < 3 || isempty(pct) pct = 90; end tmp = tcp_prob(:); thr = max(0.05, prctile(tmp, pct)); if exist('findpeaks', 'file') == 2 [~, locs] = findpeaks(tmp, ... 'MinPeakHeight', thr, ... 'MinPeakDistance', min_gap); cp = locs; else cp = find(tmp(2:end-1) > tmp(1:end-2) & ... tmp(2:end-1) >= tmp(3:end)) + 1; if ~isempty(cp) keep = [true; diff(cp) >= min_gap]; cp = cp(keep); end end end这个函数有两个关键参数。min_gap是相邻候选断点的最小时间间隔，避免同一个“宽峰”被拆成多个假变点；pct是概率分位数，默认取 90% 分位数作为阈值起点。对幅值明显的大突变，可以降到 80% 或直接看峰形；对噪声大的短序列，则提高分位数，减少误报。
调用方式也很直接：先用上一章的方式取出趋势变点概率，再传入函数得到断点位置：
cp = find_beast_changepoints(tcp_prob, 30, 90); plot(t, tcp_prob); hold on; plot(t(cp), tcp_prob(cp), 'ro');把min_gap设为周期的一半通常比较稳，例如日数据周期 365 时，不同变点之间至少隔 180 天左右。这个函数不依赖 Signal Processing Toolbox，找不到findpeaks时会自动回退到简单的一阶差分峰值判断，适合放进课程设计或工程脚本里直接复用。
本文还有配套的精品资源，点击获取
