最近在学习机器学习，发现对于与生物医学信号相关的机器学习任务，在选定特征时，各种针对时间序列的熵是绕不开的重要特征，诸如近似熵，样本熵，模糊熵等。因为它们所包含的信息要远比均值方差等特征要多得多，通过写python程序实现的过程中收获了不少，这里简单总结一下。
写好一个程序的前提一定是理解其具体的计算思路，所以在写程序之前首先需要知道那些熵到底是怎么算出来的，这里个人强烈建议直接查找文献，而不要完全依赖各种二手的博客，因为有可能描述不清楚而直接导致程序写错。
近似熵是Pincus在1991年提出的一种只需要较短数据就能表现信号的动力学参数，它具有以下特点：
它的计算方法如下：【图片来自文献1】
这里有一个点需要注意，那就是近似熵在计算相似向量的个数时，是会包含其自身的，也就是说，总的矢量个数为N-m+1，这一点在程序编写时要尤为注意。
像这种非常常用的熵，必然会有第三方库整合其算法，这里推荐使用的是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是指绝对量，这里是强制将其转化为相对量，即几倍的标准差。
打开第三方库中函数的定义，可以发现其计算熵的方式是基于循环来计算的，因此效率不是特别高。加上计算近似熵一般都有几千个点，计算起来会非常慢。而如果通过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】
这里个人觉得计算方式有点奇怪，假设
m初始值为2，根据上面的计算方式，当m=2时，将原时间序列划分为子序列时，最后一个点x(N)是不考虑的，这样就能得到N-2个子序列，而不是N-1个。但是当m加1之后，划分子序列时又要考虑最后一个点，因此最后得到的子序列还是N-2个。
关于这个问题，如果要强制理解，那只能理解为要保证两种划分模式下子序列个数相同。
这里我一开始理解错了，因为很多博客喜欢直接说“将m变成m+1，重复上述过程”，但实际上似乎不是只更换参数的意思，之所以这样理解是因为我试了好几个第三方库，它们的计算结果就是按照上面这种理解方式。
/*【2022.11.16】找到一篇文献3提到了样本熵的这个特点，并且强调这个就是样本熵相比于近似熵的一个重要改进点：*/
这里推荐使用的第三方库还是上面提到的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]
正如上面的注意部分所强调的，在写代码的时候要尤为注意，就是子序列划分那一块的代码。
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写的计算样本熵的代码，也可以看看。链接
模糊熵是在样本熵的基础上进行改进得到的。从上面对样本熵的表述来看，在判断一个序列与另一个序列是否近似时，使用的是阶跃判断，即只有相似(1)和不相似(0)之间的判断，而模糊熵则是引入了一个相似度的概念，类似于模糊控制中的隶属度。
对模糊控制不熟悉的同学也可以看一下我的另外一篇博客：模糊控制算法
关于模糊熵的计算方式，发现网上很多博客甚至很多文献（也不知道咋参考的。。。）在计算模糊熵这块有很多版本。所以为了得到正确答案，我参考了模糊熵的“鼻祖论文”——陈伟婷在2007发表在IEEE上的论文4，截图如下：
如果数据量小且不想写代码的可以参考使用第三方库。
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]
这里需要注意的就是对计算规则的理解。其代码如下所示。
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. ↩︎
刘慧, 谢洪波, 和卫星, 等. 基于模糊熵的脑电睡眠分期特征提取与分类[J]. 数据采集与处理,2010,25(4):484-489. ↩︎
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. ↩︎