样本熵理论相关知识与代码实现

关于本博客的说明： 本次博客主要分享样本熵（Sample Entropy, SampEn, SE）的理论相关知识及其代码实现.
之前分享过一篇有关近似熵的博客**“近似熵理论相关知识与代码实现”**，有兴趣的博友欢迎前往观看，不足之处还请不吝赐教！https://blog.csdn.net/cratial/article/details/79707169

一、理论基础

**样本熵(SampEn)**是基于近似熵(ApEn)的一种用于度量时间序列复杂性的改进方法，在评估生理时间序列的复杂性和诊断病理状态等方面均有应用[1].
由于样本熵是近似熵的一种改进方法，因此可以将其与近似熵联系起来理解.

算法表述如下：

1. 设存在一个以等时间间隔采样获得的$N$维的时间序列$u(1),u(2),...,u(N)$.

2. 定义算法相关参数$m$,$r$，其中，$m$为整数，表示比较向量的长度，$r$为实数，表示“相似度”的度量值.

3. 重构$m$维向量$X(1),X(2),...,X(N-m+1)$,其中$X(i)=[u(i),u(i+1),...,u(i+m-1)]$.

4. 对于$1\le i\le N-m+1$，统计满足以下条件的向量个数$$B_i^m(r)=(numberofX(j)suchthatd[X(i),X(j)]\le r)/(N-m),i\ne j$$
   其中，$d[X,X^{\ast }]$定义为$d[X,X^{\ast }]=\underset{a}{\max }|u(a)-u^{\ast }(a)|,X\ne X^{\ast }$
   $u(a)$为向量$X$的元素，$d$表示向量$X(i)$与$X(j)$的距离，由对应元素的最大差值决定，$j$的取值范围为$[1,N-m+1]$，但是$j\ne i$.

5. 求$B_i^m(r)$对所有$i$值的平均值，记为$B^m(r)$，即$$B^m(r)=(N-m+1{)}^{-1}\sum _{i=1}^{N-m+1}{B}_i^m(r)$$.

6. 令$k=m+1$，重复步骤3-4，得$A^k(r)=(N-k+1{)}^{-1}{\sum }_{i=1}^{N-k+1}{A}_i^k(r)$，其中$$A_i^k(r)=(numberofX(j)suchthatd[X(i),X(j)]\le r)/(N-k),i\ne j$$

7. 则样本熵(SampEn)定义为 S a m p E n = lim ⁡ N → ∞ { − l n [ A k ( r ) / B m ( r ) ] } SampEn=\mathop {\lim }\limits_{N \to \infty } \{ -ln[A^k(r)/B^m(r)]\} SampEn=N→∞lim​{−ln[Ak(r)/Bm(r)]}
   由于在实际计算应用过程中，$N$不可能为$\infty$，因此当$N$取有限值时，样本熵估计为$$SampEn=-ln[A^k(r)/B^m(r)]\}$$

其中， $ln$表示自然对数，$m$和$r$由第2步定义.
参数选择：嵌入维数$m$一般取1或2；相似容限$r$的选择在很大程度上取决于实际应用场景，通常选择$r=0.1\ast std\sim 0.25\ast std$，其中$std$表示原时间序列的标准差.

近似熵与样本熵理论上的对比[2]
设$B$为维数为$m$时，时间序列的自相似概率；$A$为维数为$k=m+1$时，时间序列的自相似概率，得出$CP=A/B$。近似熵的计算是以$-ln(CP)$为模型，并且计算出了所有模型的平均值。为了防止出现计算$ln(0)$的情况，近似熵在算法的第4步中包含了与自身向量的比较，这种方式与新信息观点是不相容的，所以一定会存在偏差。不同的是，样本熵计算的是和的对数，且不包含与自身向量的比较，所以其优势在于包含更大的$A、B$，以及更加准确的$CP$估计.
与近似熵相比，样本熵具有两个优势：样本熵的计算不依赖数据长度；样本熵具有更好的一致性，即参数$m$和$r$的变化对样本熵的影响程度是相同的.

二、代码实现

Python
（备注：示例程序是根据上文算法在Wikipedia关于Sample entropy的python实现代码上进行修改后得到的[1]）

# SampEn calculation

import numpy as np

def SampEn(U, m, r):

     def _maxdist(x_i, x_j):
           return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

     def _phi(m):
          x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
          B = [(len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) - 1.0) / (N - m) for x_i in x]
          return (N - m + 1.0)**(-1) * sum(B)

     N = len(U)

     return -np.log(_phi(m+1) / _phi(m))

# Usage example
U = np.array([85, 80, 89] *17)
print(SampEn(U,2,3))

令输入数据为一个以3为周期的周期时间序列，维度$N=51$，即
S N = { 85 , 80 , 89 , 85 , 80 , 89 , . . . } S_N=\left \{85,80,89,85,80,89,... \right \} SN​={85,80,89,85,80,89,...}
设定参数：嵌入维度$m=2$，相似容限$r=3$
程序输出样本熵 $SampEn=0.0008507018803128114$

MATLAB
（备注：示例程序是根据上文算法在Kijoon Lee分享的样本熵代码“SampE”[3]的基础上修改得来的，在此表示感谢！）

% notification: this is a modification work based on that of Kijoon Lee.
function saen = SampEntropy( dim, r, data, tau )
% SAMPEN Sample Entropy
%   calculates the sample entropy of a given time series data

%   SampEn is conceptually similar to approximate entropy (ApEn), but has
%   following differences:
%       1) SampEn does not count self-matching. The possible trouble of
%       having log(0) is avoided by taking logarithm at the latest step.
%       2) SampEn does not depend on the datasize as much as ApEn does. The
%       comparison is shown in the graph that is uploaded.

%   dim     : embedded dimension
%   r       : tolerance (typically 0.2 * std)
%   data    : time-series data
%   tau     : delay time for downsampling (user can omit this, in which case
%             the default value is 1)
%

if nargin < 4, tau = 1; end
if tau > 1, data = downsample(data, tau); end

N = length(data);
result = zeros(1,2);

for m = dim:dim+1
    Bi = zeros(1,N-m+1);
    dataMat = zeros(m,N-m+1);
    
    % setting up data matrix
    for i = 1:m
        dataMat(i,:) = data(i:N-m+i);
    end
    
    % counting similar patterns using distance calculation
    for j = 1:N-m+1
        % calculate Chebyshev distance, excluding self-matching case
        dist = max(abs(dataMat - repmat(dataMat(:,j),1,N-m+1)));
        % calculate Heaviside function of the distance
        % User can change it to any other function
        % for modified sample entropy (mSampEn) calculation
        D = (dist <= r);
        % excluding self-matching case
        Bi(j) = (sum(D)-1)/(N-m);
    end
    
    % summing over the counts
    result(m-dim+1) = sum(Bi)/(N-m+1);
    
end

saen = -log(result(2)/result(1));

end

采用与Python实现代码相同的计算数据与参数，得到样本熵为$SampEn=8.507018803128114e-04$，此结果与Python计算结果一致.

---

参考文献：
[1]. https://en.wikipedia.org/wiki/Sample_entropy
[2]. 李立,曹锐,相洁.脑电数据近似熵与样本熵特征对比研究[J].计算机工程与设计,2014,35(03):1021-1026.
[3]. https://cn.mathworks.com/matlabcentral/fileexchange/35784-sample-entropy?s_tid=srchtitle