【Python】近似熵,样本熵,模糊熵计算高效版

前言

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

整体思路

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

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

1 近似熵(Approximate Entropy, ApEn)

1.1 理论基础

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

  它的计算方法如下:【图片来自文献1

【Python】近似熵,样本熵,模糊熵计算高效版

【Python】近似熵,样本熵,模糊熵计算高效版

  这里有一个点需要注意,那就是近似熵在计算相似向量的个数时,是会包含其自身的,也就是说,总的矢量个数为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 理论基础

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

【Python】近似熵,样本熵,模糊熵计算高效版

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

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

请添加图片描述

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上的论文4,截图如下:

【Python】近似熵,样本熵,模糊熵计算高效版

【Python】近似熵,样本熵,模糊熵计算高效版

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. ↩︎

文章出处登录后可见!

已经登录?立即刷新

共计人评分,平均

到目前为止还没有投票!成为第一位评论此文章。

(0)
心中带点小风骚的头像心中带点小风骚普通用户
上一篇 2023年3月29日
下一篇 2023年3月29日

相关推荐