近似熵(Approximate Entropy)

• 统计学中，近似熵ApEn是一种用于量化时间序列数据波动的规律性和不可预测性¹的非线性分析技术。最初，规律性是通过精确的规律性统计来衡量的，其要集中在利用各种熵进行测量。然而，精确的熵计算需要大量的数据，并且结果会受到系统噪声的极大影响，因此将这些方法应用于实验数据是不切实际的。ApEn由Steve M. Pincus开发，通过修改精确的正则性统计量Kolmogorov-Sinai熵来处理这些限制。ApEn最初是为了分析医疗数据而开发的²，例如心率¹，后来将其应用于金融³，生理学⁴和气候科学⁵。
• 优点:
  • 只要比较短的数据就能得出比较稳健的估计值 , 所需数据点数大致是10-5000点, 一般是1000点左右；
  • 采用短数据估计特征量的可以随着时间过程的发展, 采用滑动窗来估计特征参数随时间的变化。这种动态信息往往是实际工作中更需要的。
  • 较好的抗噪和抗干扰能力。特别是对偶尔产生的瞬态强干扰有较好的承受能力；
  • 无论信号是随机的或确定性的都可以使用，因此也可以用于随机成分和确定性；

算法

1. 定义一个包含$N$个数据点时间序列: $u(1),u(2),\dots ,u(N)$
2. 设定一个$m$表示窗口的长度，生成一组维数为$m$的向量: $\mathbf{x}(1),\mathbf{x}(2),\dots ,\mathbf{x}(N-m+1)$, 其中: x ( i ) = { u i , u ( i + 1 ) , … , u ( i + m − 1 ) } , i = 1 → N − m + 1 \mathbf{x}(i)=\left\{u{i},u(i+1),\dots,u(i+m-1)\right\},i=1\to N-m+1 x(i)={ui,u(i+1),…,u(i+m−1)},i=1→N−m+1
3. 定义$\mathbf{x}(i)$和$\mathbf{x}(j)$的距离$d[\mathbf{x}(i),\mathbf{x}(j)]$为两者对应元素中差值最大的一个，即:$$d[\mathbf{x}(i),\mathbf{x}(j)]=ma{x}_{k=0\to m-1}[\left|x(i+k)-x(j+k)\right|]$$， 并对每一个$i$值计算$x(i)$与其余矢量$x(j),j=1\to N-m+1,j\ne i$间的距离
4. 给定指定阈值$r$对每一个$i$值统计距离$d$小于$r$的数目及此数目与距离总数$N-m$的比值，记为$C_i^m(r)$即: C i m ( r ) = 1 N − m { d [ x ( i ) , x ( j ) ] < r } C_{i}^{m}(r)=\frac{1}{N-m}\left\{d[\mathbf{x}(i),\mathbf{x}(j)]<r\right\} Cim​(r)=N−m1​{d[x(i),x(j)]<r}
5. 现将$C_i^m(r)$取对数，再求其对所有$i$的平均值，记作:$${\varphi }^m(r)=\frac{1}{N-m+1}\sum _{i=1}^{N-m+1}ln{C}_i^m(r)$$
6. 再讲维数加$1$，变为$m+1$，重复步骤$2\to 5$，得到$C_i^{m+1}(r)$和${\varphi }^{m+1}(r)$
7. 理论上该序列的近似熵为: $$ApEn(m,r)=li{m}_{N\to \infty }[{\varphi }^m(r)-{\varphi }^{m+1}(r)]$$

• 一般而言，此极限值以概率1存在，实际工作时$N$不可能为$\infty$。当$N$为有限值时按上述步骤得到的是序列长度为$N$时的$ApEn$估计值。记为:$$ApEn(m,r,N)={\varphi }^m(r)-{\varphi }^{m+1}(r)$$
• 根据实践，建议$m=2,r=0.1\sim 0.25S{D}_u$，$S{D}_u$为原始数据$u$的标准差

Python实现

def ApEn(time_series, m=2, r=0.15):
    """
    Approximate Entropy

    Parameters
    ----------
    time_series: {array-like}, 1D data 
    m: int, Embedding dmension
    r: float, Radius distance threshold

    Return
    ----------
    The approximate entropy estimates of the data sequence
    """
    def max_dist(x_i, x_j):
        return max([abs(ia - ja) for ia, ja in zip(x_i, x_j)])

    def phi(m):
        x = [[time_series[j] for j in range(i, i + m - 1 + 1)]
             for i in range(N - m + 1)]
        C = [
            len([1 for x_j in x if max_dist(x_i, x_j) <= r]) / (N - m + 1.0)
            for x_i in x
        ]
        return (N - m + 1.0)**(-1) * sum(np.log(C))

    N = len(time_series)
   
    return phi(m) - phi(m + 1)