{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 概要\n\n* このコンペで行った内容を共有します．\n* 提出したモデルには複数のKernelを使用しており，コードももう少し複雑です．\n* Kaggleの紹介と合わせて短時間で共有するため，完結に実行できるレベルにコードを略してまとめています．"},{"metadata":{},"cell_type":"markdown","source":"# インポートとファイル確認\n\n* test\n  * zipファイルで，拡張機能によって中のファイルを直接読み込めます．\n* train.csv\n  * 学習用に容易された巨大な時系列ファイルです．長さは6億行余り．\n* sample_submission.csv\n  * 提出用フォーマットのCSVファイルです．"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nimport os, time, gc\nfrom tqdm import tqdm_notebook\n\nimport librosa\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_absolute_error\nfrom catboost import CatBoostRegressor\n\nos.listdir('../input')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 読込\n\n* pandasという分析ライブラリを使用して学習データを読み込みます\n* 空きメモリに余裕が無くなるため，実際はpandasの機能で1500万ずつ分割して読み込んでいます．\n* 列の説明\n  * acoustic_data -> 音響列．16bitで読まないとメモリエラー\n  * time_to_feailure -> 地震までの時間がsec単位で入っています\n* 読み込むのに約2分半かかります"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n\ntrain = pd.read_csv(\n    '../input/train.csv', \n    dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\nprint('loaded:', train.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 波形を確認\n\n* サンプリングしないとグラフでメモリオーバー．\n* オレンジ波形の大きい揺れが地震で，おおよそ地震のタイミングで青い線が0になっています．\n* たった16回の地震から次の地震を予測するのは難しそうです，，"},{"metadata":{"trusted":true},"cell_type":"code","source":"data = train['acoustic_data'].values[::200]\nttf = train['time_to_failure'].values[::200]\n\nfig, ax = plt.subplots(1, 1, figsize=(16, 8))\nax.plot(data, label='acoustic_data')\nax.plot(ttf * 100, label='time_to_feailure')\nplt.legend()\nplt.show()\n\nfig, ax = plt.subplots(1, 1, figsize=(16, 8))\nax.plot(data[:200000], label='acoustic_data')\nax.plot(ttf[:200000] * 100, label='time_to_feailure')\nplt.legend()\nplt.show()\n\ndel data, ttf\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 特徴作成\n\n* 他の方のKernelsを参考しつつ，150,000ごとの波形データから特徴を作成していきました．\n* 9つの特徴だけに省略していますが，実際は沢山の特徴を作っています．\n* 実際の特徴量は下の「算出した特徴量」を参照してください．\n* ↓に挙げた以外にも色々な組み合わせを試して，数千個の中から効果がありそうなものを組み合わせて700個程度に絞りました．\n* 大体ライブラリを使って容易に計算できますが，そこから過学習しない良い特徴量を選ぶのが大変でした．\n* 各特徴量の詳細を学んでいる時間がなく，ライブラリをただ使っただけの特徴も多いです．\n* データ追加とカラム追加を同じメソッドにすることで，特徴を増減させた時のコード修正量を減らしています．\n* パフォーマンスがあまり変わらなかったので，特徴量をまとめるところはnumpy配列ではなくpythonのlistを使用しています．\n\n----\n\n### 算出した特徴量\n\n* 基本統計量\n  * 平均，標準偏差，尖度，幅，分位数，四分位範囲，ピーク数\n  * 単純移動平均，単純移動標準偏差，指数移動平均，指数移動標準偏差（各統計量を算出）\n* 周波数変換\n  * フーリエ変換の実数部・虚数部（特徴がある部分以外をカットして移動平均，移動偏差×各統計量を算出）\n  * MFCC：メル周波数ケプストラム係数（2次元データの特徴がある部分の平均値と各統計量を算出）\n  * PSD:パワースペクトル密度推定（パラメータ変えて各統計量を算出）\n* その他統計量\n  * STALTA（地震予測によく使用される古典的なアルゴリズム）\n  * 歪度、平均絶対偏差，偏回帰係数，調和平均，幾何平均，基本統計量のパラメータ変更色々\n  * 自己相関，歪度の棄却限界値，モーメント，k統計量とその分散\n  * Time reversal asymmetry statistic（時系列向け特徴抽出アルゴリズムらしい）\n* 更にデータを高周波域を削除して逆変換したものを再度↑の統計量を算出したもの\n\n左5万とか右5万を使った統計量算出は，過学習を起こすだけであまり意味が無さそうだったので採用しませんでした．"},{"metadata":{"trusted":true},"cell_type":"code","source":"def add(row, key, value, mode=False):\n    if mode:\n        row.append(key)\n    else:\n        row.append(value)\n\ndef agg(row, d, prefix='', mode=False):\n    add(row, 'mean{}'.format(prefix), d.mean(), mode)\n    add(row, 'std{}'.format(prefix), d.std(), mode)\n    add(row, 'kurt{}'.format(prefix), d.kurt(), mode)\n    add(row, 'range{}'.format(prefix), d.max() - d.min(), mode)\n\n    add(row, 'q0.01{}'.format(prefix), np.quantile(d, 0.01), mode)\n    add(row, 'q0.05{}'.format(prefix), np.quantile(d, 0.05), mode)    \n    add(row, 'q0.95{}'.format(prefix), np.quantile(d, 0.95), mode)\n    add(row, 'q0.99{}'.format(prefix), np.quantile(d, 0.99), mode)\n    \n    add(row, 'iqr{}'.format(prefix), np.subtract(*np.percentile(d, (75, 25))), mode)\n\ndef wave2row(d, mode=False):\n    row = []\n\n    agg(row, d, mode=mode)\n    \n    return row\n\ndef make_cols():\n    cols = wave2row(pd.Series(list(range(150000))), mode=True)\n    cols.extend(['seg_id', 'time_to_failure'])\n    return cols","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 例1 FFT"},{"metadata":{},"cell_type":"markdown","source":"* 省いたコードの一例としてフーリエ変換後の実数部のグラフを示します．\n* 左右対称なので左半分のみ，更にここから意味のありそうな左から20000個を使っています．"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16, 6))\nplt.plot(np.real(np.fft.fft(train[0:150000]['acoustic_data']))[1:75000])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 例2 MFCC\n\n* メル周波数ケプストラム係数は有用な特徴量になりました．\n* こちらも意味のありそうな領域に絞って使用しています．\n* この画像による深層学習（ディープラーニング）も試しましたが，時間がかかり過ぎるのと調整が難しかったので諦めました．"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10,8))\nsns.heatmap(librosa.feature.mfcc(\n    train[0:150000]['acoustic_data'].values.astype(np.float32), n_mfcc=64)[2:])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 学習データの特徴量作成\n\n* 150,000毎の特徴（説明変数）と地震までの秒数（目的変数）の表示を作成します．\n* 波形情報から扱いやすい行列情報に変換されます．\n* tqdmというライブラリを使うと，イテレーター部分をtqdm_notebookで囲むだけで進捗を表示してくれて便利です．"},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_train(i, f):\n    data = f[i:i + 150000]\n    \n    row = wave2row(data['acoustic_data'])\n    \n    add(row, 'seg_id', str(i))\n    add(row, 'time_to_failure', data[-1:]['time_to_failure'].values[0])\n    \n    return row\n\nindexes = [i for i in range(0, len(train), 150000)]\ntrain_rows = []\nfor i in tqdm_notebook(indexes):\n    train_rows.append(make_train(i, train))\n\ntrain_rows = pd.DataFrame(train_rows, columns=make_cols())\ntrain_rows.to_csv('train.csv', index=False)\nprint('saved train:', train_rows.shape)\n\ntrain_rows.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# テストデータの特徴量作成\n\n* テストデータも同様に特徴量を作成します．\n* こちらは150,000毎にファイルが分かれているので1ファイルずつ処理します．\n* 最初マルチプロセスで処理していたのですが，マルチプロセスに対応していないライブラリがあったので諦めました．"},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_test(seg_id):\n    data = pd.read_csv('../input/test/' + seg_id + '.csv', dtype={'acoustic_data': np.int16})\n    row = wave2row(data['acoustic_data'])\n    add(row, 'seg_id', seg_id)\n    add(row, 'time_to_failure', 0.0)\n    \n    return row\n\nsubmission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\ntest_rows = []\nfor seg_id in tqdm_notebook(submission.index):\n    test_rows.append(make_test(seg_id))\n\ntest_rows = pd.DataFrame(test_rows, columns=make_cols())\ntest_rows.to_csv('test.csv', index=False)\nprint('saved test', test_rows.shape)\n\ntest_rows.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 学習と推論"},{"metadata":{},"cell_type":"markdown","source":"* 学習データと検証データを組み替えて学習する，クロスバリデーションで学習と推論を行っています．\n* 実際はiterationsは2048，他あまりチューニングはしてません\n* また，以下のようなロジックで特徴選択を繰り返しました．\n  * randam stateを4つ分学習を繰り返して平均を取る\n  * 何度か学習を繰り返して，各説明変数の重要度から200程度の特徴に絞る\n  * 特徴を10個程度のブロックに分けて，ブロックごとに処理（8, 16, 2４, ３２みたいな可変長ブロック）\n  * 特徴を1つずつ増やしていって精度が上がったら追加，下がったら追加しない\n  * 特徴を1つずつ減らしていって精度が上がったら削除，下がったら削除しない\n  * ブロック分繰り返す"},{"metadata":{"trusted":true},"cell_type":"code","source":"def learn(params):\n    local_train = train_rows.copy()\n    local_test = test_rows.copy()\n    \n    y = local_train['time_to_failure']\n    x = local_train.drop(['seg_id', 'time_to_failure'], axis=1)\n    \n    cols = params['cols']\n    if len(cols) == 0:\n        cols  = list(x.columns)\n    print('cols:', len(cols))\n    \n    cols.sort()\n\n    x = x[cols]\n    test_x = local_test.drop(['seg_id', 'time_to_failure'], axis=1)[cols]\n    \n    x, y = x.values, y.values\n    \n    fold = KFold(n_splits=5, random_state=params['random_state'], shuffle=True)\n    \n    val_accs = []\n    val_preds = []\n    weights = []\n    submissions = []\n    \n    for i, (train_index, val_index) in enumerate(fold.split(x, y)):\n        train_x, train_y = x[train_index], y[train_index]\n        val_x, val_y = x[val_index], y[val_index]\n        \n        print('  fold:', i + 1, train_x.shape, val_x.shape)\n        \n        model = CatBoostRegressor(\n            iterations=256, learning_rate=0.015, verbose=32, eval_metric='MAE',\n            use_best_model=True, task_type='GPU')\n        model.fit(train_x, train_y, eval_set=(val_x, val_y))\n        \n        val_pred = model.predict(val_x)\n        val_pred = np.where(val_pred < 0, 0, val_pred)\n        val_acc = mean_absolute_error(val_y, val_pred)\n        val_accs.append(val_acc)\n        \n        val_f = pd.DataFrame()\n        val_f['val_index'] = val_index\n        val_f['time_to_failure'] = val_pred\n        val_preds.append(val_f)\n\n        weight = pd.DataFrame()\n        weight['col'] = cols\n        weight['weight'] = model.feature_importances_\n        weights.append(weight)\n\n        test_pred = model.predict(test_x)\n        test_pred = np.where(test_pred < 0, 0, test_pred)\n        submission = pd.DataFrame()\n        submission['seg_id'] = local_test['seg_id']\n        submission['time_to_failure'] = test_pred\n        submissions.append(submission)\n\n    print('  ', ['{0:.4f}'.format(acc) for acc in val_accs], '{0:.4f}'.format(np.mean(val_accs)))\n    \n    weights = pd.concat(weights) if len(weights) > 1 else weights[0]\n    val_preds = pd.concat(val_preds) if len(val_preds) > 1 else val_preds[0]\n    val_preds.sort_values('val_index', inplace=True)\n    val_preds.reset_index(drop=True, inplace=True)\n    submissions = pd.concat(submissions) if len(submissions) > 1 else submissions[0]\n    \n    return submissions, weights, val_accs, val_preds\n\ncols = []\nsubmissions, weights, val_accs, val_preds = learn({ 'cols': cols, 'random_state': 48 })","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 各特徴の重要度\n\n* 各特徴の目的変数に対する貢献を重要度として表示します．\n* 数値の絶対的な大きさに意味は無く，差を見ます．\n* eli5というライブラリによる算出も試しましたが，処理に時間がかかるので採用せず，，\n* 色々試して曖昧な指標だと判断，700 -> 200の選択にだけ使用して，細かい特徴選択には使いませんでした．"},{"metadata":{"trusted":true},"cell_type":"code","source":"means = weights.groupby('col', as_index=False).mean().rename(columns={'weight':'mean'})\nbars = pd.merge(weights, means, on='col', how='left')\nbars.sort_values('mean', ascending=False, inplace=True)\n\nbars.to_csv('weight.csv', index=False)\nprint('saved weight:', bars.shape)\n\nplt.figure(figsize=(8, 4))\nsns.barplot(x='weight', y='col', data=bars)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 特徴の精査\n\n## 目的変数と特徴の関係を確認\n\n* 特徴が想定通りのデータになっているか，確認を行います．\n* グラフで目視しながら，特徴を選んだり，調整を加えたりします．"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 1, figsize=(16, 4))\nax.plot(train_rows['time_to_failure'])\nax.plot(train_rows['q0.95'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 学習データの正解値と予想値を確認\n\n* 正解値と予測値をグラフに表示して確認します\n* ９つの簡易的な特徴でも，ある程度推論できていることが分かります．\n* グラフにもあるように，地震が長時間起こらないケースが2，３あり，少ない地震からそれを予測するのが難しいです．"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1, 1, figsize=(16, 4))\nax.plot(train_rows['time_to_failure'])\nax.plot(val_preds['time_to_failure'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 相関と分布を確認\n\n* 総当たりで各特徴間の相関を表示しています．\n* 目的変数と各特徴の相関も参考にします．\n* まったく同じ相関や似た相関を持った特徴は一方を除外すると精度が上がることがありました．"},{"metadata":{"trusted":true},"cell_type":"code","source":"colormap = plt.cm.RdBu\nplt.figure(figsize=(10,8))\n\nsns.heatmap(\n    train_rows.drop('seg_id', axis=1).corr(), linewidths=0.1, vmax=1.0, square=True,\n    cmap=colormap, linecolor='white', annot=True)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(train_rows.drop('seg_id', axis=1))\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"\n# 実際の最適化で選んだ特徴の例\n\n* 他にも学習モデルを作ってるので，これらの特徴グループだけではありません．\n* 例えば7番目の特徴は，フーリエ変換虚数部の低周波領域20000の移動平均1000の分位数0.05.\n* ewmは指数移動平均，denoiseはハイパスフィルタを通したものです．\n----\n* mean_mfcc15\n* q0.05_ewm-std30 \n* mean_mfcc4\n* range_mfcc_mean\n* q0.05_denoise_roll-std100\n* q0.05_fft-imag_low20000_roll-std1000\n* q0.05_fft-real_low20000_roll-std1000 \n* q0.05_fft-real_low20000_roll-mean1000\n* peak-count10_ewm-mean3000\n* q0.05_denoise_roll-std10\n* spkt50_denoise\n* std_mfcc3\n* peak-count50_ewm-mean30\n* q0.05_ewm-mean3000\n* mean_mfcc7\n* q0.05_fft-imag_low20000_roll-std100\n* q0.05_roll-mean10\n* q0.01_denoise_roll-mean10\n* q0.2_base\n* kurt_mfcc10\n* peak-count50_denoise_ewm-std30"},{"metadata":{},"cell_type":"markdown","source":"# 他のモデルを作成\n\n* ベースはCatBoostという決定木ベースの勾配ブースティングに基づく機械学習ライブラリで最適化しました．\n* 他にLightGBM，ニューラルネットワークを用いて別の学習モデル作成しましたが，期間が残り少なく上手く精度を上げることができませんでした．\n* 作成した3つのモデルに加えて，カーネルに上がっていた遺伝的アルゴリズムのモデルも使用しています．"},{"metadata":{},"cell_type":"markdown","source":"# 結果を提出\n\n* 以下はクロスバリデーションを平均化しただけの結果を保存しています．"},{"metadata":{"trusted":true},"cell_type":"code","source":"subs = submissions.groupby('seg_id', as_index=False).mean()\nsubs.to_csv('submission.csv', index=False)\nprint('saved submission', subs.shape)\nsubs.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 結果の分布を確認\n\n* 正解，予想，最終結果それぞれの結果の分布を確認してみます．\n* 今回の結果は分布の違いが大きく見えます..\n* 特に10秒以上の長いtime_to_failureを如何に予想するかが，よく議論に挙がっていました．"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.distplot(train_rows['time_to_failure'], kde=True, label='train')\nsns.distplot(val_preds['time_to_failure'], kde=True, label='val')\nsns.distplot(subs['time_to_failure'], kde=True, label='test')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 評価\n\n* MAE（Mean Absolute Error：平均絶対誤差）で評価され，リーダーボードにランキング表示されます．\n\n![MAE](https://wikimedia.org/api/rest_v1/media/math/render/svg/3ef87b78a9af65e308cf4aa9acf6f203efbdeded)\n\n### クロスバリデーションによる学習データでの評価\n\n* CatBoostで学習を進めるとどんどん下がっていくものの，テストデータによるランキング（公開LB）が下がってしまうので過学習と判断．\n* 過学習していると思われる結果を除くと1.97〜8程度で頭打ちに..\n\n### 公開リーダーボード（一部テストデータによる暫定ランキング）\n\n* シェイクアップ前の公開LB\n* 最初の頃からクロスバリデーションでの結果よりかなり低い値になっていました（CV2.0で公開LB1.5とか）\n* 学習データをシャッフルして推論すると稀に同じような結果になるので，評価に使われているテストデータの一部がかなり偏っていそう．\n* シングルモデルが1.43〜４辺りで頭打ちに．\n* 幾つか試すとたまに精度が上がるも，シード値やデータを少しいいじると下がってしまうのでラッキーな過学習だと判断．\n* ディスカッションでもランキング（プライベートLB）確定時に激しい順位のシェイクアップが起こりそうとの噂..\n* 異なるアプローチで最適化された学習モデルを平均して提出すると精度が上がるらしく，実践してみました．\n* 0開始と75,000開始で分割した学習データ２パターン×３モデルで６つの学習モデルを作成．\n* 遺伝的アルゴリズムモデル，少し前に作って精度が良かったモデルを合わせて8つの平均値を最終提出．\n* ベストモデルの公開LBの評価は「1.38223」で205位．\n\n### プライベートリーダーボード（シェイクアップ後のランキング確定）\n\n* 公開リーダーボードの評価は「2.4514」で77位（モデルを2つ提出できるので，公開LBとは違うモデルがランクイン）\n* 提出履歴を眺めて見ると，↑のほうに示したように特徴量を20個程度に絞った学習モデルが高い評価を出していました．\n* 提出履歴での最高は「2.4074」で順位にすると21位であるものの，そのモデルの良し悪しをコンペ中に判断できなかったわけで意味無し．\n\n# 上位入賞者のアプローチ\n\n* https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion\n* 勉強が足りず分からないことが多いものの，以下のような検証が必要だった模様．\n  * 音響データでの地震発生から時間リセットまでの外れ値の扱い\n  * 学習データとテストデータの特徴の分布や予想結果の異なりを検証\n    * 敵対的検証\n    * 学習データを繰り返しサンプリングしてKS検定\n    * 学習データにノイズを加える"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}