Kaggle比赛记：地震预测

本文前半部分是介绍Kaggle的机器学习比赛，后半则是记录今年3月到6月带领几位University of California, Irvine数学系本科生参加由美国洛斯阿拉莫斯国家实验室在Kaggle上组织的地震预测比赛。

本文并不涉及很多专业领域高深的数学公式或者机器学习模型，主要以介绍为主。一些术语的出现笔者也觉得应该不影响阅读，至于代码则可以跳过。

大型机器学习社区Kaggle

Kaggle是什么？

Kaggle是一个大数据机器学习的生态圈（现在已经被Google收购），举行机器学习比赛，同时有公有领域的很多数据集（比如UCI machine learning repo），云端即时运行的Python和R教程和kernels，还有类似于StackOverflow的找工作平台。

题外话：Kaggle的名字，自己猜测应该是源于gaggle（鹅叫象声词），因为所有新用户的默认头像是一只鹅……（不知是否创站者是《让子弹飞》的拥趸，把无知群众称作鹅）

Kaggle在开放云端运行Python和R程序的平台之前，让其名声大噪的比赛莫过于和CERN合作的寻找希格斯玻色子；和Kaggle同样名声大噪的还有古老的Gradient boosting，这次比赛也让人们知道，在后神经网络大数据时代的机器学习比赛中，用高度优化的树生成code，GB依然是山头一霸。

比赛

近几年，Kaggle的比赛根据其数据的类别主要有如下几类：

• 图像：比如自动驾驶，识别Youtube视频里的内容，卫星云图的分类，钢里面是否有瑕疵，等等等等。虽然Kaggle用户可以免费使用强大的Tesla P100 GPU，但现在每个GPU每周的使用时间不得超过30小时，所以本地如果有2080Ti这样的硬件支持，会有更多自由来实现自己的想法。
• NLP自然语言学习：比如Jigsaw举行的几次对话AI比赛，这类的比赛往往模型固定（比如BERT，U-net）,大家得分也相当接近，拿到奖牌需要一点运气。
• Tabular数据：这里主办者已经把数据中每个样本的特征给你生成好了，你需要用一个或者多个模型来根据这些特征进行某些预测；比如前两个月刚结束的IEEE组织的信用卡欺诈预测，根据每笔交易的特征来判断这笔交易是否是盗刷。还有Santander银行主办的用户回馈系列的几个比赛。
• 基于物理模型的数据：比如本文即将详述的地震预测比赛，参赛者需要根据物理（或者暴力计算）生成训练机器学习模型需要的特征（feature engineering），从对机器学习模型进行训练。

Kaggle自己也根据比赛举办的方式，把比赛分为featured（规模最大，比拼机器学习编程调参硬实力），research（会涉及某些领域的专业知识），和playground（比如一些有趣的破译密码比赛：Ciphertext Challenge II）。

论剑

大部分Kaggle的“监督学习”（supervised learning）比赛都会有如下几个元素：

• 数据：有一组训练数据（training data），训练集的目标（target）是已知的。另外有一个目标未知，需要参赛者使用模型进行预测目标的数据集（testing data）。

• 评分方式：主办方会设定一个metric，比如AUC（area under curve），MAE（mean absolute error），Dice coefficient，等等等；甚者还有组织者自己设定的metric，比如正在进行的NFL预测新赛季前进码数的比赛使用的是Continuous Ranked Probability Score。参赛者则需要参考这些metric来决定训练模型所需要的loss function。

• Leaderboard天梯：当你把通过训练集训练出的模型为测试集生成的目标预测上传到Kaggle之后，电脑会自动给你判分。接着就能即时看到你的答案和全世界其他优秀的data scientists相比有几斤几两了。为了避免作弊，测试集合又分成了两个不相交的集合（public和private），上交答案时候，你只能看到答案在public集合上的分数，作为现在leaderboard排名的依据；比赛结束的时候，每个参赛队需要从上交的几十上百个答案中选择两个，作为最终的答案去private集合上评分。开榜之后，这个分数会成为每队在leaderboard上的最后分数来决定队员是否能拿到奖牌。

• 还有少部分特殊形式的比赛，比如Two Sigma组织的股票预测，并不是一下子开榜。测试集的数据就是即时的股票交易数据，大家的模型上交之后，要接受每天的即时股票价格变化的检验而决定每个队最终的名次。

Kaggle比赛和应用数学物理研究不太一样。学术圈写论文的时候，假奇技淫巧，在同等条件下竟然能让竞争对手的方法看起来比自己类似的方法差不少。Kaggle最后开榜，大家一视同仁，90%靠的是机器学习绝对实力，10%还靠点运气（比例待商榷）。

预测地震

2019年Kaggle上规模最二大的比赛便是洛斯阿拉莫斯国家实验室（LANL），宾夕法尼亚州立大学（Penn State），和普度大学（Purdue）组织的“预测地震”。

地震的数据来源于Penn State Rock Mechanics laboratory岩石力学实验室，人工模拟岩石断裂带在不断受到压迫之后产生的剪应力（shear stress）造成的能量释放。

LANL两年前的研究中，使用了机器学习技术来预测受压迫力分布比较符合正态分布的quasi-periodic地震，成果斐然。

这次在Penn State的实验中，受力应该更加随机一些，而地震发生的频率也显得并无太多周期性（aperiodic）。其研究成果发表于Geophysics Research Letters上： Earthquake Catalog‐Based Machine Learning Identification of Laboratory Fault States and the Effects of Magnitude of Completeness，关于这篇论文对于这次比赛的重大影响，将在后记提到。

数据

训练数据集是有一整条6亿3千万个时间点的时间序列？！这次比赛算是自己真正接触“大”数据，相较之前玩过的机器学习toy模型，比如MNIST或者CIFAR几万几十万个样本相比其还真的只能算是“toy”。

这个时间序列其值是由一个piezoceramic sensor测得的电压，和地震波造成的总形变成正比，官方把这个数据以acoustic data命名，以下我们就称这个为时间（域）信号。

由于这个数据集实在太大了，用Pandas（Python Data Analysis Lib）读取的时候只能用int16精度（俗称half-precision，相比科学计算一般用的double precision，float64）：

train = pd.read_csv('../input/train.csv',
                    dtype={'acoustic_data': np.int16, 'time_to_failure':np.float32})

而另外一个变量time_to_failure，以下简称TTF，顾名思义是表示当下时间点距离下次地震还有多少时间。

采样百分之一之后，时间信号 vs TTF如下图：蓝线就是时间信号，蓝线的纵轴在左边；绿线是TTF，纵轴在右边。每次绿线归零之后，就是地震发生了。

放大来看，临近地震的时候时间信号如下图：大家可以发现，在每次地震之前，总会有一次比较大的振幅的压力波传到感应器，然后没过多久就地震了。

再继续放大，我们会发现时间信号其实大部分都是噪音，肉眼很难看到任何趋向的可寻之迹。而TTF则是接近于一个离散的阶梯函数。根据主办方的说法，时间数据的采样率是400万赫兹，也就是每$2.5\times 10^{-7}$秒有一次时间数据采集。

测试集则和训练集不太一样，是2600多个分段的时间信号（无TTF），每段15万个时间点（$0.375$秒），我们随机选取几段如下：

这次比赛的最终目的就是：

对测试集中的每一段时间信号，预测出该信号结束时的time to failure值，即是预测这些波动结束后还有多久地震。评分的标准是MAE，即是预测值和实际值的绝对值差的平均值。

参赛者想取得好成绩的话，最需要做两项工作：

1. 从这些看似杂乱无章的时间数据中的找到有“规律”的特征（feature engineering），来指引机器学习的模型预测TTF。
2. 建立一个可信可靠robust的cross-validation（以下简称CV），从而让我们敢相信一个模型的在训练集中分出来的CV集的表现可以类推至不知道目标的测试集。

听起来似乎挺容易的，但赛后学习金牌的方法，方知金牌的答案对于细节的注重远胜于其他人。引用UC Irvine数学系编程一霸陈龙老师喜欢说的一句话："The devil is in the details."

UCI Math小队参赛

本次比赛的UCI数学系小分队里面，自己算个山寨指导老师，主力队员是Ziteng，他也是UCI数学系机器学习课程的第一期学员。

Ziteng在2019冬季的Math 10课程的期末project中和队友学习了如何用Keras调用很多大型的网络，最终的得分甚至超过了日本假名数据集Kuzushiji-MNIST的作者给出的最高得分benchmark ResNet+VGG的ensemble。

另：如果读者感兴趣，自己可以写写如何利用Kaggle强大的平台，用目标已知数据集搭建自己机器学习或者科学计算的in-class competition。

UCI Math小分队最终在private leaderboard上取得铜牌的模型（public 311名，private 384名），便是主要由Ziteng筛选的特征，基于6-fold的CV使用XGboost和LightGBM训练出来的。

我们生成了很多特征，滤掉高频噪音之后，频率域下的积分，移动最大值，调和平均；时间域的移动平均，ARMA fitting，等等等等等。感觉很多是都是无效的噪音，虽然和目标TTF相关性很高，但频率域下甚至是时间域上的基于积分的特征相互之间的相关性也很高，要知道，特征之间的相关性高并不是好事。训练数据中出现了几次“雷声大，雨点小”的地震，空有超大的时间信号的振幅，而实际地震则是没有。如果简单用这些特征来建立模型，模型的训练很有可能会被这种雷大没雨的带偏，增加了筛选特征的难度。

一些特征和目标之间的对比可见下图：

除了gradient boosting模型，我们还试验了简单的feedforward的神经网络和recurrent神经网络（RNN）的模型，还试验了在loss函数中加入各种自定义的regularization，但效果并不理想。

另：虽然RNN模型并不理想，通过这次比赛，自己也学习了通过Keras前端调用tensorflow的时候可以使用一个叫fit_generator的对付超大型数据训练的犀利武器。由于数据本身已经太大，生成特征之后一次性全部读入内存根本不现实，于是应用数学工作者们发明了每次只读入少量的样本特征进行训练（mini-batch）。到了真正实现的时候，则是用过Python的generator来完成的。每个generator都含有Python特有的While True... yield桥段。

比如如下是自己基于starter code实现的generator：下面这个generator就可以把6亿个时间点的数据分成带重叠的2万小段，每次随机挑选32段15万个时间点数据，喂给create_feature函数，create_feature是任何可以由时间域数据生成特征的函数就可以了。

def generator(data, min_index=0, 
              max_index=None, 
              batch_size=32, 
              n_steps=150, 
              step_length=1000):
 
    if max_index is None:
        max_index = len(data) - 1
     
    while True:
        # Pick indices of ending positions
        seg_length = n_steps * step_length
        index_range = range(min_index + seg_length, max_index, 20000)
        rows = np.random.choice(index_range, size=batch_size, replace=False)
         
        # Initialize feature matrices and targets
        samples = np.zeros((batch_size, n_steps, n_features))
        
        ## adding a sequence of targets or a single target
        targets = np.zeros((batch_size, n_targets))
        
        for j, row in enumerate(rows):
            samples[j] = create_feature(data[:, 0], 
                                  last_index = row, 
                                  n_steps = n_steps,
                                  step_length = step_length)
            
            if n_targets == 1:
                ## single target
                targets[j] = data[row - 1, 1]
            elif n_targets == 2:
                ## here the training needs to be chosen as the ones after the first earthquake
                targets[j,0] = data[row - 1, 1]
                ## time_after_failure 
                taf_idx = index_earthquake_start[np.sum(row > index_earthquake_start) - 1]
                targets[j,1] = eps*(data[taf_idx, 1] - targets[j,0])
            elif n_targets > 2:
                ## multiple targets (preferably odd number)
                for i in range(n_targets):
                    # targets are all in one segments
                    # targets are spanning multiple neighboring segments
                    targets[j,i] = data[row - 1 - i*(seg_length//2), 1] 

        yield samples, targets

后记

本次比赛其实经历了一次重大事故，用Kaggle论坛的上的话叫做data leakage。细心的网友发现，其实本次比赛的数据是来源于Penn State岩石力学实验室的“p4677”号地震实验。虽然p4677在Penn State的网站上已经没有，但有类似的数据比如p4581，热心网友还把训练集的TTF和p4677试验在前面提到过的论文中振幅log尺度下的图进行了比较，如下图。甚至还有网友用了直尺去测量论文中的图片来验证，几位Kaggle宗师戏称这次事故为“ruler gate”。

根据组织者的描述，训练集和测试集取材于一次连续试验，那么测试集明显就是紧接着训练集之后的时间数据，切成15万的小段，随机打乱之后成为比赛的测试集数据。

那次实验在他们的论文中已经提到（如下图）。细心网友发现，通过leaderboard probing，提交一份全部为0的答案可以得知public测试集的TTF平均为4.017（小型地震），推而得知private测试集的平均应该是6.32（中大型地震）。

基于这个信息，虽然有违机器学习的最初的目标，往模型预测里面加上一个简单的后处理（post-procesing）可有奇效：把模型的预测值人工调整到6.32（全部预测值除以现在的平均然后在乘以6.32即可）。赛后开榜发现，很多金牌得主都用了这个技巧。

不过赛后学习第一名的答案还是让自己收获良多，他们根据这次leakage，建立了自己特别的robust cross-validation策略和训练集甄选方式，让训练出的模型直接就针对预测中大型地震的预测，而不是人工干预最终的预测结果。相比我们的模型用的只是普通的$k$-fold CV，用“out-of-fold”预测的平均值作为该模型CV值。

第一名的队长选出的两个在private leaderboard上排第一的答案，在public leaderboard上只能排2000多名（前450名可以拿牌）。赛后开榜自己做了个简单测试，把Ziteng训练出来private榜得分最高的模型，后处理到6.32平均TTF，居然能排在private第23名（？！）离金牌一名之遥。但是，这个答案在public上连1500名都不到，作为山寨领队这个答案肯定是不敢选为最终答案的。所以，经历了提交答案截止前五分钟的纠结，自己明白要第一名做出这个决定需要多么大的心脏，用一位宗师的话来说就是：

如果你建立了robust cross-validation，你要相信数学的力量。

有心网友甚至通过Kaggle的API制作了这次地震预测比赛的“shake-up”统计：4500多个队伍，public榜和private榜的平均名次差别是872，最大名次上升是4024（？！），最大名次下降是3674，也就是说，前一分钟可能还是leaderboard上的银牌，而开榜之后就变成了“幺鸭子”。而下图这位仁兄，最后一次提交答案是开榜之前两个月，估计觉得CV和leaderboard太没规律了，没意思不玩了，结果得了金牌。

续

本次比赛中小露身手的UCI数学系本科生Ziteng，在2019暑假参加Chemistry and Mathematics in Phase Space主办的量子化学比赛中独自拿到了铜牌（2749队中排名221）。

这次比赛的数据是13万个分子中700万个原子磁场上偶联（j-coupling）的强度，是由Density Function Theory中的Kohn-Sham模型计算得出的。如果有时间笔者也打算介绍下这次比赛，自己学习了如何在tensorflow中用Keras前端或者chainer前端实现自定义的图卷积神经网络（graph convolutional neural network），在网络里面加入基于graph Laplacian的卷积层生成特征给regressor。其间还学习使用各式各样的量子化学Python强大的开源包，修改底层code加上并行化来即时生成想要的特征，助力让自己得到了第一个Kaggle上的银牌（39名）。