演算法之python實現近似熵、互近似熵演算法

理論基礎近似熵？

• 定義：近似熵是一個隨機複雜度，反應序列相鄰的m個點所連成折線段的模式的互相近似的機率與由m+1個點所連成的折線段的模式相互近似的機率之差。

• 作用：用來描述複雜系統的不規則性，越是不規則的時間序列對應的近似熵越大。反應維數改變時產生的新的模式的可能性的大小。

對於eeg訊號來說，由於雜訊存在、和訊號的微弱性、多重訊號源疊加，反映出來的是混沌屬性，但是同一個人在大腦活動相對平穩的情況下，其eeg近似熵應該變化不大。

• 證明和對應幾何意義可參考論文：wenku.baidu.com/view/4ec89e44b307e87101f696ef.html

互近似熵

• 從近似熵定義引申出來的，近似熵描述的是一段序列的自相似程度，互近似熵比較的是兩段序列的複雜度接近程度；熵值越大越不相似，越小越相似；

近似熵演算法分析

1. 設存在一個以等時間間隔採樣獲得的m維的時間序列u(1),u(2),...,u(N).

2. 定義相關參數維數m,一般取值為2，相似容限即閥值r，其中，維數表示向量的長度；r表示“相似性”的度量值.

3. 重構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)|

4. 統計所有的d[X,X∗]<=r的個數g,則g/(N-M)就是本次的i取值對應的相似機率，計算所有i和j取值的機率對數的平均值，即熵值Φm(r)；

5. 取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代碼實現使用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倍以上。

• 參考：

• wenku.baidu.com/view/4ec89e44b307e87101f696ef.html

• http://blog.sina.com.cn/s/blog_6276ec79010118cx.html

• 79707169

• 天宇之遊
• 出處：http://www.cnblogs.com/cwp-bg/
• 本文著作權歸作者和部落格園共有，歡迎、交流，但未經作者同意必須保留此段聲明，且在文章明顯位置給出原文連結。