由于本人对于脑机接口以及脑电技术的极度爱好（其实目的是：是把U盘插到大脑里，然后就不用学习了哈哈哈哈），近几月看了较多这方面的内容，变打算写下博客总结分析一下。

目录

一、  机器学习分析简介

二、机器学习分析的脑电特征

三、机器学习训练分析

四、机器学习分析的特征选择和降维

五、机器学习分析的选择分类器

六、机器学习分析的结果评估

 七、机器学习实例分析

机器学习和模式识别已被广泛应用于脑电信号分析领域，为从高纬度脑电信号提取和描述与任务相关的大脑状态提供了新的方法。鉴于近年来机器学习技术引发的持续关注及广泛的应用，我收集了近年来常见的利用机器学习技术处理脑电信号的案例，并分析总结出了本文内容。首先，我会带领大家分析并简单介绍机器学习的基本概念，接着引出一个关于区分睁眼闭眼状态下静息脑电的科学问题，并逐步介绍如何提取特征、训练特征、选择特征、降维分析、选择分类器、评估选择，最后进行模式表达。最后我们探讨了深度学习在未来研究中的应用。最后利用MATLAB进行了详细的脑电数据处理的介绍与应用。

基于脑电信号毫秒级别的时间分辨率,脑电已被广泛地应用于研究人脑中与感觉和认知相关的神经活动（Hu et al.，2011;Tu et al.，2014，2016;Zhang et al.，2012)。通过结合相关实验，对被试施加外界刺激或者要求其进行指定任务，脑与行为间的关联可以通过脑电反映出来。然而，脑电能在何种程度上实现“大脑解码”仍然是个疑问。例如，运用脑电来检测个体的行为状态(睁眼或是闭眼)准确率能够达到多少?我们能否运用脑电实时监测个体的认知程度?除开发脑电采集系统和优化实验设计之外,创新的数据分析方法是解决上述问题的关键。

近年来,使用机器学习技术来分析脑电信号引起了广泛的关注( Blankertz et al.， 2011;Makeig et al.，2012;Müller et al.，2008 )。越来越多的研究证据表明,机器学习可以从高维度且有噪声的脑电信号中提取有意义的信息。鉴于相关技术的关注度和广泛应用,本章我们将介绍如何使用机器学习来分析脑电信号，包括机器学习的方法和相关应用，并探讨深度学习算法在未来研究中的应用前景。

一、  机器学习分析简介

详细的机器学习请参考这篇博客：机器学习系列（一）——机器学习简介_zxhohai的博客-CSDN博客_机器学习
研究者设计神经成像实验区分大脑的不同状态,而从实验中采集到的脑电信号往往包含大脑中的认知或感觉反应。传统的分析方法大多运用回归分析来建立基于特定假设的脑电信号特征(例如，ERP及其与行为指标之间的关系)，或是通过单独分析每个样本，来确定时间域或频率域中的哪些特征参与了大脑中认知状态或感觉反应的处理。我们称这些传统方法为“单变量分析”。从理论上讲，如果大脑的某个特定时间点上的活动在两种状态下存在差异，那么就有可能利用这种脑活动来解码个体的认知状态或感觉反应。例如，α频段振荡的能量与睁眼和闭眼状态密切相关。然而,大多数情况下，我们往往很难找到一个能够提供足够大的差别来区分两种状态的单独特征，从而导致无法做出准确的判断。
机器学习分类器可以被看作一个函数,其将不同状态下的大脑活动的各种特征值作为自变量，并预测未知状态属于何类（因变量)。下面介绍关于机器学习的几个基本概念。
类:一个对象所属的类别。在一个类中，一组模式共享公共属性并且通常来自同一个源。
模式 : 一个对象的特征集合，以及该对象的类信息。样本:对象的任何给定的模式都称为样本。
特征 : 一组带有区分和鉴别一个对象的信息的变量。
特征向量 : 一个样本中K个特征的集合，以某种方式排列咸K维向量。特征空间:特征向量所在的K维空间。
在一个脑电实验中，特征可以从时间域、频率域或者空间域中被提取。类可以是认知状态或感觉反应。我们可以将N个样本的K个特征（例如，一次试验或一个被试为一个样本）指定为类标签。

二、机器学习分析的脑电特征

为了在脑电研究中使用机器学习分析，我们需要确定哪些脑电特征可以用来区分实验条件。对于自发性脑电数据，我们通常在频率域通过频率和幅值来研究信号。例如，闭眼状态的静息状态脑电数据α频段震荡的能量显著高于闭眼状态。

有一些特征既可以用来分析自发性脑电，也可以用来分析ERP。

 总之，选择合适的脑电信号特征对于机器学习分析至关重要。研究人员可以根据假设或者数据驱动搜索做出决定。

三、机器学习训练分析

 机器学习分类器的目标是基于通过训练获得的一些先验知识来估计给定特征向量对应的正确类别。训练是分类器学习特征向量与其对应的类标签之间映射关系的过程。特征向量和对应类标签之间的这种关系在特征空间中形成了一个决策边界,它使得不同类的模式彼此分开。
在训练分类器之前,我们需要将样本（特征向量和类标签)整理成为训练数据集（用于训练分类器的数据;已知数据类标签）和测试数据集(用于评估分类器性能的数据;通常与训练数据同时采集或从训练数据中划分;未知数据类标签)。在某些情况下，还会有一个独立数据集（采用被训练好的分类器对未知数据进行分类;未知数据类标签)。在训练分类器过程中,我们需要在使用更多的有噪声的样本(例如,单次提取的脑电信号)或更少的干净的信号(例如，平均信号)之间权衡。一方面，拥有更多的样本有助于分类器在训练中进行更好的参数估计,并提高测试的能力,以保证测试的准确性;另一方面，噪声会影响参数估计，在某些情况下会使得分类器的参数过度拟合噪声,影响分类器的性能。
使用交叉验证的方法可以利用所有数据进行训练和测试。

对于N个样本，我们可以把它们分为N-1个训练样本和一个测试样本，并重复相同的过程N次，使得每个样本都被用作测试样本一次。这种方法被称为“留一法”交叉验证(Tu et al.，2016)。尽管在严格意义上每次迭代中训练的分类器有所不同，但由于它们共享很多训练数据，所以可以预期它们是非常相似的。除了“留一法”之外,其他种类的交叉验证，例如“K折叠”[其中K是数据集被分割成的不同折叠的数量(例如,5折或者10折)]也已被广泛应用于脑电研空中。与“留一法”相比,“K折叠”交叉验证在计算上更有效，并且训练的分类器可能不会过度拟合。如下是十折交叉验证的示意图：

 MATLAB函数”crossvalind“可以帮助我们生成N个样本的”K折叠“交叉验证的索引该索引包括相等的从1到K的整数，用来将原来的N个样本分成K个子数据集。更详细的介绍如下：

 load fisheriris
 indices = crossvalind('Kfold',species,10);
 cp = classperf(species);
 for i = 1:10
    test = (indices == i); 
    train = ~test;
    class = classify(meas(test,:),meas(train,:),species(train,:));
    classperf(cp,class,test);
end
 cp.ErrorRate

matlab--交叉验证函数crossvalind_囊萤映雪的萤的博客-CSDN博客_crossvalind

 为了通过训练建立更好的分类器，研究者需要考虑训练的效果，即分类器正确识别训练数据类别的能力以及泛化的效果,即分类器识别测试数据类别的能力。一个好的机器学习分类器应具有一个能够提供最好的泛化能力而不是完美的训练效果的决策边界。在某些情况下,分类器因为学习了数据中的噪声使得训练数据完美地被分类。我们将这些分类器学习了噪声或是随机误差而不是真实的数据与类标签之间关系的现象称为过度拟合。

四、机器学习分析的特征选择和降维

脑电数据通常会跨越多个域具有极高数据维度,但样本量非常有限,机器学习在脑电研究中的效果往往会受到“维度诅咒”的影响(Mwangi et al.,2014)。例如，一个64通道的时频图(时域中有1000个样本，频率域中有100个样本）具有6400 000个特征（假设每一个时间-频率点为一个特征)。特征选择或降维对于从高维度脑电数据中识别一小部分判别特征,以获得更高的分类准确性和更好的分类器可解释性是至关重要的(Tu et al.，2015)。
特征选择或降维既可以使用类标签(称为“有监督学习方法”)，也可以不使用类标签(称为“无监督学习方法”)。主成分分析是一种传统的降维方法,其沿着一个方向投影高维度数据。在这个方向上,数据的方差由少量潜变量在一个线性子空间内实现最大化(Jolliffe,2011)。从数学上看,主成分分析运用正交变换将一组相关变量的观测数据转换为一组称为主成分的线性不相关变量。通过省略数据中低方差的主成分,我们仅丢失少量信息。假设我们只保留L(L小于变量数）的主成分，这样新数据就只有L列,但是包含大部分数据信息。通常情况下,在脑电分析中主成分分析可以用于减少时间、频率、通道、试次和被试的维度，因为脑电数据的这些域中往往包含冗余信息(Hu et al.，2015 )。
 

在MATLAB中，我们可以用函数“pca”来实现主成分分析。代码如下。

clear all; close all;
[coef, score, latent] - pca(X)

对于一个n行、k列的矩阵X(X的每一行代表样本，每一列代表特征)，这个函数的返回主成分系数为“coef ”、主成分评分为“score”（即X在主成分空间的表达)，以及主成分方差为“latent”(X的协方差矩阵的特征值)。在“coef”中,每一列代表了一个主成分的系数,这些列按照主成分的方差降序排列。在“score” 中每一行代表了一 一个样本，而每一-列代表了一个主成分。
然而，无监督学习方法不利用类标签来减少维度，不能够保证类在低维度空间中被很好地区分。相比之下，有监督学习方法利用类标签来确保高维数据可以被映射到低位空间，且不同的类可以在这个空间中被很好地区分( Mwangietal, 2014)。 在这里，我们简要介绍三类有监督的特征选择方法。
最直接的方法是使用单变量统计技术的过滤类方法，包括t检验、方差分析和皮尔逊相关性。这些方法根据在区分类别之间差异中特征的关联性来对特征进行排序，通常对于高维度脑电数据表现不佳。包装类方法使用分类或回归目标函数根据模型中的分类权重对特征进行排序。递归特征消除是最流行的包装类特征选择方法( Gysels et al., 2005 )。它从训练集中的所有特征开始，迭代地消除特征，直到找到最佳数量的特征。嵌人类方法通过对机器学习模型实施某些“惩罚”来选择特征。例如，最小绝对值收敛和选择算子(或称为套索算法LASSO),通过在回归系数上加上增强稀疏性的L1正则化来最小化均方误差，把较小的回归系数缩小到零来实现特征选择( Tibshirani,1996 )。在处理强相关的特征的情况下(例如，相邻的脑电时间点), LASSO算法从--组强相关变量中任意选择-个变量,即降低了机器学习模型的可解释性。因此，弹性网( elastic net)算法被开发并运用在脑成像研究中。弹性网算法通过结合L1 正则化(强制稀疏性)和L2正则化(分组约束)来提高LASSO算法的效果( Zou & Hastie, 2005 )。

五、机器学习分析的选择分类器

在准备完样本和选择信息量最大的特征之后，接下来是选择分类器(Lotte et al.，2007 )。学习特征和类标签之间的函数f的分类器分为判别模型和生成模型两种。判别模型的目标是通过设置参数来直接学习具有给定参数形式的预测函数。一个函数最简单的形式即是分类依赖于以下特征的线性组合:

Xw=x1w1+x2w2+...+xkwk

当Xw>0时，分类器判别为A类别;当Xw<0时，分类器判别为B类别。因此,学习一个线性分类器相当于学习一条能够将样本点区分成两个类别的线，我们称这条线为决策边界。线性判别分析（linear discriminant analysis,LDA）和线性支持向量机（support vector machine，SVM）是常用的学习决策边界的方法（图14.6) ( Subasi & Gursoy，2010)。

 线性判别分析通过找到一个投影方向使得不向类的所有样本都被很好地区分开。线性支持向量机在高维空间中构造超平面，以便在该超平面中实现最佳决策边界（极大化不同类之间的边界距离)。支持向量机的优点在于它能够使用核函数将原始特征空间（类不可以线性分离）映射到更高维空间（类可以线性分离)，并找到更高的支持向量维度空间。在使用支持向量机时，需要调整几个参数来实现一个良好的分类。对于线性支持向量机，需要在训练分类器之前定义成本参数（C)。这个参数告诉我们支持向量机在优化过程中希望避免每个训练样本被错误分类的程度。对于较大的成本参数值,优化过程会选择一个较小距离的超平面，这个超平面可以使得训练样本被正确分类。反之，一个较小的成本参数值会使得优化器去寻找一个较大边界的超平面，即使这个超平面错误地分类了很多样本。对于非线性支持向量机，通过内部交叉验证的网格搜索将有助于定义最佳的参数组合。
相比之下,生成模型,如高斯朴素贝叶斯(Gaussian Naive Bayes,GNB)分类器通过学习一个可以生成给定类别的样本的统计模型来进行分类(Huang et al.,2013)。其中,建模的是类条件特征值的分布，即 p(X|y=A4)和p(X|y=B),然后根据贝叶斯规则进行转换,通过确定哪个概率最大来进行分类。在实际操作中,朴素贝叶斯仅需要少量的训练数据来估计分类所需的参数便可以获得良好的分类性能。
分类器可以使用MATLAB内置函数“classify”或其他工具包来实现。使用“classifv”函数,我们可以运行以下 MATLAB代码:

clear all;
close all;
class=classify(sample,training,group,'tpye');

 这个函数将“sample”中的数据的每一行分类到“training”中的每一个“group”。“sample”和“training”必须是有相同列数的矩阵，“group”是用来定义“training”中每一行的样本分别属于哪一类的向量。我们可以通过指定“type”来使用不同种类的分类器,包括“linear”"diaglinear"" quadratic”“diagquadratic”“mahalanobis”。关于“classify”函数的详细说明，可以参考MATLAB的说明文件。
如果要运行支持向量机,推荐使用LIBSVM工具包( http://www.csie.ntu.edu.twl~cjlin/libsvm/ )。LIBSVM 集成了支持向量分类、回归和分布估计。在MATLAB中使用LIBSVM有以下几步:首先,下载和解压完工具包之后，需要将工具包添加到MATLAB路径中:主页一设置路径—添加并包含子文件夹—选择 LIBSVM 的路径。接下来需要选择并配置一个编译器，使用代码mex -setup。

六、机器学习分析的结果评估

训练结束之后,分类器的泛化效果需要在测试数据或者独立数据上进行评价。评价分类器的最直接标准是准确率，即被正确分类样本数在所有样本数中的比例。需要注意的是，准确率不适用于不平衡的数据集（两个类的样本数量有很大不同）。更重要的是，由于并非所有错误都有相同的成本，准确率并不能够反映某些分类错误的成本。例如,在疾病诊断中有两种类型的错误:①病人被误诊为健康人，可能导致其死亡;②健康人被误诊为病人，导致其服用不必要的药物。
因此,我们需要其他的标准来评价二分类的效果。阳性和阴性分别代表是否拒绝零假设。在典型的医学诊断中，真阳性代表病人被正确诊断为患病;假阳性代表健康人被错误诊断为患病;真阴性代表健康人被正确诊断为健康;假阴性代表病人被错误诊断为健康。
更多可以参考：理解准确率、精确率、召回率等评价指标含义以及在SVM模型中的应用_是故里吖的博客-CSDN博客_svm准确率

机器学习笔记14-----SVM实践和分类器的性能的评价指标(了解python画图的技巧) - 雨后观山色 - 博客园

此外,还有两个评价标准经常在实际中得以应用:灵取度衡量了阳性样本中判断为阳性的比例（真正患病的人被诊断为患病的比例)，计算方式是真阳性除以真阳性+假阴性的比值。特异性衡量了阴性样本中判断为阴性的比例(真正健康的人被诊断为健康的比例)计算方式是真阴性除以真阴性+假阳性的比例。
对于大多数分类器,使用相同的参数不可能同时实现最大的灵敏度和最大的特异性。在这种情况下，必须决定两者中哪一个更重要。研究这个问题的一个重要工具是接受者操作特征（receiver operating characteristic，ROC)曲线。在 ROC曲线中，灵敏度被绘制成一个参数的不同标准值的(1-特异性）的函数。ROC曲线下的面积( area under curve，AUC）是评价一个参数区分两个诊断组程度的衡量标准。
在衡量分类结果的时候,只提供准确率、灵敏度或者特异性的数值是不够的。统计上显著的分类结果是我们可以拒绝分类器随机决定的零假设。如果分类器没有关于使用脑电数据预测的类的任何信息,我们可以预期分类器只能随机猜测类标签。也就是说,如果测试样本是相互独立的,并且两个类的样本大小相等,那么基于二项分布,我们可以预期平均准确率为50%。因此,需要通过确定当满足零假设时所获得的准确率可能高于50%的概率来计算统计显著性。例如，我们可以使用针对50%的单样本t检验来表明分类准确率显著高于随机水平,或者我们可以使用配对t检验来比较不同特征集或分类器的性能。

然而,如果脑电数据的样本不满足独立性,则建议使用非参数置换检验来评估显著性(Maris & Oostenveld，2007)。在置换检验中，我们在训练之前随机置换数据的类标签。然后,对置换数据集执行交叉验证，并重复该过程1000次或更多次。如果在真实数据标签上训练的分类器的准确率统计上超过随机重新标记的数据标签训练的分类器的准确率（95%置信区间外),则认为分类器表现良好。

 七、机器学习实例分析

在本部分,我们将演示如何在一个简单的数据中进行机器学习分析。数据来源于一名被试在两种状态（睁眼和闭眼）下的在Oz电极上采集的静息态脑电信号。在睁眼和闭眼条件下，分别有60个数据作为训练样本，同时有另外60个未标记的数据作为测试样本。这个示例的目的是通过训练样本训练一个分类器来区分睁眼和闭眼条件，并在测试数据上进行标签识别。
在这个示例中,我们总结了以下概念。①类:睁眼和闭眼。②样本:脑电数据。③特征:α频段振荡能量。在前面的章节中,我们展示了在闭眼状态下静息态脑电的α频段振荡能量显著比睁眼状态下高。④特征矩阵:在0z电极上采集的信号的α频段振荡能量（维度为1)。⑤特征空间:α频段振荡能量（一维空间)。⑥模态:α频段振荡能量在闭眼状态下大于睁眼状态。
接下来,我们提供 MATLAB代码,包括:①在训练和测试样本中提取特征;②训练线性判别分析（ linear discriminant analysis, LDA）分类器，并用十折交叉验证评价训练效果;③在测试样本上评价LDA分类器的泛化效果。
第一步，读取 MATLAB文件“data_classification.mat”进入工作区。此文件中包含:①ec（60个闭眼状态下的脑电数据样本;每个样本包含2000个时间点);②eo （ 60个睁眼状态下的脑电数据样本;每个样本包含2000个时间点);3test_samples(60个测试样本;每个样本包含2000个时间点);④test_labels（ 60个二进制数字，1或者0;代表了测试样本的类标签)。脑电数据的采样率是200Hz,每个样本的采集时间长度是10s。代码如下：

clear all; close all;
 
load data_classification.mat
 
%% parameters
fs_test = 200; % sampling rate
N_Train = size(ec,2); % number of training trials
N_Test = size(test_samples,2); % number of test trials
 
 
%% PSD
nfft = 256; % Point of FFT
for n=1:N_Train
    [P_ec(:,n),f] = pwelch(detrend(ec(:,n)),[],[],nfft,fs_test); % calculate PSD for ec condition
    [P_eo(:,n),f] = pwelch(detrend(eo(:,n)),[],[],nfft,fs_test); % calculate PSD for eo condition
end
for n=1:N_Test
    [P_test(:,n),f] = pwelch(detrend(test_samples(:,n)),[],[],nfft,fs_test); % calculate PSD for test samples
end
 
%% feature extraction (alpha-band power in ec is significantly larger than in eo)
alpha_idx = find((f<=12)&(f>=8));  % frequency index of alpha band power
a_ec_train = mean(P_ec(alpha_idx,:)); % extract alpha band power from eo
a_eo_train = mean(P_eo(alpha_idx,:)); % extract alpha band power from ec
a_test = mean(P_test(alpha_idx,:)); % extract alpha band power from test data
 
%% 10-fold CV on training data
all_samples = [a_eo_train,a_ec_train]'; % all samples
all_labels = [ones(size(a_eo_train,2),1);zeros(size(a_ec_train,2),1)]; % labels of all samples: 1 for eo; 0 for ec
K = 10; % K-fold CV
indices = crossvalind('Kfold',all_labels,K); % generate indices for CV
for k = 1:K % K iterations
    cv_test_idx = find(indices == k); % indices for test samples in one trial of validation
    cv_train_idx = find(indices ~= k); % indices for training samples in one trial of validation
    cv_classout = classify(all_samples(cv_test_idx,:),all_samples(cv_train_idx,:),all_labels(cv_train_idx));
    cv_acc(k) = mean(cv_classout==all_labels(cv_test_idx)); % compute accuracy
    TP = sum((cv_classout==all_labels(cv_test_idx))&(cv_classout==1));
    TN = sum((cv_classout==all_labels(cv_test_idx))&(cv_classout==0));
    FP = sum((cv_classout~=all_labels(cv_test_idx))&(cv_classout==1));
    FN = sum((cv_classout~=all_labels(cv_test_idx))&(cv_classout==0));
    cv_sensitivity(k) = TP/(TP+FN); % compute specificity
    cv_specificity(k) = TN/(TN+FP); % compute sensitivity
end
cv_acc_avg = mean(cv_acc); % averaged accuracy
cv_sensitivity_avg = mean(cv_sensitivity);  % averaged sensitivity
cv_specificity_avg = mean(cv_specificity);  % averaged specificity
    
%% test on test data
% Concatenate training/test data and specify the labels
train_samples = [a_eo_train';a_ec_train']; %  training samples
train_labels = [ones(N_Train,1);zeros(N_Train,1)]; % labels of training samples: 1 for eo; 0 for ec
test_samples = [a_test']; % test samples
classout = classify(test_samples,train_samples,train_labels,'linear');
TP_test = sum((classout==test_labels)&(classout==1));
TN_test = sum((classout==test_labels)&(classout==0));
FP_test = sum((classout~=test_labels)&(classout==1));
FN_test = sum((classout~=test_labels)&(classout==0));
test_acc = sum(classout==test_labels)/N_Test; % compute accuracy
test_sensitivity = TP_test/(TP_test+FN_test); % compute specificity
test_specificity = TN_test/(TN_test+FP_test); % compute sensitivity

喜欢的小伙伴点赞加关注哦！！！