理論基礎
近似熵？
-
定義：近似熵是一個隨機復雜度，反應序列相鄰的m個點所連成折線段的模式的互相近似的概率與由m+1個點所連成的折線段的模式相互近似的概率之差。
-
作用：用來描述復雜系統的不規則性，越是不規則的時間序列對應的近似熵越大。反應維數改變時產生的新的模式的可能性的大小。
對於eeg信號來說，由於噪聲存在、和信號的微弱性、多重信號源疊加，反映出來的是混沌屬性，但是同一個人在大腦活動相對平穩的情況下，其eeg近似熵應該變化不大。
- 證明和對應幾何意義可參考論文：https://wenku.baidu.com/view/4ec89e44b307e87101f696ef.html
互近似熵
- 從近似熵定義引申出來的，近似熵描述的是一段序列的自相似程度，互近似熵比較的是兩段序列的復雜度接近程度；熵值越大越不相似，越小越相似；
近似熵算法分析
-
設存在一個以等時間間隔采樣獲得的m維的時間序列u(1),u(2),...,u(N).
-
定義相關參數維數m,一般取值為2，相似容限即閥值r，其中，維數表示向量的長度；r表示“相似度”的度量值.
-
重構m維向量X(1),X(2),...,X(N−m+1),其中X(i)=[u(i),u(i+1),...,u(i+m−1)]，X(j)=[u(j),u(j+1),...,u(j+m−1)];計算X(i)和X(j)之間的距離，由對應元素的最大差值決定；d[X,X∗]=maxa|u(a)−u∗(a)|d[X,X∗]=maxa|u(a)−u∗(a)|
-
統計所有的d[X,X∗]<=r的個數g,則g/(N-M)就是本次的i取值對應的相似概率，計算所有i和j取值的概率對數的平均值，即熵值Φm(r)；
-
取m+1重復3、4過程，計算近似熵：
ApEn=Φm(r)−Φm+1(r)
參數選擇：通常選擇參數m=2或m=3；通常選擇r=0.2∗std，其中std表示原時間序列的標准差.
- 互近似熵計算和近似熵的步驟一樣，把計算X(i)和X(j)之間的距離改為計算序列a的向量X(i)和序列b的向量Y(j)的距離；相似容限r為兩個原序列的0.2倍協方差；
python代碼實現
github源碼
使用Pincus提出的近似熵定義計算近似熵
class BaseApEn(object):
"""
近似熵基礎類
"""
def __init__(self, m, r):
"""
初始化
:param U:一個矩陣列表，for example:
U = np.array([85, 80, 89] * 17)
:param m: 子集的大小，int
:param r: 閥值基數，0.1---0.2
"""
self.m = m
self.r = r
@staticmethod
def _maxdist(x_i, x_j):
"""計算矢量之間的距離"""
return np.max([np.abs(np.array(x_i) - np.array(x_j))])
@staticmethod
def _biaozhuncha(U):
"""
計算標准差的函數
:param U:
:return:
"""
if not isinstance(U, np.ndarray):
U = np.array(U)
return np.std(U, ddof=1)
class ApEn(BaseApEn):
"""
Pincus提出的算法，計算近似熵的類
"""
def _biaozhunhua(self, U):
"""
將數據標准化，
獲取平均值
所有值減去平均值除以標准差
"""
self.me = np.mean(U)
self.biao = self._biaozhuncha(U)
return np.array([(x - self.me) / self.biao for x in U])
def _dazhi(self, U):
"""
獲取閥值
:param U:
:return:
"""
if not hasattr(self, "f"):
self.f = self._biaozhuncha(U) * self.r
return self.f
def _phi(self, m, U):
"""
計算熵值
:param U:
:param m:
:return:
"""
# 獲取矢量列表
x = [U[i:i + m] for i in range(len(U) - m + 1)]
# 獲取所有的比值列表
C = [len([1 for x_j in x if self._maxdist(x_i, x_j) <= self._dazhi(U)]) / (len(U) - m + 1.0) for x_i in x]
# 計算熵
return np.sum(np.log(list(filter(lambda a: a, C)))) / (len(U) - m + 1.0)
def _phi_b(self, m, U):
"""
標准化數據計算熵值
:param m:
:param U:
:return:
"""
# 獲取矢量列表
x = [U[i:i + m] for i in range(len(U) - m + 1)]
# 獲取所有的比值列表
C = [len([1 for x_j in x if self._maxdist(x_i, x_j) <= self.r]) / (len(U) - m + 1.0) for x_i in x]
# 計算熵
return np.sum(np.log(list(filter(lambda x: x, C)))) / (len(U) - m + 1.0)
def jinshishang(self, U):
"""
計算近似熵
:return:
"""
return np.abs(self._phi(self.m + 1, U) - self._phi(self.m, U))
def jinshishangbiao(self, U):
"""
將原始數據標准化后的近似熵
:param U:
:return:
"""
eeg = self._biaozhunhua(U)
return np.abs(self._phi_b(self.m + 1, eeg) - self._phi_b(self.m, eeg))
if __name__ == "__main__":
U = np.array([2, 4, 6, 8, 10] * 17)
G = np.array([3, 4, 5, 6, 7] * 17)
ap = ApEn(2, 0.2)
ap.jinshishang(U) # 計算近似熵
說明:
-
jinshishang函數直接計算近似熵
-
jinshishangbiao函數將原始數據標准化后計算近似熵
使用Pincus提出的近似熵定義計算互近似熵
class HuApEn(BaseApEn):
def _xiefangcha(self, U, G):
"""
計算協方差的函數
:param U: 序列1，矩陣
:param G: 序列2，矩陣
:return: 協方差，float
"""
if not isinstance(U, np.ndarray):
U = np.array(U)
if not isinstance(G, np.ndarray):
G = np.array(G)
if len(U) != len(G):
raise AttributeError('參數錯誤！')
return np.cov(U, G, ddof=1)[0, 1]
def _biaozhunhua(self, U, G):
"""
對數據進行標准化
"""
self.me_u = np.mean(U)
self.me_g = np.mean(G)
self.biao_u = self._biaozhuncha(U)
self.biao_g = self._biaozhuncha(G)
# self.biao_u = self._xiefangcha(U, G)
# self.biao_g = self._xiefangcha(U, G)
return np.array([(x - self.me_u) / self.biao_u for x in U]), np.array(
[(x - self.me_g) / self.biao_g for x in U])
def _dazhi(self, U, G):
"""
獲取閥值
:param r:
:return:
"""
if not hasattr(self, "f"):
self.f = self._xiefangcha(U, G) * self.r
return self.f
def _phi(self, m, U, G):
"""
計算熵值
:param m:
:return:
"""
# 獲取X矢量列表
x = [U[i:i + m] for i in range(len(U) - m + 1)]
# 獲取y矢量列表
y = [G[g:g + m] for g in range(len(G) - m + 1)]
# 獲取所有的條件概率列表
C = [len([1 for y_k in y if self._maxdist(x_i, y_k) <= self._dazhi(U, G)]) / (len(U) - m + 1.0) for x_i in x]
# 計算熵
return np.sum(np.log(list(filter(lambda x_1: x_1, C)))) / (len(U) - m + 1.0)
def _phi_b(self, m, U, G):
"""
標准化數據計算熵值
:param m:
:param U:
:return:
"""
# 獲取X矢量列表
x = [U[i:i + m] for i in range(len(U) - m + 1)]
# 獲取y矢量列表
y = [G[g:g + m] for g in range(len(G) - m + 1)]
# 獲取所有的條件概率列表
C = [len([1 for y_k in y if self._maxdist(x_i, y_k) <= self.r]) / (len(U) - m + 1.0) for x_i in x]
# 計算熵
return np.sum(np.log(list(filter(lambda x: x, C)))) / (len(U) - m + 1.0)
def hujinshishang(self, U, G):
"""
計算互近似熵
:return:
"""
return np.abs(self._phi(self.m + 1, U, G) - self._phi(self.m, U, G))
def hujinshishangbiao(self, U, G):
"""
將原始數據標准化后的互近似熵
:param U:
:param G:
:return:
"""
u, g = self._biaozhunhua(U, G)
return np.abs(self._phi_b(self.m + 1, u, g) - self._phi_b(self.m, u, g))
使用洪波提出的快速實用算法計算近似熵
class NewBaseApen(object):
"""新算法基類"""
@staticmethod
def _get_array_zeros(x):
"""
創建N*N的0矩陣
:param U:
:return:
"""
N = np.size(x, 0)
return np.zeros((N, N), dtype=int)
@staticmethod
def _get_c(z, m):
"""
計算熵值的算法
:param z:
:param m:
:return:
"""
N = len(z[0])
# 概率矩陣C計算
c = np.zeros((1, N - m + 1))
if m == 2:
for j in range(N - m + 1):
for i in range(N - m + 1):
c[0, j] += z[j, i] & z[j + 1, i + 1]
if m == 3:
for j in range(N - m + 1):
for i in range(N - m + 1):
c[0, j] += z[j, i] & z[j + 1, i + 1] & z[j + 2, i + 2]
if m != 2 and m != 3:
raise AttributeError('m的取值不正確！')
data = list(filter(lambda x:x, c[0]/(N - m + 1.0)))
if not all(data):
return 0
return np.sum(np.log(data)) / (N - m + 1.0)
class NewApEn(ApEn, NewBaseApen):
"""
洪波等人提出的快速實用算法計算近似熵
"""
def _get_distance_array(self, U):
"""
獲取距離矩陣
:param U:
:return:
"""
z = self._get_array_zeros(U)
fa = self._dazhi(U)
for i in range(len(z[0])):
z[i, :] = (np.abs(U - U[i]) <= fa) + 0
return z
def _get_shang(self, m, U):
"""
計算熵值
:param U:
:return:
"""
# 獲取距離矩陣
Z = self._get_distance_array(U)
return self._get_c(Z, m)
def hongbo_jinshishang(self, U):
"""
計算近似熵
:param U:
:return:
"""
return np.abs(self._get_shang(self.m + 1, U) - self._get_shang(self.m, U))
使用洪波提出的快速實用算法計算互近似熵
class NewHuApEn(HuApEn, NewBaseApen):
"""
洪波等人提出的快速實用算法計算互近似熵
"""
def _get_distance_array(self, U, G):
"""
獲取距離矩陣
:param U:模板數據
:return:比較數據
"""
z = self._get_array_zeros(U)
fa = self._dazhi(U, G)
for i in range(len(z[0])):
z[i, :] = (np.abs(G - U[i]) <= fa) + 0
return z
def _get_shang(self, m, U, G):
"""
計算熵值
:param U:
:return:
"""
# 獲取距離矩陣
Z = self._get_distance_array(U, G)
return self._get_c(Z, m)
def hongbo_hujinshishang(self, U, G):
"""
對外的計算互近似熵的接口
:param U:
:param G:
:return:
"""
return np.abs(self._get_shang(self.m + 1, U, G) - self._get_shang(self.m, U, G))
簡單測試
if __name__ == "__main__":
import time
import random
U = np.array([random.randint(0, 100) for i in range(1000)])
G = np.array([random.randint(0, 100) for i in range(1000)])
ap = NewApEn(2, 0.2)
ap1 = NewHuApEn(2, 0.2)
t = time.time()
print(ap.jinshishang(U))
t1 = time.time()
print(ap.hongbo_jinshishang(U))
t2 = time.time()
print(ap1.hujinshishang(U, G))
t3 = time.time()
print(ap1.hongbo_hujinshishang(U, G))
t4 = time.time()
print(t1-t)
print(t2-t1)
print(t3-t2)
print(t4-t3)
-
測試后發現使用快速算法比使用定義算法的計算效率提高了6倍以上。
-
參考：
- 作者：天宇之游
- 出處：http://www.cnblogs.com/cwp-bg/
- 本文版權歸作者和博客園共有，歡迎轉載、交流，但未經作者同意必須保留此段聲明，且在文章明顯位置給出原文鏈接。