近似熵(Approximate Entropy)

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

算法

定义一个包含 $N$ 个数据点时间序列: $u(1),u(2),\dots,u(N)$
设定一个 $m$ 表示窗口的长度，生成一组维数为 $m$ 的向量: $\mathbf{x}(1),\mathbf{x}(2),\dots,\mathbf{x}(N-m+1)$ , 其中: $\mathbf{x}(i)=\left\{u{i},u(i+1),\dots,u(i+m-1)\right\},i=1\to N-m+1$
定义 $\mathbf{x}(i)$ 和 $\mathbf{x}(j)$ 的距离 $d[\mathbf{x}(i),\mathbf{x}(j)]$ 为两者对应元素中差值最大的一个，即: $d[\mathbf{x}(i),\mathbf{x}(j)]=max_{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$ 间的距离
给定指定阈值 $r$ 对每一个 $i$ 值统计距离 $d$ 小于 $r$ 的数目及此数目与距离总数 $N - m$ 的比值，记为 $C_{i}^{m}(r)$ 即: $C_{i}^{m}(r)=\frac{1}{N-m}\left\{d[\mathbf{x}(i),\mathbf{x}(j)]<r\right\}$
现将 $C_{i}^{m}(r)$ 取对数，再求其对所有 $i$ 的平均值，记作: $\phi^{m}(r)=\frac{1}{N-m+1}\sum^{N-m+1}_{i=1}lnC_{i}^{m}(r)$
再讲维数加 $1$ ，变为 $m + 1$ ，重复步骤 $2\to5$ ，得到 $C_{i}^{m+1}(r)$ 和 $\phi^{m+1}(r)$
理论上该序列的近似熵为: $ApEn(m,r)=lim_{N\to \infty}[\phi^{m}(r)-\phi^{m+1}(r)]$

一般而言，此极限值以概率1存在，实际工作时 $N$ 不可能为 $\infty$ 。当 $N$ 为有限值时按上述步骤得到的是序列长度为 $N$ 时的 $A p E n$ 估计值。记为: $ApEn(m,r,N)=\phi^{m}(r)-\phi^{m+1}(r)$
根据实践，建议 $m=2,r=0.1\sim 0.25SD_{u}$ ， $SD_{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)