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

Posted 记录无知岁月

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了Python近似熵,样本熵,模糊熵计算高效版相关的知识,希望对你有一定的参考价值。

文章目录

前言

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

整体思路

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

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

1 近似熵(Approximate Entropy, ApEn)

1.1 理论基础

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

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

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

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

  这里个人觉得计算方式有点奇怪,假设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,截图如下:

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.简介

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

2.基本原理

设有长度为 N N N 的时间序列 X = [ x 1 , x 2 , … , x N ] X= [x_1, x_2,…,x_N ] X=[x1,x2,,xN],其近似熵的计算步骤如下:
Step1:将时间序列 X X X 的元素按顺序排列为具有 m m m 维数的向量,即
X i = [ x ( i ) , x ( i + 1 ) , . . . , x ( i + m − 1 ) ] X_i=[x(i),x(i+1),...,x(i+m-1)] Xi=[x(i),x(i+1),...,x(i+m1)]
式中, i = 1 , 2 , … , N − m + 1 i=1, 2 , … , N-m+1 i=1,2,,Nm+1
Step2:定义 d [ X i , X j ] d [X_i, X_j] d[Xi,Xj] 为向量 X i X_i Xi X j X_j Xj 的距离,则:
d [ X i , X j ] = m a x ∣ x ( i + k ) − x ( j + k ) ∣ , k ∈ ( 0 , m − 1 ) d [X_i, X_j]=max|x(i+k)-x(j+k)|,k∈(0,m-1) d[Xi,Xj]=maxx(i+k)x(j+k)k(0,m1)
Step3:记 B i B_i Bi d [ X i , X j ] ≤ r d [X_i, X_j] ≤ r d[Xi,Xj]r的个数( r r r 为相似容限),并计算 B i B_i Bi 与全部矢量数 N − m + 1 N-m+1 Nm+1 的比值,即:
B i m ( r ) = B i N − m + 1 B^m_i(r)=\\fracB_iN-m+1 Bim(r)=Nm+1Bi
Step4:对 B i m ( r ) B^m_i(r) Bim(r) 进行取对数运算,再求其对所有 i i i 的平均值,记作 B m ( r ) B^m(r) Bm(r) ,则有:
B m ( r ) = 1 N − m + 1 ∑ i = 1 N − m + 1 l n B i m ( r ) B^m(r)=\\frac1N-m+1\\sum_i=1^N-m+1lnB^m_i(r) Bm(r)=Nm+11i=1Nm+1lnBim(r)
Step5:令 m = m + 1 m=m+1 m=m+1,并重复 Step1~Step4,即可得到 B m + 1 ( r ) B^m+1(r) Bm+1(r)
Step6:理论上,此序列的近似熵为:
A p E n ( m , r ) = l i m [ B m ( r ) − B m + 1 ( r ) ] , N → ∞ ApEn(m,r)=lim[B^m(r)-B^m+1(r)],N→∞ ApEn(m,r)=lim[Bm(r)Bm+1(r)]N
对于实际的序列, N N N 不可能趋近无穷大,因此近似熵可表示为:
A p E n ( m , r , N ) = B m ( r ) − B m + 1 ( r ) ApEn(m,r,N)=B^m(r)-B^m+1(r) ApEn(m,r,N)=Bm(r)Bm+1(r)
近似熵本质上是一个关于序列和参数的统计值,它的大小与数据长度 N N N 、嵌入维数 m m m 和相似容限 r r r 有关。为了得到较好的统计特性以及较小的误差, 数据长度 N N N 通常在100~5000 取值,嵌入维数 m m m 一般取 1 或 2,相似容限 r r r 取(0.1~0.25)* s t d std std s t d 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):22972301.
%   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计算样本熵

信息熵 增益

通俗的解释交叉熵与相对熵