一、简介

比赛地址：传送门

描述：该比赛目标是使用地震信号来预测地震发生的时间。

训练数据：两列，第一列是地震的信号数值acoustic_data，第二列是还有多久发生地震time_to_failure（记为TTF）。数据量是629145480行，即629145480个acoustic_data和TTF。

测试数据：有2624个测试文件，每个文件有150,000行数据，只有acoustic_data，没有TTF，目标是预测这15w个acoustic_data后还有多久发生地震，即最后数据对应的TTF。

训练数据和测试数据来源相同，只是训练数据是整个数据集的前60%，测试数据是后40%。所以测试数据是将后40%按15w划分，一共生成了2624条数据。

二、分析数据特点

1、首先对training set中的acoustic_data和TTF整体分布画图分析：

图片来自kernel：传送门

可以看出，训练数据中一共有15个完整的地震，开始和结尾还有两个不完整部分。每个地震的数据长度不同，最长的是16s+的数据，最短的不到8s。

2、然后分析一个地震的波形图，对第一个地震画图：

（图片还是来自上面的kernel）

结合两个图，可以看到在每次地震发生前的某个时刻，地震波会有个很大的波动。

3、数据处理

训练数据629145480 / 1500000 = 4194.3032，即以15w为块划分的话可以得到4194条训练数据。

4、数据评测：

public榜是13%的数据，private榜是83%的数据。public榜数据太少，基本没法反应算法的好坏，所以还是要相信CV的。（这场比赛最后结果出来后，波动非常大，一共4500+队伍，很多人波动都在1000,2000以上。我认为主要是有两个原因：一时public数据太少，很多人对public榜过拟合了；二是训练数据和测试数据分布不同，这个原因应该是主要的吧，后面会介绍大佬们是如何处理数据分布的）

三、我的方案

（可以跳过，直接看下面大佬们的各种操作）

（一）第一阶段，使用time series model（rnn+cnn）

1、特征提取

训练数据以15w划分的话可以得到4194条数据，但是这样数据量比较小。所以通过overlap来扩充数据，比如每间隔7.5w条数据采样15w条数据，那么数据量就会翻倍。为防止过拟合，也可以随机overlap。

对于采样的一条数据即15w个信号，进行数据量缩减提特征，比如设置窗口大小为500，求窗口内数据的mean, std, max, min等，那么每个特征的长度就是15w/500=300，下面是我使用的特征（data即一条150000的数据）：

def extract_feature(data):_x = []_x.append(np.std(data))  # std_x.append(np.max(data))  # max_x.append(np.min(data))  # min_x.append((np.mean(data)+np.median(data))/2)  # (mean+median)/2_x.append(np.max(np.abs(data)))  # max_abs_x.append(np.mean(np.abs(data)))  # abs mean_x.append(np.percentile(np.abs(data), 90))  # abs 90%_x.append(np.mean(np.abs(np.diff(data))))  # diff_abs_mean_x.append(crests(data))  # crest_x.append(stats.skew(data))  # skew_x.append(stats.kurtosis(data))  # kurtosis_x.append(np.percentile(data, 75) - np.percentile(data, 25))  # 75% - 25%_x.append(classic_sta_lta(data, 100, 500).mean())  # sta lta_x.append(add_trend_feature(data, abs_values=True))  # treadreturn _x

其中里面用的几个函数（有些可以直接调库，不必自己写）：

def classic_sta_lta(x, length_sta, length_lta):sta = np.cumsum(x ** 2)# Convert to floatsta = np.require(sta, dtype=np.float)# Copy for LTAlta = sta.copy()# Compute the STA and the LTAsta[length_sta:] = sta[length_sta:] - sta[:-length_sta]sta /= length_stalta[length_lta:] = lta[length_lta:] - lta[:-length_lta]lta /= length_lta# Pad zerossta[:length_lta - 1] = 0# Avoid division by zero by setting zero values to tiny floatdtiny = np.finfo(0.0).tinyidx = lta < dtinylta[idx] = dtinyreturn sta / ltadef crests(x):doublediff = np.diff(np.sign(np.diff(x)))_crests = np.sum(doublediff == -2)return _crestsdef add_trend_feature(arr, abs_values=False):idx = np.array(range(len(arr)))if abs_values:arr = np.abs(arr)lr = LinearRegression()lr.fit(idx.reshape(-1, 1), arr)return lr.coef_[0]

这一阶段的特征比较少，也比较简单，也并不是所有特征都用了，有的提交只是用了里面的一部分。

2、构建模型

模型使用CNN+LSTM，输入数据是8387x14x300，即一共8000+条数据（用overlap扩充），14个特征序列，每个序列的长度为300。模型使用两层1D卷积和两层LSTM，网络结构如下（可能不是最优结构，仅供参考）：

def build_net():l = 150000 // windowprint(l)model = Sequential()model.add(Conv1D(64, 8, strides=8, input_shape=(l+1, n_feature)))model.add(Conv1D(64, 2, strides=2))model.add(BatchNormalization())model.add(Activation('relu'))model.add(LSTM(80, return_sequences=True))model.add(LSTM(80))model.add(Dense(1))print(model.summary())adam = keras.optimizers.Adam(lr=0.005)  # , beta_1=0.9, beta_2=0.999, epsilon=1e-08)model.compile(optimizer=adam, loss='mae')return model

这一阶段是刚开始自己摸索的阶段，特征和模型都比较弱，没有做CV，提交的public榜mae损失也比较大。但其实private榜揭开之后发现，最好评分的提交还是来自这一阶段（能在private榜100以内），但由于太依赖public榜而没有选择该方案作为最终结果。分析原因的话：我认为是该阶段特征和模型相对于第二阶段来说比较弱，泛化能力强一些，还没有对public榜过拟合。

（二）第二阶段，使用feature model（各种Gradient Boosting，RF，DNN，SVR）

这一阶段的主要思路来自kernel中各种方案，包括特征提取、特征选择，模型构建，CV验证等。kernel中比较常见的思路是用Python自带的一些工具包生成大量的特征，然后根据特征重要度选择一小部分特征，去掉作用不大的冗余特征。模型一般使用boosting包括lightgbm, xgboost, catboost, GBDT，其他像RandomForest，DNN，SVR也有，不过从public榜来看，没有boosting算法得分高。

1、首先，特征提取

主要来源：特征来源1，特征来源2，特征来源3

其中仅来源3就有1500+个特征，这几个来源中有很多特征是相同和相似的，而且很多特征是没有用的，需要筛选一下。

2、特征选择

这里我使用的选择方法有两种，第一是boosting算法自带的特征重要度选择，第二是pearson相关系数分析。其中来源2和来源3用boosting进行特征选择的图如下：

来源2：

来源3：

对于来源1中的特征，我使用lgbm算法对特征重要度排序，然后选择前30个左右重要特征，然后从来源2和来源3的特征中也选择重要度高的一部分（至于数目多少，我是不断尝试增加的），最终特征在100个左右。

其中用lgbm分析特征重要性：传送门（待补充）

使用pearson系数分析相关性代码：

from scipy.stats import pearsonr
def calc_pear(X_train, Y_train):pear = []for i in range(X_train.shape[1]):pear.append(pearsonr(X_train[:, i], Y_train)[0])return np.nan_to_num(np.array(pear))

3、可视化

画图观察某个特征和TTF的变化趋势，讨论区里大佬们的特征：传送门，下图是我的一个特征图：

代码：

def plot2():training_data = pd.read_csv('../input/train_features_8k_100+.csv')X = training_data.drop(['target', 'seg_id'], axis=1)y = training_data.target.valuesprint(X.shape)name = 'cid_ce_1'x = X[name].valuesx_axis = np.arange(0, len(y))fig, ax1 = plt.subplots(figsize=(10, 5))ax2 = ax1.twinx()ax1.plot(x_axis, x, c='green')ax2.plot(x_axis, y, c='red')plt.title(name)plt.show()

因为我的数据是overlap扩充的，所以有8k多条。通过观察图，可以发现特征和TTF之间有一定的正相关性。另外，通过图还可以看出很多特征是非常相似的，波形图基本是完全一样的，多个这样的特征可能会造成重复，所以还需要把共线的特征去掉。

4、构建模型，CV验证

这里以lgbm为例，代码框架来自kernel。

参数配置：

def build_model():params = {# 'task': 'train','boosting_type': 'gbdt',  # default: gbdt'objective': 'gamma','metric': 'mae',  # huber gamma'reg_alpha': 0.1302650970728192,  # def = 0, L1'reg_lambda': 0.3603427518866501,  # def = 1, L2# 'lambda_l1': 0.1,  # default: 0'num_leaves': 32,  # default: 31'min_data_in_leaf': 15,  # default: 20'max_depth': -1,  # default: -1'learning_rate': 0.005,'feature_fraction': 0.8,'bagging_fraction': 0.8,'bagging_freq': 2,'bagging_seed': seed,'verbose': 0,'num_threads': 8,'random_state': 42}model = lgb.LGBMRegressor(**params, n_estimators=10000)return model

主要运行框架，其中使用5_fold CV验证，数据是已经生成好保存在本地的。

def train_and_predict():training_data = pd.read_csv('../input/train_features_100.csv')training_data = training_data.sample(frac=1)X = training_data.drop(['target', 'seg_id'], axis=1)y = training_data.targetfolds = KFold(n_splits=5, shuffle=True, random_state=seed)oof_preds = np.zeros((len(X), 1))test_data = pd.read_csv('../input/test_features_100.csv')X_test = test_data.drop(['target', 'seg_id'], axis=1)test_preds = np.zeros((len(X_test), 1))model = Nonefor _fold, (_train, _val) in enumerate(folds.split(X)):print('Current fold:', _fold)x_train, y_train = X.iloc[_train], y.iloc[_train]x_val, y_val = X.iloc[_val], y.iloc[_val]model = build_model()model.fit(x_train, y_train,eval_set=[(x_train, y_train), (x_val, y_val)],eval_metric='mae',verbose=1,early_stopping_rounds=500)val_pred = model.predict(x_val, num_iteration=model.best_iteration_)test_fold_pred = model.predict(X_test, num_iteration=model.best_iteration_)print("Val MAE = {}".format(metrics.mean_absolute_error(y_val, val_pred)))oof_preds[_val, :] = val_pred.reshape((-1, 1))test_preds += test_fold_pred.reshape((-1, 1))test_preds /= 5oof_score = metrics.mean_absolute_error(y, oof_preds)print("Val Mean MAE = {}".format(oof_score))return test_preds, model

lgbm是我在public榜上效果最佳的算法（也可能是对public榜过拟合最严重的一个算法），其他方法还有XGB、GBDT、RF、DNN、GP、KNN等（相关链接之后慢慢补充）

GP方法：符号回归，使用遗传算法实现，即对所有特征进行各种组合，选出拟合最好的表达式。有大佬在kernel里使用这种方法，在LB上简直是奇效，比lgbm效果还要好，kernel地址：地址1，地址2。其关于gp方法的kernel是直接调用gplearn库函数：地址1，地址2，但是效果没有前面的好，我自己跑的时候效果更差，CV误差很大，提交结果也不好，可能是没调好GP超参数（猜测大佬的GP可能是自己实现的，加了一些优化）。

DNN方法：地址

四、冠军及金牌方案

总结发现，金牌区的方案基本都会对直接或间接对数据分布做处理。而且这个比赛数据有leak，使用的数据集是主办方在p4677论文中的数据集，虽然无法知道每个测试数据的结果，但可以通过计算得到均值。（我试了一下，将之前提交的结果均值平衡到6.31后，会有很大的提升，而且基本都在private的金牌区）

（以下总结中个人观点可能有误，请大佬及时指出）

1、冠军方案

传送门

其中有两个idea我感觉很重要：

（1）训练数据集选择：由上面分析可知一共有17个地震（加上前后两个不完整的），从17个地震中采样10个（进行10000次），和test比较KS statistic（检验是不是同分布），选择分布最相近的10个地震作为训练数据。

（2）特征选择：使用了4个不想管的特征，而且选择特征的规则是检测该特征在train和test上的ks statistic，如果p-value > 0，即同分布，则接受该特征。

冠军队的另外一种方案：传送门，private得分更好一些

2、另一种超越冠军的方案

传送门

比赛完后，giba大佬在kernel里放出的方案。

使用了全部的训练数据；使用wavlet去噪；特征数目为673个；去除基于mean/min/max/std等的特征（为什么？）；根据p4677将测试结果的mean均衡到6.31（这个提分很大）

3、public 2nd / private 7th 方案

传送门（大佬就是大佬，不管怎么波动都稳居top）

（1）16_fold CV

（2）每5w数据作为起点采样训练数据，即3倍数据

（3）outliers: 由于15w数据是按顺序截取的，所以在15w中间可能会出现发生地震时间重置的情况。所以大佬训练了一个binary模型来识别是否出现这种情况，然后针对这种情况再训练额外的模型

（4）特征12个，具体见kernel链接

（5）模型使用了lgb, pygam, knn，其中gam > lgb> knn，然后stacking，lgb was used for the sencond level model.

（6）根据分布选择训练数据：

根据p4677得出test的均值，From it I estimated its mean to be 6.35 which is quite close to the actual 6.32. From this we can estimate density function of TTF for train and test. We then used sample weights to map the train distribution to the test distribution.

（7）对test结果再做处理

set all high variance segments TTF to 0.31; combine binary models to fix segments TTF at 0.2.

还有一点是，每条训练和测试数据都减去本条数据的均值。

4、其他方案

第二名根据p4677只使用部分训练数据；第三名也是将均值均衡到6.31；第五名也对predict均值做了处理。

第三名介绍了如何根据p4677得出测试数据的均值：传送门

第六名kernel介绍了如何得出p4677的数据就是本比赛的数据：传送门

五、经验/教训总结

第一次参加比赛，完全小白，一路基本是跟着kernel和discussion中大佬们的思路，学了不少东西。之前只是对boosting算法原理了解一点，通过比赛对这些算法更了解一些了吧。其中还学了一些基本的CV验证，特征工程包括generate feature、feature selection等，模型ensembling，收获还是挺大的。但是结果还是挺遗憾的，比赛结束前public榜最好到30+的名次，然后一路掉到100多名，然而private榜揭开后血崩！！一下掉了700多名次。也收获了不少教训：

1、相信CV，别太关注LB，尤其像这个比赛数据量分布很不平衡的情况(13 : 87)

2、对问题和数据要深入的研究，虽然这个比赛有数据的leak，但只有下功夫去研究才会发现

LANL Earthquake Prediction收获相关推荐

kaggle-地震预测-LANL Earthquake Prediction
kaggle-地震预测-LANL Earthquake Prediction比赛记录 1.baseline 加入这次比赛的时间比较晚,所以直接先follow了别人的一个baseline来start找找 ...
label y 训练集测试集x_Adversarial validation-对抗验证| 一种解决训练集与测试集分布不一致的方法...
导语: 马上就要五一了,祝全世界人民五一快乐!在这之前,想过好几个准备这些天可以完成的专题,比如Boosting系列在搞点最近几年的新玩意,或者开一个新专题,如心心念念的GNN/GCN(主要是又可以去 ...
机器学习里的bias的一些处理方法和思考
点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达作者:砍手豪链接:https://zhuanlan.zhihu. ...
【深度学习】【物联网】深度解读：深度学习在IoT大数据和流分析中的应用
作者|Natalie 编辑|Emily AI 前线导读:在物联网时代,大量的感知器每天都在收集并产生着涉及各个领域的数据.由于商业和生活质量提升方面的诉求,应用物联网(IoT)技术对大数据流进行分析是 ...
时间序列分析（TSA）
文章目录问题分类应用领域数据集博文论文待阅问题分类根据数据格式离散值序列连续值序列根据采样频率均匀采样.等间隔采样不均匀采样根据变量数目单变量时间序列 UTS 多变量时间 ...
德国公民信用相关数据_作为公民数据科学家，没有任何事
德国公民信用相关数据数据科学,意见(Data Science, Opinion) Dear Aspiring Data Scientist, 亲爱的有抱负的数据科学家, Before you sta ...
【深度学习】深度解读：深度学习在IoT大数据和流分析中的应用
来源:网络大数据(ID:raincent_com) 摘要:这篇论文对于使用深度学习来改进IoT领域的数据分析和学习方法进行了详细的综述. 在物联网时代,大量的感知器每天都在收集并产生着涉及各个领域的数 ...
make用法总结（基于材料：“老托福听力93篇”）
在"老托福听力93篇"材料中,make的几种用法总结如下: (1)vt. 制造,制作: (2)vt. 建造: (3)make + n.1 + n.2, 使n.1成为n.2; (4) ...
首次kaggle竞赛（Riiid! Answer Correctness Prediction），收获1银1铜，分享自己的baseline：Riiid基于Catboost的准确率预测。
文章目录一.defaultdict 1.认识defaultdict: 2.如何使用defaultdict: 二.datatable包 1.datatable使用 2.数据转换三.特征构造 1.函数 ...

LANL Earthquake Prediction收获