Recognition of Epileptic EEG Signals Using a Novel Multi-View TSK Fuzzy System论文阅读

注释

大多数的超参数的选择是通过使用交叉验证获得，将结果最优的超参数作为最终参数。

Abstrct

利用基于机器学习的方法对eeg信号进行处理从而识别癫痫疾病一般包含两个步骤。
       1：特征提取，获得具有代表性的特征。2：特征分类，建立有效的模型根据获得的特征进行分类。
存在问题：1、合适的特征提取方法难以获得，2、大多数分类方法为黑盒方法，需要更多的可解释性。
       解决方案：首先通过不同的特征提取方法生成不同视角的eeg信号数据，然后使用TSK模糊系统进行分类。基于以上统称为MV-TSK-FS。另外可通过特征提取方法识别每一个特征的权重并输出。

Introsuction

概述癫痫检测的机器学习方法以及所解决的问题。本文提出的多视角特征提取方法可以避免相似的特征提取，同时采用可视化权重的方法选择最合适的特征提取方法。
主要贡献：1、癫痫检测增加了全面性和合理性。
                  2、提出了多视角模糊系统方法，对于特征的权重进行展示，增加了可解释性
                  3、方法有更强的可解释性和生成能力，优于现有技术。

Background

Data and Feature extract methods

数据：五个组，每个组包含100个单通道eeg，即100个持续时间为23.6s的样本，采样率为173.6Hz。各个数据描述如下
特征提取解决问题：不能使用原始信号，因为他维度高、噪音较多，不平稳，非线性。
特征提取方法：主成分分析（时域特征）、傅里叶变换（频域特征）、小波分析（时频特征）。
本文主要使用：1、小波包分解；2、短时傅里叶变换；3、核主成分分析。通过以上三种方法对数据进行特征提取，生成三个视角的结果。
例如：健康组A组的数据进行三种方法后的结果如图        如图a所示，通过小波包分解为6个波段，最后一个A5是近似值，b和c同理，分为多个波段。       图b是短时傅里叶变换，采用一个小的滑动窗口进行特征提取，谱值计算是通过1秒计算每0.5s的数据。首先通过STFT转化为局部平稳段，再使用傅里叶变换获得谱值，最后划分其为六个频率带宽（根据频率）
       图c是核主成分分析图，这个方法实现了复杂的非线性映射，采用高斯核函数得到六个主成分以及特征值。

TSK fuzzy system

       本文以TSK为基准模型进行扩充权重机制以及多视角机制。TSK介绍如下：
模糊推理规则一般为与规则
这里采用了if…then规则，可以参考书籍或论文。并对其进行了解释。这里k是模糊规则数量，一般设置为5，10，15视具体情况而定，也可以自学习。AK表示第k个模糊规则的子集（输入向量的一部分），将输入向量x映射到模糊集合中在进行模糊规则操作得到fk(x)在以前的TSK模型中，其输出y可以表示为隶属度乘以函数之和。在隶属度函数上可以选择高斯径向基函数，其c的选择可以通过FCM算法进行获得，具体步骤如下：将以上公式以及结果进行简化
总结得到的公式如上。将各个输入向量以及计算过程用矩阵表示，再利用矩阵特性从而得到结果。
以上便是TSK模糊系统的过程，在传统的TSK模糊系统学习过程，是通过一个目标函数且最小化成功互获得的。前面这个是一个正则参数，主要作用是用来提高泛化能力，后面这个是误差，最小化他们的和即得到最优的结果。
p和x是之前通过模糊规则以及原始数据变化后得到，最后y是c维向量，即各个特征的结果表示。具体来讲c个值表示x属于各个特征的隶属度，如果第j个值为1其他值为0，表示xi属于第c类。相当于预测结果x属于哪一类，y是真实类别，最小化两者差距。
为了得到最优的参数pg，最小化目标函数，我们定义初始导数为0则优化结果p可以表示为

这样一个TSK模糊系统久搭建完了。其余参数回到前面去看，FCM算法以及高斯径向基函数。
代码块

import numpy as np
import math
from sklearn.preprocessing import OneHotEncoder
from sklearn.cluster import KMeans


class BaseCluster:
    def __init__(self, n_cluster, m, scale=1.):
        self.n_cluster = n_cluster
        self.m = m
        self.scale = scale

    def _euclidean_distance(self, X, Y=None):
        """
        return the element-wise euclidean distance between X and Y
        :param X: [n_samples_X, n_features]
        :param Y: if None, return the element-wise distance between X and X, else [n_samples_Y, n_features]
        :return: [n_samples_X, n_samples_Y]
        """
        if Y is None:
            Y = X.copy()
        Y = np.expand_dims(np.transpose(Y), 0)
        X = np.expand_dims(X, 2)
        D = np.sum((X - Y)**2, axis=1)
        return np.sqrt(D)

    def predict(self, X, y=None):
        """
        predict membership grad using fuzzy rules
        :param X: [n_samples, n_features]
        :param y: None
        :return: Mem [n_samples, n_clusters]
        """
        X = np.array(X, dtype=np.float64)
        if y is not None:
            y = np.array(y, dtype=np.float64)

        assert hasattr(self, 'variance_'), "Model not fitted yet."
        assert hasattr(self, 'center_'), "Model not fitted yet."
        d = -(np.expand_dims(X, axis=2) - np.expand_dims(self.center_.T, axis=0))**2 \
            / (2 * self.variance_.T)
        d = np.exp(np.sum(d, axis=1))
        d = np.fmax(d, np.finfo(np.float64).eps)
        return d / np.sum(d, axis=1, keepdims=True)

    def fit(self, X, y=None):
        raise NotImplementedError('Function fit is not implemented yet.')

class FuzzyCMeans(BaseCluster):
    def __init__(self, n_cluster, scale=1., m='auto', error=1e-5, tol_iter=200, verbose=0):
        """
        Implantation of fuzzy c-means
        :param n_cluster: number of clusters
        :param m: fuzzy index
        :param error: max error for u_old - u_new to break the iteration
        :param tol_iter: total iteration number
        :param verbose: whether to print loss infomation during iteration
        """
        self.error = error
        self.tol_iter = tol_iter
        self.n_dim = None
        self.verbose = verbose
        self.fitted = False


        super(FuzzyCMeans, self).__init__(n_cluster, m, scale)

    def get_params(self, deep=True):
        return {
            'n_cluster': self.n_cluster,
            'error': self.error,
            'tol_iter': self.tol_iter,
            'scale': self.scale,
            'm': self.m,
            'verbose': self.verbose
        }

    def set_params(self, **params):
        for p, v in params.items():
            setattr(self, p, v)
        return self

    def fit(self, X, y=None):
        X = np.array(X, dtype=np.float64)
        if y is not None:
            y = np.array(y, dtype=np.float64)

        if self.m == 'auto':
            if min(X.shape[0], X.shape[1]-1) >=3:
                self.m = min(X.shape[0], X.shape[1]-1) / (min(X.shape[0], X.shape[1]-1) - 2)
            else:
                self.m = 2

        N = X.shape[0]
        self.n_dim = X.shape[1]

        # init U
        U = np.random.rand(self.n_cluster, N)

        self.loss_hist = []
        for t in range(self.tol_iter):
            U, V, loss, signal = self._cmean_update(X, U)
            self.loss_hist.append(loss)
            if self.verbose > 0:
                print('[FCM Iter {}] Loss: {:.4f}'.format(t, loss))
            if signal:
                break
        self.fitted = True
        self.center_ = V
        self.train_u = U
        self.variance_ = np.zeros(self.center_.shape)  # n_clusters * n_dim

        for i in range(self.n_dim):
            self.variance_[:, i] = np.sum(
                U * ((X[:, i][:, np.newaxis] - self.center_[:, i].transpose())**2).T, axis=1
            ) / np.sum(U, axis=1)
        self.variance_ *= self.scale
        self.variance_ = np.fmax(self.variance_, np.finfo(np.float64).eps)

        return self

    def _cmean_update(self, X, U):
        old_U = U.copy()
        old_U = np.fmax(old_U, np.finfo(np.float64).eps)
        old_U_unexp = old_U.copy()
        old_U = self.normalize_column(old_U)**self.m

        # compute V
        V = np.dot(old_U, X) / old_U.sum(axis=1, keepdims=True)

        # compute U
        dist = self._euclidean_distance(X, V).T  # n_clusters * n_samples
        dist = np.fmax(dist, np.finfo(np.float64).eps)

        loss = (old_U * dist ** 2).sum()
        dist = dist ** (2/(1-self.m))
        dist = np.fmax(dist, np.finfo(np.float64).eps)
        U = self.normalize_column(dist)
        if np.linalg.norm(U - old_U_unexp) < self.error:
            signal = True
        else:
            signal = False
        return U, V, loss, signal

    def normalize_column(self, U):
        return U/np.sum(U, axis=0, keepdims=True)

    def __str__(self):
        return "FCM"

    def fs_complexity(self):
        return self.n_cluster * self.n_dim


class ESSC(BaseCluster):
    """
    Implementation of the enhanced soft subspace cluster ESSC in the paper "Enhanced soft subspace clustering
    integrating within-cluster and between-cluster information".
    """
    def __init__(self, n_cluster, scale=1., m='auto', eta=0.1, gamma=0.1,
                 error=1e-5, tol_iter=50, verbose=0, init='kmean', sparse_thres=0.0):
        """
        :param n_cluster:
        :param scale:
        :param m:
        :param eta:
        :param gamma:
        :param error:
        :param tol_iter:
        :param verbose:
        :param init:
        :param sparse_thres: percentile for dropping attributes
        """
        super(ESSC, self).__init__(n_cluster, m, scale)
        assert gamma > 0, "gamma must be larger than 0"
        assert (eta < 1) and (eta > 0), "eta must be in the range of [0, 1]"
        self.eta = eta
        self.gamma = gamma
        self.error = error
        self.tol_iter = tol_iter
        self.n_dim = None
        self.verbose = verbose
        self.fitted = False
        self.init_method = init
        self.sparse_thres = sparse_thres

        self.U, self.weight_, self.center_, self.v0 = None, None, None, None

    def get_params(self, deep=True):
        return {
            'n_cluster': self.n_cluster,
            'scale': self.scale,
            'eta': self.eta,
            'gamma': self.gamma,
            'max_iter': self.tol_iter,
            'error': self.error,
            'verbose': self.verbose,
        }

    def set_params(self, **params):
        for p, v in params.items():
            setattr(self, p, v)
        return self

    def fit(self, X, y=None):
        """
        :param X: shape: [n_samples, n_features]
        :param y:
        :return:
        """
        X = np.array(X, dtype=np.float64)
        if y is not None:
            y = np.array(y, dtype=np.float64)
        if self.m == 'auto':
            if min(X.shape[0], X.shape[1]-1) >=3:
                self.m = min(X.shape[0], X.shape[1]-1) / (min(X.shape[0], X.shape[1]-1) - 2)
            else:
                self.m = 2
        self.n_dim = X.shape[1]
        self.v0 = np.mean(X, axis=0, keepdims=True)  # v0: data center
        self.weight_ = np.ones([self.n_cluster, self.n_dim])/self.n_dim
        if self.init_method == 'random':
            self.center_ = X[np.random.choice(np.arange(X.shape[0]), replace=False, size=self.n_cluster), :]  # init by k-mean
        elif self.init_method == 'kmean':
            from sklearn.cluster import KMeans
            km = KMeans(n_clusters=self.n_cluster)
            km.fit(X)
            self.center_ = km.cluster_centers_
        else:
            raise ValueError('init method only supports [random, kmean]')
        # self.visual(X)

        loss = []

        old_V = None

        for i in range(self.tol_iter):
            self.U = self.update_u(X, self.weight_, self.center_, self.v0)  # U: [n_clusters, n_samples]

            self.center_ = self.update_v(X, self.U, self.v0)
            self.weight_ = self.update_w(X, self.U, self.center_, self.v0)

            loss.append(self.overall_loss(X))
            if i >= 1 and math.sqrt(np.sum((self.center_ - old_V)**2)) < self.error:
                break
            else:
                old_V = self.center_.copy()
        self.loss = loss
        self.fitted = True

        self.variance_ = np.zeros(self.center_.shape)  # n_clusters * n_dim

        for i in range(self.n_dim):
            self.variance_[:, i] = np.sum(
                self.U * ((X[:, i][:, np.newaxis] - self.center_[:, i].T) ** 2).T, axis=1
            ) / np.sum(self.U, axis=1)
        self.variance_ *= self.scale
        self.variance_ = np.fmax(self.variance_, np.finfo(np.float64).eps)
        self.norm_sparse_thres = np.percentile(self.weight_.reshape([-1]), self.sparse_thres)

        return self

    def update_u(self, X, W, V, v0):
        v = np.expand_dims(V.T, axis=0)
        x = np.expand_dims(X, axis=2)
        d1 = np.sum((x-v)**2 * W.T, axis=1)
        d2 = np.sum((V-v0)**2 * W, axis=1)[np.newaxis, :]
        d2 = np.repeat(d2, repeats=d1.shape[0], axis=0)
        d2 = np.fmax(d2, np.finfo(np.float64).eps)
        min_eta = np.min(d1/d2, axis=1, keepdims=True)
        self.sample_eta = np.minimum(self.eta, min_eta)
        d = d1 - d2*self.sample_eta
        d = np.fmax(d, np.finfo(np.float64).eps)
        d = d ** (1/(1-self.m))
        d = d / np.sum(d, axis=1, keepdims=True)
        return d.T

    def update_v(self, X, U, v0):
        u = np.expand_dims(U**self.m, axis=2)
        wd = u * (X - self.sample_eta * v0)[np.newaxis, :, :]  # weighted distance
        wd = np.sum(wd, axis=1) / np.sum(u*(1-self.sample_eta), axis=1)
        return wd

    def update_w(self, X, U, V, v0):
        v = np.expand_dims(V.T, axis=0)
        x = np.expand_dims(X, axis=2)
        u = np.expand_dims(U.T, axis=1)
        d1 = np.sum((u**self.m) * ((x-v)**2), axis=0).T
        u = np.expand_dims(U.T, axis=2)
        d2 = np.sum((V-v0)**2 * (u**self.m), axis=0)
        sig = d1 - self.eta * d2
        sig = -sig/self.gamma
        sig = np.fmin(sig, 700)  # avoid overflow
        # print(np.max(sig))
        sig = np.exp(sig)
        sig = np.fmax(sig, np.finfo(np.float64).eps)  # avoid underflow
        sig = np.fmin(sig, np.finfo(np.float64).max)
        return sig / np.sum(sig, axis=1, keepdims=True)

    def overall_loss(self, X):
        v = np.expand_dims(self.center_.T, axis=0)
        x = np.expand_dims(X, axis=2)
        d1 = np.sum(np.sum((x-v)**2 * self.weight_.T, axis=1) * (self.U.T**self.m))
        return d1

    def __str__(self):
        return "ESSC"

Mul-view TSK fuzzy system

以前多个特征提取和分类的步骤以及缺点：首先，使用不同的特征提取方法获得的不同数据集训练多个系统的计算成本较大。其次，很难确定哪种特征提取方法最适合于给定的癫痫性脑电图识别任务。第三，当基于不同视图的数据集对多个系统进行训练时，还不清楚如何从不同系统产生的多个决策中得出最终结论。
传统的多视角特征分类本文提出的模型框架。通过多视图学习过程，可以评估不同视图的重要性，并可以很容易地获得多视图决策。

基于香农熵的视角加权机制

主要是这个公式：这里的xgik表示第k个视角的第i个样本值，y是c维标签向量的第i个样本，c是所有的分类数。pgk表示第k个视角的pg，w表示第k个视角的权重，k是所有视角数（在这里是3），N是样本数，入w是香农熵的正则参数。通过以上两个公式确定wk，此时的入w也很重要，可以手工调整和交叉验证。

Multi-view collaborative learning mechanism

这里的主要公式如上，
关于解释没有太多需要自己理解的，读完即可。主要是先验知识的不同，这里先验知识需要FCM算法获得。通过最小化7.a，同理，也需要交叉验证或者手工调整

Objective function

关于以上的关系解释
对8-a采用条件极值，即拉格朗日乘除法，添加入为了解决优化问题，采用交叉迭代策略，这在FCM中也有使用，迭代过程包含两个步骤，更新p和更新权重w，采用控制变量法。
第一步：固定权重，求此时的pk极值
求偏导数，求此时的值。
固定pk求此时的w，更新规则如上。