时序分解 | Matlab实现贝叶斯变化点检测与时间序列分解

效果一览

在这里插入图片描述

基本介绍

Matlab实现贝叶斯变化点检测与时间序列分解
1.Matlab实现贝叶斯变化点检测与时间序列分解，完整源码和数据；
BEAST（突变、季节性和趋势的贝叶斯估计）是一种快速、通用的贝叶斯模型平均算法，用于将时间序列或1D序列数据分解为单个分量，如突变、趋势和周期性/季节性变化，如赵等人（2019）所述。BEAST可用于变化点检测（例如，断点、结构中断、状态变化或异常）、趋势分析、时间序列分解（例如，趋势与季节性）、时间序列分割和中断时间序列分析。
2.运行主程序main即可，其余为函数，无需运行，运行环境matlab2020及以上。
贝叶斯变化点检测和时间序列分解是两种在时间序列分析中常用的技术。

贝叶斯变化点检测（Bayesian Change Point Detection）是一种用于检测时间序列中突变点或结构变化的方法。它基于贝叶斯统计方法，通过考虑数据的先验分布和后验分布来确定变化点的位置和数量。该方法可以应用于多种类型的时间序列。
时间序列分解（Time Series Decomposition）是将时间序列分解为不同组成部分的过程。通常，一个时间序列可以分解为趋势（Trend）、季节性（Seasonality）和残差（Residual）三个部分。趋势表示时间序列的长期趋势变化，季节性表示时间序列在固定周期内的重复模式，而残差则表示无法由趋势和季节性解释的随机波动。时间序列分解可以帮助我们更好地理解时间序列的结构和特征，以及对序列进行预测和分析。

程序设计

完整源码和数据获取方式：Matlab实现贝叶斯变化点检测与时间序列分解。

%% get values from keys. The last arg is the default value if the key is missing from varagin/KeyList
  
   start           = GetValueByKey(KeyList, ValList, 'start',  []);
   deltat          = GetValueByKey(KeyList, ValList, 'deltat', []);
   time            = GetValueByKey(KeyList, ValList, 'time',   []);    
   period          = GetValueByKey(KeyList, ValList, 'period',  []); 
   nsamples_per_period  = GetValueByKey(KeyList, ValList, 'freq',  []); 
    
   season          = GetValueByKey(KeyList, ValList, 'season',        'harmonic'); 
   sorder_minmax   = GetValueByKey(KeyList, ValList, 'sorder.minmax', [1,5]); 
   scp_minmax      = GetValueByKey(KeyList, ValList, 'scp.minmax',    [0,10]); 
   sseg_min        = GetValueByKey(KeyList, ValList, 'sseg.min',      []); 
   sseg_leftmargin = GetValueByKey(KeyList, ValList, 'sseg.leftmargin',  []); 
   sseg_rightmargin= GetValueByKey(KeyList, ValList, 'sseg.rightmargin', []); 
   
   deseasonalize   = GetValueByKey(KeyList, ValList, 'deseasonalize', false); 
   detrend         = GetValueByKey(KeyList, ValList, 'detrend', false); 
   
   torder_minmax   = GetValueByKey(KeyList, ValList, 'torder.minmax', [0,1]); 
   tcp_minmax      = GetValueByKey(KeyList, ValList, 'tcp.minmax',    [0,10]); 
   tseg_min        = GetValueByKey(KeyList, ValList, 'tseg.min',      []);
   tseg_leftmargin = GetValueByKey(KeyList, ValList, 'tseg.leftmargin',  []); 
   tseg_rightmargin= GetValueByKey(KeyList, ValList, 'tseg.rightmargin', []); 

   precValue       = GetValueByKey(KeyList, ValList, 'precValue',       1.5); 
   precPriorType   = GetValueByKey(KeyList, ValList, 'precPriorType',   'componentwise');    
   
   hasOutlierCmpnt = GetValueByKey(KeyList, ValList, 'hasOutlier',        []); 
   ocp_max         = GetValueByKey(KeyList, ValList, 'ocp.max',           10); 
      
   mcmc_seed       = GetValueByKey(KeyList, ValList, 'mcmc.seed',     0);         
   mcmc_samples    = GetValueByKey(KeyList, ValList, 'mcmc.samples',  8000);
   mcmc_thin       = GetValueByKey(KeyList, ValList, 'mcmc.thin',     5); 
   mcmc_burnin     = GetValueByKey(KeyList, ValList, 'mcmc.burnin',   200);
   mcmc_chainNumber= GetValueByKey(KeyList, ValList, 'mcmc.chains',   3);  
   
   ci               = GetValueByKey(KeyList, ValList, 'ci',             false);   
   printProgressBar = GetValueByKey(KeyList, ValList, 'print.progress', true);     
   printOptions     = GetValueByKey(KeyList, ValList, 'print.options',  true);    
   quiet            = GetValueByKey(KeyList, ValList, 'quiet',          false);   
   gui              = GetValueByKey(KeyList, ValList, 'gui',            false); 
   methods          = GetValueByKey(KeyList, ValList, 'method',        'bayes'); 
   
%% Convert the opt parameters to the individual option parameters (e.g., metadata, prior, mcmc, and extra)

   %......Start of displaying 'MetaData' ......
   metadata = [];
   metadata.isRegularOrdered = true;
   metadata.season           = season;
   metadata.time             = time;
   metadata.startTime        = start;
   metadata.deltaTime        = deltat;
   if isempty(period) && ~isempty(deltat) && ~isempty(nsamples_per_period) && ~strcmp(season, 'none')
       period=nsamples_per_period*deltat;
   end   
   metadata.period           = period;

 
   if strcmp(metadata.season, 'svd')
      % if isempty(freq)|| freq <= 1.1 || isnan(freq)
      %     error("When season=svd, freq must be specified and larger than 1.");
      % end
      % metadata.svdTerms = svdbasis(y, freq, deseasonalize);
   end
   metadata.missingValue     = NaN;
   metadata.maxMissingRate   = 0.75;
   metadata.deseasonalize    = deseasonalize;
   metadata.detrend          = detrend;
   metadata.hasOutlierCmpnt  = hasOutlierCmpnt;
%........End of displaying MetaData ........

%......Start of displaying 'prior' ......
   prior = [];
   prior.modelPriorType	  = 1;
   if ~strcmp(metadata.season, 'none')              
       prior.seasonMinOrder   = sorder_minmax(1);
       prior.seasonMaxOrder   = sorder_minmax(2);
       prior.seasonMinKnotNum = scp_minmax(1);
       prior.seasonMaxKnotNum = scp_minmax(2);   
       prior.seasonMinSepDist = sseg_min;
	   prior.seasonLeftMargin  = sseg_leftmargin;
	   prior.seasonRightMargin = sseg_rightmargin;
   end   
   prior.trendMinOrder	  = torder_minmax(1);
   prior.trendMaxOrder	  = torder_minmax(2);
   prior.trendMinKnotNum  = tcp_minmax(1);
   prior.trendMaxKnotNum  = tcp_minmax(2);
   prior.trendMinSepDist  = tseg_min;
   prior.trendLeftMargin  = tseg_leftmargin;
   prior.trendRightMargin = tseg_rightmargin;

   if hasOutlierCmpnt
      prior.outlierMaxKnotNum = ocp_max;
   end
        
   prior.precValue        = precValue;
   prior.precPriorType    = precPriorType;
%......End of displaying pripr ......

%......Start of displaying 'mcmc' ......
   mcmc = [];
   mcmc.seed                      = mcmc_seed;
   mcmc.samples                   = mcmc_samples;
   mcmc.thinningFactor            = mcmc_thin;
   mcmc.burnin                    = mcmc_burnin;
   mcmc.chainNumber               = mcmc_chainNumber;
   
   %mcmc.maxMoveStepSize           = 28
   mcmc.trendResamplingOrderProb  = 0.1000;
   mcmc.seasonResamplingOrderProb = 0.1700;
   mcmc.credIntervalAlphaLevel    = 0.950;
%......End of displaying mcmc ......

%......Start of displaying 'extra' ......
   extra = [];
   extra.dumpInputData        = true;
   extra.whichOutputDimIsTime = 1;
   extra.computeCredible      = ci;
   extra.fastCIComputation    = true;
   extra.computeSeasonOrder   = true;
   extra.computeTrendOrder    = true;
   extra.computeSeasonChngpt  = true;
   extra.computeTrendChngpt   = true;
   extra.computeSeasonAmp     = ~strcmp(metadata.season, 'svd');
   extra.computeTrendSlope    = true;
   extra.tallyPosNegSeasonJump= false;
   extra.tallyPosNegTrendJump = false;
   extra.tallyIncDecTrendJump = false;
   extra.printProgressBar     = printProgressBar;
   extra.printOptions         = printOptions;
   extra.quiet                = quiet;
   extra.consoleWidth         = 70;
   extra.numThreadsPerCPU     = 2;
   extra.numParThreads        = 0;
%......End of displaying extra ......


 if (gui)
    out=Rbeast(' beastv4demo',            y, metadata, prior, mcmc, extra);
 else
    out=Rbeast( strcat('beast_',methods), y, metadata, prior, mcmc, extra);
 end
 
end