基于有序模式的度量对多变量时间序列进行非线性分析（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客 

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

OPA（有序模式分析）工具箱旨在对多变量时间序列进行非线性分析，随着越来越流行的基于顺序模式的度量[1-5]，这些度量可以有效地计算[6，7]并可视化：
- 排列熵 （cfg.method = 'PE'） [2] - 具有平秩的有序模式的排列熵 （cfg.method = 'eqPE'） [4，8] - 排列熵和有序模式分布 （cfg.method = 'opdPE'） [3]

- 有序模式的条件熵 （cfg.method = 'cePE'） [6] - 鲁棒排列熵 （cfg.Method = 'rePE'） [4，7]

- 将添加新的基于有序模式的措施（欢迎您的建议和反馈！

输入
OPA 工具箱的接口由函数 outdata = OPanalysis（ cfg， indata ） 提供，其中
- cfg 是具有方法参数的配置结构;
- INDATA 是要分析的数据。

输出
- 输出数据 - 基于序数模式的复杂性度量的计算值

📚2 运行结果

部分代码：

function ax4 = OPplot( x, order, delay, windowSize, forAxis, timeUnits )
% @brief plots the sequence of ordinal patterns for time series x
% @author Valentina Unakafova
% @date 18.07.2017

    
nPoints = max( size( x ) ); % length of time series x
load( ['table' num2str( order ) '.mat'] );
patTable    = eval(['table' num2str(order)]);
dTau        = order*delay;
order1      = order + 1;
ordNumArray = zeros( 1, windowSize ); % ordinal numbers for all tau
nPatterns   = factorial( order1 );      % the number of ordinal patterns of order d    
ancNumArray = nPatterns./factorial( 2:order1 ); % the ancillary numbers for computation ordinal numbers
ordDist     = zeros( 1, nPatterns ); % the distribution of the ordinal patterns
previousPatternAr = zeros( 1, delay );
% for the first window
for iTau = 1:delay  
    counter = iTau; 
    % computation of the first ordinal number
    invNum      = zeros( 1, order );    
    invNum( 1 ) = ( x( dTau + iTau - delay ) >= x( dTau + iTau ) );
    for j = 2:order
        invNum( j ) = sum( x( ( order - j )*delay + iTau ) >= ...
          x( ( order1 - j )*delay + iTau:delay:dTau + iTau ) );
    end        
    ordNumArray( counter ) = sum( invNum.*ancNumArray ); % the first ordinal pattern
    ordDist( ordNumArray( counter ) + 1 ) = ordDist( ordNumArray( counter ) + 1 ) + 1;
    % the ordinal distribution for the first window
    for j = dTau+delay+iTau:delay:windowSize+dTau
        counter = counter+delay;                               
        pos = sum( x( j - dTau:delay:j ) >= x( j ) );
        ordNumArray( counter ) = patTable( ordNumArray( counter - delay )*order1 + pos );
        ordDist( ordNumArray( counter ) + 1 ) = ordDist( ordNumArray( counter ) + 1 ) + 1;            
    end  
    previousPatternAr( iTau ) = ordNumArray( counter );
 end
       
 ordPatPlot( 1:windowSize ) = ordNumArray;% for plot    
 clear ancNumArray invNum pointPos pointPos;
 
 iTau = delay;  % current shift 1:TAU
 iPattern = 1;% position of the current pattern in the window
 % a cycle over the all points in the time series
 x = x( windowSize:nPoints );
 for iPos = 2:nPoints-windowSize-dTau
   % next ordinal number computation                   
   pos = 1;
   for j = iPos:delay:iPos+dTau-delay
       pos = pos + ( x( j ) >= x( iPos + dTau ) );         
   end                  
        
   nNew = patTable( previousPatternAr( iTau )*order1 + pos );% the next ordinal pattern inside the window  
   % permutation entropy for current window computation
   nOut = ordNumArray( iPattern );          
   previousPatternAr( iTau ) = nNew;
   ordNumArray( iPattern )   = nNew; 
   
   nNew = nNew + 1;
   nOut = nOut + 1;       
   if nNew ~= nOut
            % update the ordinal distribution  
      ordDist( nNew ) = ordDist( nNew ) + 1; % new pattern
      ordDist( nOut ) = ordDist( nOut ) - 1; % old pattern
   end
        
   iTau     = iTau + 1;
   iPattern = iPattern + 1;
   if ( iTau > delay ) 
     iTau = 1; 
   end
   if ( iPattern > windowSize ) 
     iPattern = 1; 
   end
   ordPatPlot( iPos + windowSize - 1 ) = nNew; % for plots
 end
    
 % PLOTTING
 plotFontSize   = 8;
 % please, change MarkerSize to 0.1 when saving figure to EPS file
 plotMarkerSize = 0.5;
 ordPatPlot     = ordPatPlot/nPatterns;  
 ordPatPlot     = 1 - ordPatPlot;
 ax4            = subplot( 4, 1, 2 ); 
 set( gca, 'FontSize', plotFontSize ); 
    
 plotLength = length( ordPatPlot );
 plot( forAxis( 1:plotLength ), ordPatPlot( 1:plotLength ), 'ko', ...
                                          'MarkerSize', plotMarkerSize ); 
 grid on;
 axis( [ 0 forAxis( plotLength ) 0 1.0 ] );
 if ( exist( 'timeUnits' ) )
    xlabel( timeUnits, 'FontSize', plotFontSize );
 else
    xlabel( 'Time', 'FontSize', plotFontSize );
 end
 title( '(b) ORDINAL PATTERNS', 'FontSize', plotFontSize );   

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1] Amigo, J.M., Keller, K. and Unakafova, V.A., 2015. On entropy, entropy-like quantities, and applications. Discrete & Continuous Dynamical Systems-Series B, 20(10). 
[2] Bandt C., Pompe B., Permutation entropy: a natural complexity measure for time series. Physical review letters, 2002, APS 
[3] Keller, K., and M. Sinn. Ordinal analysis of time series. Physica A: Statistical Mechanics and its Applications 356.1 (2005): 114--120 
[4] Keller, K., Unakafov, A.M. and Unakafova, V.A., 2014. Ordinal patterns, entropy, and EEG. Entropy, 16(12), pp.6212-6239. 
[5] Zanin, M., Zunino, L., Rosso, O.A. and Papo, D., 2012. 
Permutation entropy and its main biomedical and econophysics applications: a review. Entropy, 14(8), pp.1553-1577. 
[6] Unakafova, V.A., Keller, K., 2013. Efficiently measuring complexity on the basis of real-world Data. Entropy, 15(10), 4392-4415. 
[7] Unakafova, V.A., 2015. Investigating measures of complexity for dynamical systems and for time series (Doctoral dissertation, University of Luebeck). 
[8] Bian, C., Qin, C., Ma, Q.D. and Shen, Q., 2012. Modified permutation-entropy analysis of heartbeat dynamics. Physical Review E, 85(2), p.021906. 
[9] Amigo, J.M., Zambrano, S. and Sanjuan, M.A., 2008. Combinatorial detection of determinism in noisy time series. EPL (Europhysics Letters), 83(6), p.60005. 
[10] Cao, Y., Tung, W.W., Gao, J.B. et al., 2004. Detecting dynamical changes in time series using the permutation entropy. Physical Review E, 70(4), p.046217. 
[11] Riedl, M., Muller, A. and Wessel, N., 2013. Practical considerations of permutation entropy. The European Physical Journal Special Topics, 222(2), pp.249-262.

🌈4 Matlab代码实现