2020年华为杯数学建模竞赛C题论文和代码

面向康复工程的脑电信号分析和判别模型

    摘           要：

脑电信号的识别和分类是脑机接口技术中非常重要的一环，使用者无需通过复杂训练就可以获得较高的识别准确率，具有稳定的锁时性和高时间精度特性。本文使用基于监督学习的随机森林，SVM算法，半监督学习S3VM 等算法，以较高的分类准确度，完成了对P300脑电信号的分析和判别。

针对问题一，我们设计了0.1~20Hz，阻带增益-80dB的46阶带通巴特沃斯IIR 滤波器对冗余数据做预处理，以是否出现P300脑电信号作为输入标签。考虑到发生P300事件和不发生的样本数不平衡，采用随机排列抽取法生成训练集，并建立了基于随机森林算法的监督学习分类模型。在python环境下搭建算法，预测5轮测试数据，采取“投票”机制得到不同受试者的10个测试结果(共48个目标字符)， 目标字符识别准确率为92%，相比KNN 模型84%的准确率，识别效果更好。

针对问题二，在问题一数据处理的基础上，从闪烁时刻起记录有效训练数据，并生成一个维度为7500×20的矩阵，同样地对P300 电位做标签化处理，得到一个7500维的向量。将通道选择问题转换成回归模型，进一步引入“ spike and slab”先验分布，建立贝叶斯稀疏回归模型。使用python3.8实现EP算法，推导出模型后验分布近似分布，依据支撑向量后验分布参数选择出最优通道组合，基于通道选择结果，再次采用问题一的算法和模型，找到了对所有被试者都适用的一组最优通道组合：Fz，F3，F4，C3，Cz，C4，CP3,CP4,CP5,CP6,P3,P4, P7。

针对问题三，在问题一数据处理和问题二通道选择的基础上，选择10个字符的数据作为训练集(共12×5×10×5=3000个样本), 其中前 30%作为有标签数据集,后70%作为无标签数据集，剩下作为验证集(共12×5×2×5=600个样本) 。建立半监督学习的S3VM 分类模型，并引入正负样本比例因子r排除正负样本数目不平衡的问题。模型通过带标记样本预训练分类器，用拉格朗日启发式方法对无标签样本进行标签预测，扩充训练集。通过模拟退火的最优化方法进行迭代，并在训练时对初始无标签样本和有标签样本赋予不同的权重。最终在验证集上的平均准确率达到75.3%，在测试数据中得到的结果基本吻合。

针对问题四，我们根据附件2所给的精确无噪声的样本标签数据，以及四个分区 Alpha, Beta, Theta, Delta所占能量强度不同的特点, 分别建立随机森林, 支持向量机(SVM) ，以及LightGBM 三种监督学习算法训练睡眠预测模型。在测试集和训练集的选取

上，采用“交叉验证法”，克服因为训练集测试集不同分割方法引起模型性能不稳定。测试选取不同比例训练集数据对模型预测精度影响。综合比较SVM，随机森林和LGBM算法预测精度和时间复杂度，得出，LGBM算法得到的精度为89%，且算法测试运行时间较短，其综合运行性能最优。

综上所述，本文建立的脑电信号分析和判别模型地实现了P300脑-机接口对目标字符的识别，最优通道选择，解决了带标签样本较难获得的问题。基于睡眠脑电信号频谱能量不同特点，建立模型实现了对睡眠状态的预测，有望成为评估睡眠质量、诊断和治疗睡眠相关疾病的重要辅助工具。

关键词：P300脑机接口； 监督学习； 随机森林算法； 通道选择； 贝叶斯稀疏回归模型；EP; spike and slab; 半监督学习S3VM 算法

神经康复和智能机器人等领域具有重要的研究意义的巨大的应用研发潜力。

P300是一种由小概率事件刺激后约300ms产生的诱发事件相关电位，是研究最多、应用最广泛的脑电信号，具有特定的目标性和可解释性，也是测试认知功能最常用的指标之一。P300 电位作为一种内源性成分，它不受刺激物理特性影响，与知觉或认知心理活动有关，与注意、记忆、智能等加工过程密切相关。基于P300的脑-机接口优点是使用者无需通过复杂训练就可以获得较高的识别准确率，具有稳定的锁时性和高时间精度特性，因此，P300-BCI是最适合严重残疾患者独立长期在家庭环境下做康复训练使用的BCI系统。由于P300脑机接口存在字符分类识别正确率不高以及信息传输率较低的问题，所以创建合适的模型算法来提高系统的精度、速度和能力是研究的核心目标。

1.2 问题提出

问题一：数据处理与目标分类问题。在考虑目标的分类准确率并保证信息传输速率的情况下，根据附件1所给数据，使用小于等于5轮次测试数据，找出附件1中5个被试测试集中的9或10个待识别目标，并给出具体的分类识别过程，可与几种方法进行对比说明设计方法的合理性。

问题二：数据选择与通道优化问题。设计一个通道选择算法处理原始冗余信息和通道数据，在满足通道组合数量的基础上给出针对每个单独被试、更有利于分类的通道名称组合。基于通道选择的结果，进一步分析对于所有被试都较适用的一组最优通道名称组合，并给出具体分析过程。

问题三：半监督学习分类问题。为了减少训练时间，请根据附件1所给数据，选择适量的样本作为有标签样本，其余训练样本作为无标签样本，在问题二所得一组最优通道组合的基础上，设计一种学习的方法，并利用问题二的测试数据(char13-char17)检验方法的有效性，同时利用所设计的学习方法找出测试集中的其余待识别目标(char18-char22) 。

问题四：睡眠分期预测问题。根据附件2中所给的特征样本，用尽可能少的训练样本设计一个睡眠分期预测模型，并得到相对较高的预测准确率，给出训练数据和测试数据的选取方式和分配比例，说明具体的分类识别过程，并结合分类性能指标对预测的效果进行分析。

2. 问题假设

(1) 假设五名被测试成年人身体健康，无精神疾病史、脑部疾病史且视力正常；

(2) 假设设备确定能够采集到P300 电位

(3) 假设两次目标刺激产生的正电位波峰不会重合；

(4) 假设个体间的差异不会超过设定值；

3. 符号说明

Pi,j	判断是否出现P300 电位的输入值
E	波峰出现的区间集合
Fs	采样频率
ti,j	第i个目标字符在第j轮的起始采样时间
d₁	20个通道一次采样数据的样本向量
tk	闪烁起始时刻
N	一组样本中每个通道的采样点数
D	维度为7500×20的有效数据矩阵
z	支撑向量
μ	支撑向量的后验分布不等于0的概率
{xᵢ,yᵢ}ι=1	带标记的脑电数据样本集
{x₃}x=l+1	未标记的脑电数据样本集
l	有标记样本集大小
u	无标记样本及大小
k	第二问中选取的k个最优通道
N	一组样本中每个通道的采样点数
r	正负样本个数比
xi	第i行睡眠特征数据或通道权重
φ(x₂)	映射后的特征向量

4. 问题一的模型建立与求解

4.1问题一的分析

从问题的描述可知，在小概率刺激发生后 300 毫秒范围左右 P300 将会出现一个正向的波峰(根据P300的时域特征)。考虑到个体间的差异性，我们首先假设P300的波峰特征出现在目标符号所在行列闪现后的300毫秒左右的可偏差范围内。本题在预处理不同通道采集信号的基础上，结合是否在规定区间内出现P300脑电信号波峰特征为特征向量(分类标准)，考虑到训练数据种类的不平衡，需要平衡P300和非 P300 电位的数据样本，并以此形成训练集并构建决策树，训练随机森林算法。使用算法预测 P300 信号出现的行和列，通过行列位置最终确定出S1，S4，S5的10个和S2，S3的9个待识别字符。

4.2 模型建立

4.2.1 数据降噪

本题附件给出数据是由传感器采集的原始的信号，由于电子元器件自身噪声，环境的影响以及被试者动作状态的影响，比如眨眼，皱眉头等，都会给采集数据带来一定的干扰，这些干扰大多为高频信号干扰。不正确的处理干扰数据会给接下来的模型训练带来极大的

影响。一般情况下，人体P300脑电信号的能量集中在0.1Hz 到20Hz频带内，将原始数据经过通带为0.1Hz到20Hz的带通滤波器既可以滤除掉大多数高频噪声并且保留有效信息。我们基于 Matlab2018b 设计通带为 0.1Hz 到20Hz,阻带增益-80dB 的 46 阶带通巴特沃斯IIR 滤波器，其幅频特性如图4-1所示：

选取1号被试者要求注视字符‘B’时，1号通道第250个采样点往后250个采样序列为例，图4-2为未经滤波的信号，图4-3为经过带通滤波器滤除干扰后的信号。可以明显的看出原始信号“杂乱无章”的干扰被去除，经过滤波后的信号在第46个采样点附近出现疑似P300信号，查阅1号测试者的 train   event数据发现，第46个采样点附近正好是字符‘B’的所在行闪烁的时候，此刻被试者1大脑中产生P300信号。说明设计的滤波器很好的滤除掉了干扰，保留了有效信息。

4.2.2 冗余数据去除

滤波虽然可以消除随机干扰带来的影响，但题目所给的数据是由20通道，250Hz的采样频率得到，数据量仍然十分庞大，原始数据中必然包含大量的冗余信息，冗余信息不仅会加大系统传输处理数据的负担，甚至会对后期字符预测模型的训练产生严重影响，降低模型的预测精度。题干中描述在小概率刺激发生后 300 毫秒范围左右 P300 将会出现的一个正向波峰，考虑到个体间差异性，且为了保证能够完整包含该P300波峰特征成分，假设 P300 的波峰特征出现在目标符号所在行列闪现时刻tA后的( tₖ500+tₖ区间内，又采样点数和时间间隔成正比，即

N= Fs·△t

即有效信号仅会发生在目标所在行和列闪烁后的125个采样点内，而目标信号闪烁之前人脑并不会产生 P300 信号，因此采样闪烁之前的信号是没有意义的，不仅增加了数据冗余，还会干扰随后预测模型的精度。当然，我们不仅要采集目标字符出现后125个采样点这些有效数据，也需要非目标字符所在行或列闪烁后的125个采样点后的信号，这个数据作为人体不出现 P300 信号时候的正常脑电波信号，以此作为负的训练样本。每个行和列的闪烁时间为80ms，闪烁后间隔80ms进行下一次闪烁，因此数据选取窗口的滑动点数大小应为 250Hz*80+80ms=40,数据的选取如下图所示：

4.2.3 数据标签的制作

预处理的另外一个重要的问题就是标签的制作。使用4.2.1以及4.2.2数据预处理及冗余去除过程后，我们仅仅需要判断当一个行或列闪烁时所采集的信号，判断 P300 信号是否出现，这里P300信号只有两种状态，有和无，所以将脑电信号的输入标签化处理：

区间E中出现波峰特征

10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

区间E中没有波峰特征

其中，P₄，用于判断是否出现P300 电位， E∈ti,jti,j+500,ti,j表示第i个目标字符在第j轮的起始采样时间。

4.2.4基于随机森林的分类问题算法

随机森林( Random Forest，RF) 算法是一种组合分类器算法，包含许多独立的决策树( Decision Tree, DT) [6]。决策树算法以树状结构表示数据分类的结果，是一种典型的数据分类器。决策树是随机森林的核心部分，依据条件熵和信息增益保证决策树的左右生成，其基本原理为[7]：

(1)采用 Bagging的方法从样本数N的数据集D中随机地选取不同变量数，用以形成众多的随机决策树并决定每棵树的分类节点。

(2)计算数据集D的信息熵，以及各个变量Xᵢ和D之间的条件熵。信息熵H(D)体现了数据发生的不确定性，信息熵的值越高，表明越混乱，分类效果不好。条件熵H(D|X₄)是在数据集D发生的前提下，X；变量所带来的熵值变化。用公式可表示为：

                         HD=E-logDi=-∑i=1nDilog2D,                     (2)

                              HD|Xᵢ=HDXᵢ-HXᵢ                     (3)

其中，H(D，X₃)为数据集D同变量之间的联合熵：

                         HDXi=-∑i=1n∑j=1npDXilogDXi                   (4)

(3) 计算信息增益g(D，X₄)，以表示数据不确定性的减小程度，计算公式为：

                            gDXᵢ=HD-HD|Xᵢ                         (5)

(5) 以信息熵为度量生成熵值下降最快的树，同时选择信息增益率和 Gini系数最大的变量作为决策树分裂节点的变量。Gini系数的计算公式为：

                                Ginip=1-∑k=1n|Ck||D|2                        (6)

(6) 每一棵决策树会根据待分类的数据给出一个基于该分类树的分类结果，并最终汇总全部森林里重复度最高的决策树结果(取众数) 作为随机森林算法的分类结果。

用随机森林算法解决分类问题的流程如下：

4.3模型的求解

4.3.1平衡训练数据种类

使用4.1预处理数据训练模型时，存在一个非常严重的问题。每个被试者进行的每个字符的每轮实验中，每个行和列都会闪烁一次，即一共闪烁12次，但在每轮实验中，出现目标字符的行和列一共只有2个。因此，经过预处理后标签为1的数据，即出现 P300 波形的数据仅有16.7%，其余84%的数据的标签均为0，即没有P300信号。用这样一个标签极其不平衡的数据训练的模型，最后模型必然出现预测结果大部分为0的情况。模型只要预测结果都为0，训练精度就能达到 84%。

为了解决不同标签数据量不平衡问题，我们采用随机排列抽取法来降低标签为0 数据量，保留标签为1数据量的方法。具体做法如下：每个被试者在进行同一个字符的5轮实验中，每轮实验目标字符所在行和列闪烁时所采集的数据，即标签为1的数据保留。剩下的行或列闪烁时，在每轮实验均分别取两个不重复的行或列作为有效训练数据。比如：第一个被试者在注视字符‘B’时，5轮实验中，每次行1或者列8闪烁时的数据保留，然后取其采样点后 125个数据作为有效数据，经过滤波，标签为1。如预处理部分所述，共得到10组标签为1的有效信号。剩下共10个行或者列闪烁时分别在不同组取，如2，3行闪烁时采集的标签为0的信号取与实验1，行4，5闪烁时候，经过预处理后标签为0的数据取于实验 2，以此类推。其余字符闪烁时也做同样的处理。这样每个实验者，在进行每个不同字符的实验时候，均能得到 10个标签为0 和标签为1 的有效训练数据。它们数量相等，不会出现直接将所有实验数据做同样处理后出现的不平衡情况。

另外，除有效字符所在位置外，余下的 10 个行和列在 5 个实验中的取法共有 C102C82⋅⋅⋅C42种方法，并且不同的字符进行不同取法进行组合，会产生大量不同的数据集。为了消除不同取法所带来的随机影响，每次采用不同取法整合出标签0的数据，与标签为1 的数据合成一个有效训练数据训练模型，将每种取法均实验一次显然是不可行，也是没有必要的。我们采取下列这样方法：将实验字符所在行或者列外的其他10个行或列随机重排两次，每次依次从排列的行或列中选择2个，这两个行和列的数据按顺序采集于实验1，2，…。比如1号被试者在进行字符‘B’实验时，将‘B’所在行或列去除，剩余10个行或列的一次随机排列为：

9, 5, 3, 7, 2, 4, 12, 10, 11

这些行或者列闪烁后采集数据标签一定为0，其中列9行5 闪烁时的数据采集于实验1，行3列7采集于实验2，以此类推。这样的取法每个被试者每个数据采集2次。所以每个被试者均可以得到两组有效训练数据。

4.3.2实验结果

基于 Python3.8搭建出随机森林字符预测模型，使用上述的数据处理得到训练数据训练模型，对测试集数据做同样的预处理后，送入训练出的模型进行预测计算。最为理想的数据预测情况应该是每个被试者的每个字符，行和列均出现一个P300信号，取出行和列，即可定位出预测字符。但模型不可能具有 100%精确度，有可能在其他位置也预测出 P300信号，我们采用“少数服从多数”的投票方法进行结果的提取。每个字符的5组测试数据中，模型预测出出现P300信号位置，就给该位置加1，即投给该位置一票，5组预测实验结束后，分别取得票最高的行和列作为最终结果，并且以此定位字符作为预测的最终结果。五个被试者的“投票”表格情况如下：

表 4-1 被测者 S1 P300信号预测位置统计表

char 行列	13	14	15	16	17	18	19	20	21	22
1	0	3	1	0	0	0	1	0	3	1
2	0	1	0	0	5	0	2	0	0	0
3	2	0	0	0	0	0	0	1	0	0
4	0	0	0	0	0	2	0	2	1	0
5	0	0	0	1	3	0	0	0	0	1
6	0	0	2	0	0	0	0	0	0	2
7	1	0	2	0	0	0	1	0	1	0
8	0	0	0	0	0	1	0	0	0	0
9	2	1	0	0	3	0	0	1	0	0
10	0	0	1	1	1	0	0	0	0	0
11	0	0	1	0	2	0	3	0	0	0
12	0	2	0	0	0	0	0	2	0	1

表4-2 被测者 S2 P300信号预测位置统计表

Char 行列	13	14	15	16	17	18	19	20	21	22
1	0	3	1	0	0	0	1	1	3	0
2	0	1	0	0	4	0	2	0	0	0
3	2	0	0	0	0	0	0	1	0	0
4	0	0	0	0	0	2	0	2	1	0
5	0	0	0	2	3	0	0	0	0	1
6	0	0	2	0	0	0	0	0	0	2
7	2	0	2	0	0	0	1	0	2	0
8	0	0	0	0	0	2	0	0	0	0
9	0	2	0	0	3	0	1	1	0	0
10	0	0	1	1	0	1	1	1	0	0
11	1	0	1	0	0	0	3	0	1	0
12	0	3	0	0	0	0	0	2	0	1

表 4-3 被测者 S3 P300信号预测位置统计表

Char 行列	13	14	15	16	17	18	19	20	21	22
1	0	2	1	0	0	0	0	1	2	0
2	0	0	0	1	4	1	2	0	1	1
3	3	0	0	0	0	0	0	0	0	0
4	0	1	0	0	0	2	0	2	0	0
5	0	0	0	2	3	0	0	0	0	1
6	1	0	3	0	0	0	0	0	0	2
7	2	0	2	0	0	0	1	0	2	0
8	0	0	0	0	0	2	0	0	0	0
9	0	0	0	0	3	0	1	1	1	0

10	0	1	3	1	1	1	1	0	0	0
11	0	0	1	0	0	0	2	0	0	0
12	0	3	0	0	0	0	0	2	0	1

表4-4 被测者 S4 P300信号预测位置统计表

Char 行列	13	14	15	16	17	18	19	20	21	22
1	0	2	0	0	0	0	0	0	3	0
2	0	1	0	1	4	0	2	0	0	0
3	3	1	0	0	1	0	0	0	0	0
4	0	0	0	0	0	3	0	1	1	1
5	0	0	0	2	2	0	1	0	0	1
6	0	0	2	0	0	0	0	0	0	2
7	2	0	2	0	0	0	1	0	2	0
8	0	1	0	1	0	2	0	0	0	0
9	1	0	0	0	2	0	0	1	1	1
10	0	0	1	2	0	1	1	0	0	0
11	0	0	0	0	0	0	2	0	0	0
12	0	2	0	0	0	0	0	2	0	2

表 4-5 被测者 S5 P300信号预测位置统计表

Char 行列	13	14	15	16	17	18	19	20	21	22
1	0	3	1	0	0	0	1	1	2	0
2	0	1	0	0	3	0	2	0	0	0
3	1	0	0	1	1	1	0	0	0	0
4	0	0	0	0	1	2	0	2	1	0
5	0	1	1	2	2	0	0	0	0	1
6	0	0	2	0	0	0	0	0	0	2
7	2	0	3	0	0	1	1	0	2	0
8	0	0	0	0	0	2	0	0	0	0
9	0	2	0	0	2	0	0	1	1	1
10	0	1	1	2	1	0	1	0	0	0
11	1	0	0	0	0	0	2	0	1	0
12	0	3	0	0	0	0	0	2	0	2

依据每个字符的5 轮实验中，模型预测的各位置出现 P300 信号次数，分别取出现次数最高的行和列序号，以此定位出需要识别的字符，其结果如下：

表4-6待识别字符统计表

被试者 字符	S1	S2	S3	S4	S5
Char13	M	0	M	M	M
Char14	F	F	F	F	F
Char15	5	5	8	5	5

Char16	2	2	2	2	2
Char17	I	I	I	I	I
Char18	T	T	T	T	T
Char19	K	K	K	K	K
Char20	x	x	X	x	x
Char21	A	A	A	A	A
Char22	0	/	/	0	0

在给出的五个测试结果中，随机森林算法对每个受试者的识别准确度情况如下表：

表4-7随机森林算法前五个测试数据精度表

受试者	S1	S2	S3	S4	S5
识别准确度	100%	80%	80%	100%	100%

4.3.3 简单的对比实验

为了检验随机森林算法的优越性，我们使用同样的数据训练最经典，最简单的机器学习模型KNN，采用同样“少数服从多数”方法，模型的字符预测结果如下：

表4-8 KNN模型字符预测表

被试者 字符	S1	S2	S3	S4	S5
Char13	M	N	M	O	M
Char14	F	F	F	F	F
Char15	Y	5	5	5	5
Char16	2	2	2	2	2
Char17	I	J	I	I	I
Char18	T	T	K	T	T
Char19	K	K	K	K	K
Char20	x	X	x	x	x
Char21	A	A	A	A	A
Char22	L	/	/	0	0

在给出的五个测试结果中，KNN算法对每个受试者的识别准确度情况如下表：

表4-9KNN算法前五个测试数据精度表

受试者	S1	S2	S3	S4	S5
识别准确度	80%	60%	100%	80%	100%

与随机森林法对比，KNN在应对高维二分类数据时性能较低，但其大部分结果与随机森林法相同。与题目提供的前五个字符的识别结果相比，随机森林法仅有两个结果预测错误，正确率高达92%，而KNN算法正确率仅为80%。在算法运行时间上看，KNN 是明显优于随机森林的。二者的正确率和运行时间如下表所示：

表4-10 随机森林与KNN对比表

算 法	KNN	随机森林
时 间(s)	0.0529	0.1730
精 度(%)	80	92

5. 问题二的模型建立与求解

5.1问题二的分析

题目给出的原始的20个通道采集的脑电数据非常庞大，其中参杂着很多冗余信息，不仅会带来系统的复杂性，增大系统数据传输和处理的压力，甚至会造成模型的分类和识别精度的下降，因此我们将剔除影响较小甚至没有影响的通道以选取最优的通道组合。传统的通道选择方法分为两类，一类是经验选择法，另外一种是机器学习法。经验选择法需要扎实的生物学相关的知识和前期大量的实验来论证各个通道的比例权重，而且其选择的算法不能满足每个被试者的特点，会造成重要数据丢失，影响实验结果。机器学习算法则是通过挖掘数据内在的特点，能够根据不同的被试者自动调整各个通道的采集数据的重要性。所以我们选择了机器学习算法来解决通道选择问题。

考虑到脑电信号子集的空间稀疏性， 目前已经存在多种稀疏贝叶斯回归方法，比如对所有脑电通道施加L₁范数约束，选择出具有空间稀疏性的最优通道子集 LASSO 方法[8], 对通道不同组合施加L₁范数约束的组GROUP-LASSO法[9]。本文以“ spike and slab”先验建立稀疏贝叶斯回归模型，以此训练出每个被试者的最佳通道组合。

5.2模型建立

不失一般性，我们以第一个被试者为例来推导我们的模型。首先，将某个采样点20个通道采集的数据记作d₁，d₄为一个长度为20的向量，即为给的表格数据中的一行。将每行或每列闪烁时刻起N个采样数据作为有效数据，每个被试者均进行5轮实验，每轮实验行和列一共闪烁12次，所有每个被试者有60N有效数据。和问题一保持一致，仍然选取N=125。将被试者所有的有效数据组合成一个矩阵，并且用字母D表示， D=d1d2d60kT.D是维度为7500×20的矩阵，在压缩感知领域也称为“字典”。同样地，我们也将P300信号进行标签化处理，记作y，y是维度为7500的向量。如果相应的采样时刻出现P300信号，那么y对应的位置为1，否则为0。20个通道的权重向量为x。那么权系数求解问题就转换成一个线性回归模型：

y= Dx+ε                                             (7)

为了方便模型的推导，将字典D的维度记作P×K.ε为具有精度为σ的高斯零均值的随机变量。为了得到关于x的稀疏解, 我们引入“ spike and slab”稀疏促进函数[5]。

                           x|z-∏j=1K1-zjδxj+zjNxj|0τ0                   (8)

其中x₁是K维权重系数向量x的第i个分量，他们是K个满足均值为0，方差为τ₀且互相独立高斯变量。z称为支撑向量，是K个i. i. d伯努利分布的二维随机变量。 zᵢ=1意味着 xᵢ≠0,，进一步意味着通道i 是需要选择的，是对睡眠状态判定有着至关重要的通道； zᵢ=0意味着 xᵢ=0则意味这条通道是多余的。δ(·)代表冲激函数。支撑向量满足的分布为

                                  z|π-∏j=1KBerπj                              (9)

根据贝叶斯准则，隐藏变量x，z的后验概率分布为

                        pxz|y=py|xpx|zpz∫py|xpx|zpzdxdz                    (10)

上式分母为一个归一化常数，不影响后验分布的具体形式。观察上式，分子为K维高斯分布与K维伯努利分布以及K维“ spike and slab”分布的乘积形式，想要得到其准确的闭式数学表达式是不可能的，进一步分母是分子的复杂的多元积分，分子闭式表示式的难以求解必然导致计算积分值也同样是不可行的，所以式(10) 是无法求得解析解的，即无法精确计算权重系数以及其支撑向量的精确后验分布。但在实际工程问题中，可以用来得到上式的近似解。常见的此类近似算法有MCMC，Gibbs sampler，VB等等，但上述算法均存在效率低，收敛情况难以检测等问题。本题中，我们采取期望传播算法解决上式( Expectation Propagation, 简称EP算法)。EP 算法高效, 并且收敛情况易于检测优点是我们选择的主要原因[4]。同样通过贝叶斯准则得到参数的后验分布天生具有防止过拟合的优点，这点是最大似然估计和最大后验准则所不具备的。

5.3模型的求解

 Expectation Propagation, 简称 EP 算法是一种稳定,快速的分布近似算法。它的核心思想是在指数分布族( Exponential Family) 概率分布中选取与似然函数耦合的分布，最小化两者的KL散度，从而求得模型的近似分布的解析表达式。将式(7) ~(9) 带入(10) 得到权重系数及其支撑变量的后验分布表达式为：

                           pxz|y=1Zpy|xpx|zpz                        (11)

 =1ZNy|Dxσ02I∏j=1κ1-zjδxj+zjNxj|0τ0p(z.

注意到似然函数p(y|x)为多维高斯分布，并且支撑集满足伯努利分布，而高斯-伯努利分布本身就是耦合的分布，由此可以推断出x，z的联合后验分布也一定是高斯-伯努利分布的。为了避免整体近似的复杂度，进一步观察上式，其可以分解为三项

f₁(x) = N(y| Dx, 21

 f2xz=∏j=1K1-zjδxj+zjNxj|0τ0=∏j=1Kf2jxjzj(12)

 f3zi=∏j=1KBerπj

第二项又可以进一步分解为K项，可以使用 EP 算法分别推导上三式的近似分布从而得到整体的近似解。EP算法的第一步就是选择合适的近似因子，近似因子的选择一方面要保证简单性，另外一方面又要保证模型具有足够的复杂度以能够很好的拟合数据内在的分布。注意到f₁只与x有关，f₂与x，z有关，f₃与z有关。针对本模型，我们选择3个属于指数分布族的概率分布形式如下：

代码

1.  function Hd = bandpass

2.  %BANDPASS Returns a discrete- time filter object.

3.  % Butterworth Bandpass filter designed using FDESIGN. BANDPASS.

5.  % All frequency values are in Hz.

6.  Fs = 250;  % Sampling Frequency

7.

8.  Fstop1 = 0.1;                           % First Stopband Frequency

9.  Fpass1 = 0.2;                           % First Passband Frequency

10. Fpass2 = 20;                           % Second Passband Frequency

11. Fstop2 = 30;                           % Second Stopband Frequency

12. Astop1 = 80;                           % First Stopband Attenuation (dB)

13. Apass  = 1;                            % Passband Ripple (dB)

14. Astop2 = 80;                           % Second Stopband Attenuation (dB)

15. match  = ' stopband';                     % Band to match exactly

16.

17. % Construct an FDESIGN object and call its BUTTER method.

18. h  = fdesign. bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, …

19.                                     Astop2, Fs);

20. Hd = design(h, ' butter', 'MatchExactly', match);

21.

22.% [EOF]

2.随机森林预测模型 ( python 3.8) :

1.  from sklearn. ensemble import RandomForestClassifier

2.  from sklearn. model   selection import train   test   split

3.  from sklearn. metrics import accuracy   score

4.  from sklearn. preprocessing import StandardScaler

5.  import nump y as np

6.  import scipy. io as sio

7.

8. #一号被试预测模型

9. #数据导入

10. file = sio. loadmat(' train. mat')

11. s1   tem = file['S1  train']

12.

13. f ile1 = sio. loadmat(' trainevent. mat')

14. s1   eve = file1['S1  train   event']

15.

16. file2 = sio. loadmat('testEVENT. mat')

17. s1   test   tem = file2['S1  test']

18.

19. file3 = sio. loadmat(' test. mat')

20. s1 _ test _ eve = file3['S1_ test _ event']

21.

22. # 训练1号被试者

23. dict={'0':[1,8],'2':[1,10],'4':[2,7],'6':[2,12],'8':[3,9],'10':[3,11],'12':[4,7],'14':[

4,10],'16':[5,8

24. ],'18':[5,12],'20':[6,9],'22':[6,11]}

25. labels= []

26. train _ data = np. zeros((60*12,125))

27. index = 0

28. for i in range(0,24,2): #列

59.     for j in range(66): #行

30.         h = s1 _ eve[j][i]

31.          if h<100:

32.              be = s1 _ eve[j][i+1]

33.              train _ data[ index,:] = s1 _ tem[ be-1: be+124,i//2]

34.              if h in dict[ str(i)]:

35.                  labels. append(1)

36.              else:

37.                  labels. append(0)

38.              index+=1

39.

40. labels= np. array( labels)

41. scaler = StandardScaler()# 标准化转换

42. scaler. fit( train _ data)# 训练标准化对象

43. train _ data= scaler. transform( train _ data)

44.

45. x_ train, x_ test, y_ train, y_ test = train _ test _ split( train _ data, labels)

46. clf = RandomForestClassifier()

47. clf. fit(x_ train, y_ train)

48. pre s = clf. predict(x_ test)

49. print( accuracy _ score( pres,y_ test))

50. #预测一号被试者

51. predict _ data = np. zeros((60*10,125))

52. index = 0

53. for i in range(0,20,2): #列

54.      for j in range(66):#行

55.         h = s1 _ test _ eve[j][i]

56.          if h<100:

57.              be = s1 _ test _ eve[j][i+1]

58.              predict _ data[ index,:] = s1 _ test _ tem[ be-1: be+124,i//2]

59.              index+=1

60.

61. scaler = StandardScaler() # 标准化转换

62. scaler. fit( predict   data)  # 训练标准化对象

63. predict   data= scaler. transform( predict   data)

64.

65. res1 = clf. predict( predict   data)

66. vote = np. zeros((10, 13))

67. index = 0

68. for i in range(0,20,2): #列

69.              for j in range(66): #行

70.                  h = s1   test   eve[j][i]

71.        ifh<100:

72.                         if res1[ index]:

73.                               vote[i//2,h] +=1

74.                         index+=1

3.基于EP的稀疏贝叶斯回归代码实现(python3.8)

1.  import nump y as np

2.

3.  import matplotlib. pyplot as plt

4.  import scipy. io as sio

5.  def reconstruct(A, y, P, K, sigma02, nvi2, tau0, alpha):

6.             ATA, ATy = np. dot(A. T, A), np. dot(A. T, y)

7.             J, H = ATA / sigma02, ATy / sigma02

8.             # Initialize Global Approximation

9.              mi, sigmai, nvi = np. zeros(A. shape[1]), np. ones(A. shape[1]), 0.5 *

10. np. ones(A. shape[1])

11.             mi   old = mi

12.             nvi   old = nvi

13.

14.            # Initialize Site Approximation

15.            mi2 sigmai2  inv, sigmai2  inv, nvi1 = np. zeros(A. shape[1]), 1e-6 *

16. np. ones(A. shape[1]), 0.5 * np. ones(K)

17.             for ittc s in range( max   itt):

18.                  # compute cavity

19.                   mi   sigmai   inv   bar, sigmai   inv   bar = mi / sigmai - mi2 sigmai2  inv, 1 /

20. sigmai - sigmai2  inv

21.                   nvi   bar = nvi2

22.

23.                  # positive cavity variance

24.                   index = sigmai   inv   bar < 0

25.                   update   index = sigmai   inv   bar > 0

26.                   sigmai   inv   bar[ index] = 1

27.

28.                  # Compute matching moments

29.         mi _ new, vi _ new, nvi1_ new

30. =

31.compute _ cs _ spike _ and _ slab _ moments( mi _ sigmai _ inv _ bar, sigmai _ inv _ bar,

32. nvi _ bar, tau0)

33.

34.        # Compute Site Updates

35.        sigmai2_ new _ inv = 1 / vi _ new - sigmai _ inv _ bar

36.        mi2_sigmai2_ inv _ new = mi _ new / vi _ new - mi _ sigmai _ inv _ bar

37.          idx = sigmai2_ new _ inv <= 0

38.         sigmai2_ new _ inv[ idx] = 1e-2

39.         sigmai2_ new _ inv = bound _ values(sigmai2_ new _ inv, MINIMUM_VARIANCE

40. , MAXIMUM_VARIANCE)

41.

42.        # adjust vi _ new and sigmai2_ new _ inv

43.         vi _ new[ idx] = 1 / (sigmai2_ new _ inv[ idx] + sigmai _ inv _ bar[ idx])

44.        mi2_sigmai2_ inv _ new[ idx] = mi _ new[ idx] / vi _ new[ idx]

45. - mi _ sigmai _ inv _ bar[ idx]

46.

47.        # update with damping

48.        nvi1[ update _ index] = (1 - alpha) * nvi1[ update _ index] + alpha *

49.  nvil _ new[ update _ index]

50.         sigmai2_ inv[ update _ index] = (1 - alpha) * sigmai2_ inv[ update _ index] +

51. alpha * sigmai2_ new _ inv[ update _ index]

52.         mi2_sigmai2_ inv[ update _ index] = (1 - alpha) *

53. mi2_sigmai2_ inv[ update _ index] + alpha * mi2_sigmai2_ inv _ new[

54.              update _ index]

55.

56.         # Update global approximation

57.          mi, sigmai, nvi = update _ global(A, P, K, sigma02, mi2_sigmai2_ inv,

58. sigmai2_ inv, J, H, nvi1, nvi2)

59.         # Check for EP convergence

60.          ep _ diff = np. linalg. norm( mi - mi _ old) / np. linalg. norm( mi)

61.          nvi _ dif = np. linalg. norm( nvi - nvi _ old) / np. linalg. norm( nvi)

62.          if ep _ diff < 1e-5 and nvi _ dif < 1e-5:

63.              break

64.          mi _ old = mi

65.          nvi _ old = nvi

66. return mi

67. def update _ global(A,P, K, sigma02, mi2_sigmai2_ inv, sigmai2_ inv, J, h, nvi1, nvi2,

68. diagnoal _ only= True):

69.      if P > K: # compute directly

70.         p1 = J + np. diag(sigmai2_ inv)

71.          mi, sigmai = np. linalg. solve(p1, h + mi2_sigmai2_ inv), np. linalg. inv(p1)

72.