Python近似熵，样本熵，模糊熵计算高效版

Posted 2023-03-28 记录无知岁月

文章目录

• 前言
• 整体思路
• 1 近似熵(Approximate Entropy, ApEn)
• • 1.1 理论基础
  • 1.2 python第三方库实现
  • 1.3 基于多线程numpy矩阵运算实现
• 2 样本熵 (Sample Entropy, SampEn)
• • 2.1 理论基础
  • 2.2 python第三方库实现
  • 2.3 基于多线程numpy矩阵运算实现
• 3 模糊熵
• • 3.1 理论基础
  • 3.2 python第三方库实现
  • 3.3 基于numpy实现
• 总结
• 参考文献

前言

  最近在学习机器学习，发现对于与生物医学信号相关的机器学习任务，在选定特征时，各种针对时间序列的熵是绕不开的重要特征，诸如近似熵，样本熵，模糊熵等。因为它们所包含的信息要远比均值方差等特征要多得多，通过写python程序实现的过程中收获了不少，这里简单总结一下。

整体思路

  写好一个程序的前提一定是理解其具体的计算思路，所以在写程序之前首先需要知道那些熵到底是怎么算出来的，这里个人强烈建议直接查找文献，而不要完全依赖各种二手的博客，因为有可能描述不清楚而直接导致程序写错。

一个讲得很全面，但程序编写不友好的教程

1 近似熵(Approximate Entropy, ApEn)

1.1 理论基础

  近似熵是Pincus在1991年提出的一种只需要较短数据就能表现信号的动力学参数，它具有以下特点：

• 只需要比较短的数据就能得到比较稳健的估计值，所需要的数据点数大致是100~5000点，一般是1000点左右；
• 有较好的抗噪和抗干扰能力。特别是对偶尔产生的瞬态强干扰有较好的承受能力。

  它的计算方法如下：【图片来自文献¹】

  这里有一个点需要注意，那就是近似熵在计算相似向量的个数时，是会包含其自身的，也就是说，总的矢量个数为N-m+1，这一点在程序编写时要尤为注意。

1.2 python第三方库实现

  像这种非常常用的熵，必然会有第三方库整合其算法，这里推荐使用的是EntropyHub这个库。里面包含了多种熵的计算方式。其使用方式如下。

import EntropyHub as EH
import numpy as np
def ApEn (Datalist, r=0.2, m=2):
	th = r * np.std(Datalist)
	return EH.ApEn(Datalist,m,r=th)[0][-1]

需要注意的是，里面的阈值容限r是指绝对量，这里是强制将其转化为相对量，即几倍的标准差。

1.3 基于多线程numpy矩阵运算实现

  打开第三方库中函数的定义，可以发现其计算熵的方式是基于循环来计算的，因此效率不是特别高。加上计算近似熵一般都有几千个点，计算起来会非常慢。而如果通过numpy矩阵运算实现，再加上多线程，其速度会快很多。

不熟悉numpy使用的可以看一下我另外一篇博客

其代码如下所示。

from pathos.multiprocessing import ThreadPool as Pool #多线程
import numpy as np
def ApEn2 (s :list|np.ndarray, r:float, m :int =2):
	s = np.squeeze(s)
	th = r * np.std(s) #容限阈值
	def phi (m):
		n = len(s)
		x = s[ np.arange(n-m+1).reshape(-1,1) + np.arange(m) ]
		ci = lambda xi: (( np.abs(x-xi).max(1) <=th).sum()) / (n-m+1) # 构建一个匿名函数
		c = Pool().map (ci, x) #所传递的参数格式: 函数名,函数参数
		return np.sum(np.log(c)) /(n-m+1)

  值得一提的是，这里用到了numpy的广播模式，即如果两个不同型的矩阵相加减，其会自动复制矩阵内的数值，使其成为同型矩阵，然后再加减。举个例子

import numpy as np

A = np.array([1,2])
B = np.array([[1,2],
              [2,3]])

C = B - A  # type: ignore
print(C)

#输出:
>>> [[0 0]
 [1 1]]

---

2 样本熵 (Sample Entropy, SampEn)

2.1 理论基础

  样本熵是基于近似熵的改进，计算方式非常类似，但是也有一些差别。其计算方式如下图所示，注意红字哦~ 【截图来自文献²】

  这里个人觉得计算方式有点奇怪，假设m初始值为2，根据上面的计算方式，当m=2时，将原时间序列划分为子序列时，最后一个点x(N)是不考虑的，这样就能得到N-2个子序列，而不是N-1个。但是当m加1之后，划分子序列时又要考虑最后一个点，因此最后得到的子序列还是N-2个。
  关于这个问题，如果要强制理解，那只能理解为要保证两种划分模式下子序列个数相同。
  这里我一开始理解错了，因为很多博客喜欢直接说“将m变成m+1，重复上述过程”，但实际上似乎不是只更换参数的意思，之所以这样理解是因为我试了好几个第三方库，它们的计算结果就是按照上面这种理解方式。

/*【2022.11.16】找到一篇文献³提到了样本熵的这个特点，并且强调这个就是样本熵相比于近似熵的一个重要改进点：*/

2.2 python第三方库实现

  这里推荐使用的第三方库还是上面提到的EntropyHub，它里面也有计算样本熵的函数。如下所示。

import EntropyHub as EH
import numpy as np
def SampleEntropy2(Datalist, r, m=2):
	th = r * np.std(Datalist) #容限阈值
	return EH.SampEn(Datalist,m,r=th)[0][-1]

2.3 基于多线程numpy矩阵运算实现

  正如上面的注意部分所强调的，在写代码的时候要尤为注意，就是子序列划分那一块的代码。

from pathos.multiprocessing import ThreadPool as Pool #多线程
import numpy as np
def SampleEntropy(Datalist, r, m=2):
	list_len = len(Datalist)  #总长度
	th = r * np.std(Datalist) #容限阈值

	def Phi(k):
		list_split = [Datalist[i:i+k] for i in range(0,list_len-k+(k-m))] #将其拆分成多个子列表
		#这里需要注意，2维和3维分解向量时的方式是不一样的！！！
		Bm = 0.0
		for i in range(0, len(list_split)): #遍历每个子向量
			Bm += ((np.abs(list_split[i] - list_split).max(1) <= th).sum()-1) / (len(list_split)-1) #注意分子和分母都要减1
		return Bm
	## 多线程
	# x = Pool().map(Phi, [m,m+1])
	# H = - math.log(x[1] / x[0]) 
	H = - math.log(Phi(m+1) / Phi(m))
	return H

除用python实现外，这里还有一个是用MATLAB写的计算样本熵的代码，也可以看看。链接

---

3 模糊熵

3.1 理论基础

  模糊熵是在样本熵的基础上进行改进得到的。从上面对样本熵的表述来看，在判断一个序列与另一个序列是否近似时，使用的是阶跃判断，即只有相似(1)和不相似(0)之间的判断，而模糊熵则是引入了一个相似度的概念，类似于模糊控制中的隶属度。

对模糊控制不熟悉的同学也可以看一下我的另外一篇博客：模糊控制算法

  关于模糊熵的计算方式，发现网上很多博客甚至很多文献（也不知道咋参考的。。。）在计算模糊熵这块有很多版本。所以为了得到正确答案，我参考了模糊熵的“鼻祖论文”——陈伟婷在2007发表在IEEE上的论文⁴，截图如下：

3.2 python第三方库实现

  如果数据量小且不想写代码的可以参考使用第三方库。

import EntropyHub as EH
import numpy as np
def FuzzyEn2(s:np.ndarray, r=0.2, m=2, n=2):
	th = r * np.std(s)
	return EH.FuzzEn(s, 2, r=(th, n))[0][-1]

3.3 基于numpy实现

  这里需要注意的就是对计算规则的理解。其代码如下所示。

import numpy as np
import math
def FuzzyEn(s, r = 0.2, m = 2, n = 2):
	'''s:需要计算熵的向量; r:阈值容限(标准差的系数); m:向量维数; n:模糊函数的指数
	'''
	N = len(s)  #总长度
	th = r * np.std(s) #容限阈值

	def Phi(k):
		list_split = [s[i:i+k] for i in range(0,N-k+(k-m))] #将其拆分成多个子列表
		#这里需要注意，2维和3维分解向量时的方式是不一样的！！！
		B = np.zeros(len(list_split))
		for i in range(0, len(list_split)): #遍历每个子向量
			di = np.abs(list_split[i] - np.mean(list_split[i]) - list_split + np.mean(list_split,1).reshape(-1,1)).max(1)
			Di = np.exp(- np.power(di,n) / th)
			B[i] = (np.sum(Di) - 1) / (len(list_split)-1) #这里减1是因为要除去其本身，即exp(0)
		return np.sum(B) / len(list_split)
	H = - math.log(Phi(m+1) / Phi(m))
	return H

总结

  计算各种熵的关键还是在于对计算方式的理解，如果博客说法不一，那就去查找文献，如果文献说法不一，那就去找提出这个熵的论文。
  计算速度方面，发现对于较大的数据量，如3000，第三方库计算近似熵和样本熵的速度比numpy矩阵运算速度慢，但模糊熵计算速度却比numpy矩阵运算速度快很多。

但是按理说，模糊熵的复杂度至少是样本熵的两倍，但是模糊熵的计算速度却比样本熵快，我估计是第三方库作者可能是觉得样本熵的代码太简单，没有必要进行优化。

参考文献

---

1. 杨福生,廖旺才.近似熵:一种适用于短数据的复杂性度量[J].中国医疗器械杂志,1997(05):283-286. ↩︎

2. 赵志宏, 杨绍普.一种基于样本熵的轴承故障诊断方法[J].2012-06. ↩︎

3. 刘慧, 谢洪波, 和卫星, 等. 基于模糊熵的脑电睡眠分期特征提取与分类[J]. 数据采集与处理,2010,25(4):484-489. ↩︎

4. Chen W, Wang Z, Xie H, Yu W. Characterization of surface EMG signal based on fuzzy entropy. IEEE Trans Neural Syst Rehabil Eng. 2007 Jun;15(2):266-72. doi: 10.1109/TNSRE.2007.897025. PMID: 17601197. ↩︎

排列熵模糊熵近似熵样本熵的原理及MATLAB实现之近似熵

说明：“本博文为排列熵、模糊熵、近似熵、样本熵的原理及MATLAB实现”系列博文的最后一篇，关于排列熵、模糊熵、样本熵的内容请阅读博客：

排列熵

模糊熵

样本熵

近似熵

• • 四、近似熵
  • • 1.简介
    • 2.基本原理
    • 3.MATLAB代码
    • 参考文献

四、近似熵

1.简介

近似熵（approximate entropy，ApEn）可以定量描述时间序列的复杂程度，序列的复杂性越大，相应的近似熵也越大。近似熵的值受数据量的影响较小，对于非平稳、非线性序列的量化结果稳定，在实际工程中得以广泛应用。

2.基本原理

设有长度为 $N$ 的时间序列$X=[x_1,x_2,\dots ,x_N]$，其近似熵的计算步骤如下：
Step1：将时间序列 $X$ 的元素按顺序排列为具有 $m$ 维数的向量，即
$$X_i=[x(i),x(i+1),...,x(i+m-1)]$$
式中，$i=1,2,\dots ,N-m+1$。
Step2：定义 $d[X_i,X_j]$ 为向量 $X_i$ 与 $X_j$ 的距离，则：
$$d[X_i,X_j]=max|x(i+k)-x(j+k)|，k\in (0,m-1)$$
Step3：记 $B_i$ 为 $d[X_i,X_j]\le r$的个数( $r$ 为相似容限)，并计算 $B_i$ 与全部矢量数 $N-m+1$ 的比值，即：
$$B_i^m(r)=\frac{B_i}{N-m+1}$$
Step4：对 $B_i^m(r)$ 进行取对数运算，再求其对所有 $i$ 的平均值，记作 $B^m(r)$ ，则有：
$$B^m(r)=\frac{1}{N-m+1}\sum _{i=1}^{N-m+1}ln{B}_i^m(r)$$
Step5：令 $m=m+1$，并重复 Step1~Step4，即可得到 $B^{m+1}(r)$ 。
Step6：理论上，此序列的近似熵为：
$$ApEn(m,r)=lim[B^m(r)-B^{m+1}(r)]，N\to \infty$$
对于实际的序列，$N$ 不可能趋近无穷大，因此近似熵可表示为：
$$ApEn(m,r,N)=B^m(r)-B^{m+1}(r)$$
近似熵本质上是一个关于序列和参数的统计值，它的大小与数据长度 $N$ 、嵌入维数 $m$ 和相似容限 $r$ 有关。为了得到较好的统计特性以及较小的误差， 数据长度$N$ 通常在100～5000 取值，嵌入维数 $m$ 一般取 1 或 2，相似容限 $r$ 取(0.1～0.25)*$std$，$std$ 为序列的标准差。

3.MATLAB代码

% 主程序
clc;
clear;
close all;

%% 产生仿真信号
fs = 1000;	% 数据采样率
t = (0:1/fs:(1-1/fs));          % 时间
x = cos(50*pi*t+sin(5*pi*t));   % 数据

%% 画图
figure;
plot(t,x);
xlabel('t/s');ylabel('幅值');title('信号的时域波形');

%% 求仿真信号x的近似熵
m = 2;                                  % 嵌入维数
r0 = 0.2;                               % 相似容限的系数
r = r0*std(x);                          % 相似容限
appEn = ApproximateEntropy(m,r,x);   % 近似熵

% 求近似熵的函数
function appEn = ApproximateEntropy(dim, r, data, tau)
%   近似熵算法的提出者：Pincus S M . Approximate entropy as a measure of system complexity[J]. Proceedings of the National Academy of Sciences ,1991,88(6):2297—2301.
%   Input:
%       dim：嵌入维数(一般取1或者2)
%       r：相似容限( 通常取0.1*Std(data)~0.25*Std(data) )
%       data：时间序列数据，data须为1xN的矩阵
%       tau：下采样延迟时间（在默认值为1的情况下，用户可以忽略此项）
%   Output:
%       appEn：所求数据的近似熵
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(N-m+1,1);
    dataMat = zeros(N-m+1,m);
    
    % 设置数据矩阵，构造成m维的矢量
    for i = 1:N-m+1
        dataMat(i,:) = data(1,i:i+m-1)排列熵模糊熵近似熵样本熵的原理及MATLAB实现之近似熵

 排列熵模糊熵近似熵样本熵的原理及MATLAB实现
 样本熵（SampEn）计算时间序列复杂度（Python程序）
 matlab计算样本熵
 信息熵 增益
 通俗的解释交叉熵与相对熵