Kaggle比赛系列：（1）LANL Earthquake Prediction

比赛介绍

地震预测主要是为了解决3个问题：什么时候会发生？（when）会在哪里发生？（where）会发生多大规模的地震？(what)。这个竞赛主要是预测地震什么时候会发生，具体来说，将根据实时地震数据预测实验室地震发生前的剩余时间（地震是模拟的，数据是真实的）
数据集：链接: https://pan.baidu.com/s/1oOZK6jc8BzxKUcp8OsgHMQ 提取码: 9i8c
数据集描述：
(1) train.csv (9G):训练数据集，只包含信号和地震发生剩余时间两个字段，无空值；
(2) test folder :测试数据集，共包含2624个文件，对应2624次预测结果；
(3) sample_submission.csv:提交模板，每个预测字段对应一个预测时间。

参赛大神的比赛记录

原文地址：Earthquakes FE. More features and samples（Andrew Lukyanenko）

1.导入数据处理和模型包

import numpy as np
import pandas as pd
import os

import matplotlib.pyplot as plt
%matplotlib inline
from tqdm import tqdm_notebook
from sklearn.preprocessing import StandardScaler
from sklearn.svm import NuSVR, SVR#NuSVR, SVR是SVM的变种方法
from sklearn.metrics import mean_absolute_error
pd.options.display.precision = 15

import lightgbm as lgb#新安装
import xgboost as xgb#新安装
import time
import datetime
from catboost import CatBoostRegressor#新安装
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold
from sklearn.metrics import mean_absolute_error
from sklearn.linear_model import LinearRegression
import gc
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')#注意这里不能用双引号

from scipy.signal import hilbert
from scipy.signal import hann
from scipy.signal import convolve
from scipy import stats
from sklearn.kernel_ridge import KernelRidge

在Anaconda环境下运行，基本上只有lightgbm、xgboost、catboost三个包需要安装
安装方法

pip install lightgbm xgboost catboost

2.数据加载

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

CPU times: user 2min 7s, sys: 15.7 s, total: 2min 22s
Wall time: 2min 20s
%%time 将会给出cell的代码运行一次所花费的时间

3.观察训练数据集

train_acoustic_data_small = train['acoustic_data'].values[::50]#切片抽样，抽样间隔为50
train_time_to_failure_small = train['time_to_failure'].values[::50]

fig, ax1 = plt.subplots(figsize=(16, 8))
plt.title("Trends of acoustic_data and time_to_failure. 2% of data (sampled)")
plt.plot(train_acoustic_data_small, color='b')
ax1.set_ylabel('acoustic_data', color='b')
plt.legend(['acoustic_data'])
ax2 = ax1.twinx()#twinx（）函数表示共用横坐标
plt.plot(train_time_to_failure_small, color='g')
ax2.set_ylabel('time_to_failure', color='g')
plt.legend(['time_to_failure'], loc=(0.875, 0.9))#不设置loc参数，图注会被覆盖
plt.grid(False)

del train_acoustic_data_small
del train_time_to_failure_small

当在一个figure上面出现两个纵坐标的时候，注意绘图操作都采用plt,只有在设置纵坐标时才分别使用ax1和ax0。
从上图中观察可以发现：
（1）地震发生前的信号会出现巨大波动，且自然数据呈现出周期性；
（2）从可视化结果可以看出，当信号出现巨大波动，并伴随着一系列微动信号的时候可以预测地震发生时间。
（3）作者尝试在一定区间内选择一个阈值（1000或2000）与信号的最大值作比较，以此作为一个预测特征，但是发现没有作用。

4. 特征工程

# Create a training file with simple derived features
rows = 150_000#与150000相同，每次采样数量
segments = int(np.floor(train.shape[0] / rows))# 下取整.抽样频率,总数是629145480,超过4000次抽样
#特征1:趋势特征
def add_trend_feature(arr, abs_values=False):#添加信号的趋势特征(自然趋势,绝对值趋势),特征6
    idx = np.array(range(len(arr)))# 将Series类型转化为ndarray类型,现在idx是从零开始,arr是从1开始
    if abs_values:
        arr = np.abs(arr)#当abs_values=True时,取绝对值
    lr = LinearRegression()
    lr.fit(idx.reshape(-1, 1), arr)
    return lr.coef_[0]#逻辑回归的特征就是Θ的取值:lr.coef_指的是各个特征对结果影响的大小(特征系数)θ>0, 正相关；θ<0，负相关；θ值绝对值越大，特征值对标记影响越大；
#特征2:震相特征
def classic_sta_lta(x, length_sta, length_lta):#震相自动识别算法STA/LTA
    sta = np.cumsum(x ** 2)#cumsum函数是累加函数axis=0,按行累加;axis=1,按列累加
    # Convert to float
    sta = np.require(sta, dtype=np.float)
    # Copy for LTA
    lta = sta.copy()
    # Compute the STA and the LTA
    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]
    sta /= length_sta
    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]
    lta /= length_lta

    # Pad zeros
    sta[:length_lta - 1] = 0

    # Avoid division by zero by setting zero values to tiny float
    dtiny = np.finfo(0.0).tiny
    idx = lta < dtiny
    lta[idx] = dtiny

    return sta / lta
#制作批次训练数据集
X_tr = pd.DataFrame(index=range(segments), dtype=np.float64)#这里的segment相当于深度学习经常使用的batch

y_tr = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])
#特征3:聚合特征
total_mean = train['acoustic_data'].mean()#平均值
total_std = train['acoustic_data'].std()#标准差
total_max = train['acoustic_data'].max()#最大值
total_min = train['acoustic_data'].min()#最小值
total_sum = train['acoustic_data'].sum()#累加和
total_abs_sum = np.abs(train['acoustic_data']).sum()#绝对值累加和
#特征4:差值比率平均值(类似比率这种数值型特征最好使用平均值)
def calc_change_rate(x):
    change = (np.diff(x) / x[:-1]).values#np.diff是累减,这里得到了差值的比率
    change = change[np.nonzero(change)[0]]#np.nonzero函数是numpy中用于得到数组array中非零元素的位置（数组索引）的函数
    change = change[~np.isnan(change)]#排除空值	
    change = change[change != -np.inf]#排除负无穷
    change = change[change != np.inf]#排除正无穷
    return np.mean(change)#获得比率的平均值

for segment in tqdm_notebook(range(segments)):#tqdm_notebook模块是用来显示进度条的
    seg = train.iloc[segment*rows:segment*rows+rows]#iloc是根据索引定位记录的,loc则可以使用column名和index名进行定位
    x = pd.Series(seg['acoustic_data'].values)#将DataFrame中的一列(signal)转化为Series
    y = seg['time_to_failure'].values[-1]#这里取的是quaketime的最后一个值,以这个值作为预测对像
    #添加特征
    y_tr.loc[segment, 'time_to_failure'] = y
    X_tr.loc[segment, 'mean'] = x.mean()
    X_tr.loc[segment, 'std'] = x.std()
    X_tr.loc[segment, 'max'] = x.max()
    X_tr.loc[segment, 'min'] = x.min()
    
    X_tr.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x))#这里并没有采用绝对值,只是使用了临值差
    X_tr.loc[segment, 'mean_change_rate'] = calc_change_rate(x)#对整体segment应用
    X_tr.loc[segment, 'abs_max'] = np.abs(x).max()
    X_tr.loc[segment, 'abs_min'] = np.abs(x).min()
    
    X_tr.loc[segment, 'std_first_50000'] = x[:50000].std()#前50000条数据的标准差
    X_tr.loc[segment, 'std_last_50000'] = x[-50000:].std()#后50000条数据的标准差
    X_tr.loc[segment, 'std_first_10000'] = x[:10000].std()#前10000条数据的标准差
    X_tr.loc[segment, 'std_last_10000'] = x[-10000:].std()#后10000条数据的标准差
    
    X_tr.loc[segment, 'avg_first_50000'] = x[:50000].mean()
    X_tr.loc[segment, 'avg_last_50000'] = x[-50000:].mean()
    X_tr.loc[segment, 'avg_first_10000'] = x[:10000].mean()
    X_tr.loc[segment, 'avg_last_10000'] = x[-10000:].mean()
    
    X_tr.loc[segment, 'min_first_50000'] = x[:50000].min()
    X_tr.loc[segment, 'min_last_50000'] = x[-50000:].min()
    X_tr.loc[segment, 'min_first_10000'] = x[:10000].min()
    X_tr.loc[segment, 'min_last_10000'] = x[-10000:].min()
    
    X_tr.loc[segment, 'max_first_50000'] = x[:50000].max()
    X_tr.loc[segment, 'max_last_50000'] = x[-50000:].max()
    X_tr.loc[segment, 'max_first_10000'] = x[:10000].max()
    X_tr.loc[segment, 'max_last_10000'] = x[-10000:].max()
    
    X_tr.loc[segment, 'max_to_min'] = x.max() / np.abs(x.min())#比值
    X_tr.loc[segment, 'max_to_min_diff'] = x.max() - np.abs(x.min())#差值
    X_tr.loc[segment, 'count_big'] = len(x[np.abs(x) > 500])#超过500阈值的值
    X_tr.loc[segment, 'sum'] = x.sum()
    
    X_tr.loc[segment, 'mean_change_rate_first_50000'] = calc_change_rate(x[:50000])#对特殊区域segment应用
    X_tr.loc[segment, 'mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])
    X_tr.loc[segment, 'mean_change_rate_first_10000'] = calc_change_rate(x[:10000])
    X_tr.loc[segment, 'mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])
    
    X_tr.loc[segment, 'q95'] = np.quantile(x, 0.95)#
    X_tr.loc[segment, 'q99'] = np.quantile(x, 0.99)
    X_tr.loc[segment, 'q05'] = np.quantile(