文章目录

前言
整体思路
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，这一点在程序编写时要尤为注意。

python_26">1.2 python第三方库实现

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

python">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是指绝对量，这里是强制将其转化为相对量，即几倍的标准差。

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

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

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

其代码如下所示。

python">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的广播模式，即如果两个不同型的矩阵相加减，其会自动复制矩阵内的数值，使其成为同型矩阵，然后再加减。举个例子

python">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，重复上述过程”，但实际上似乎不是只更换参数的意思，之所以这样理解是因为我试了好几个第三方库，它们的计算结果就是按照上面这种理解方式。

python_92">2.2 python第三方库实现

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

python">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]

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

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

python">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上的论文³，截图如下：

在这里插入图片描述

python_145">3.2 python第三方库实现

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

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]

numpy_156">3.3 基于numpy实现

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

python">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矩阵运算速度快很多。

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

参考文献

杨福生,廖旺才.近似熵:一种适用于短数据的复杂性度量[J].中国医疗器械杂志,1997(05):283-286. ↩︎
赵志宏, 杨绍普.一种基于样本熵的轴承故障诊断方法[J].2012-06. ↩︎
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. ↩︎