{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Baseline RF model: Reproducing the 2017 Rouet-Leduc et al. paper\n\n_ONGOING WORK, FEATURE ENGINEERING STILL INCOMPLETE_\n\nWe read: \"_The competition builds on initial work from Bertrand Rouet-Leduc, Claudia Hulbert, and Paul Johnson. B. Rouet-Leduc prepared the data for the competition._\" Their original article is:\n\nRouet-leduc, B., Hulbert, C., Lubbers, N., Barros, K., Humphreys, C.J., Johnson, P.A. (2017), Machine Learning Predicts Laboratory Earthquakes. Geophys. Res. Lett., 44, 9276-9282, doi: [10.1002/2017GL074677](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017GL074677)\n\nIn the introduction of the competition's discussion section, Bertrand Rouet-Leduc gives us more details and links to three additional articles. Compared to what has been published, he mentions that, for this Kaggle competition, they selected \"_an experiment that exhibits a **very aperiodic** and more realistic behavior compared to the data we studied in our early work, with **earthquakes occurring very irregularly**._\" He adds on a second post that \"_you will have to build a model that predicts the time remaining before failure from a chunk of seismic data, like we have done in our first paper above on easier data._\"\n\nAfter a quick search on the discussion and kernel list, I could not find any obvious link to a work trying to reproduce the 2017 Rouet-Leduc article. I will therefore summarize their work below and try to apply their method to the LANL Kaggle dataset. It is very likely that others already did it (please comment if so) but I will here focus on this 2017 study and describe it in detail, which I believe might be of interest to some Kagglers.\n\n\n\n## I. The Rouet-Leduc et al. (2017) study\n\n### 1. Abstract\n\nI here simply copy-paste their abstract from the GRL website:\n\n\"_We apply machine learning to data sets from shear laboratory experiments, with the goal of identifying hidden signals that precede earthquakes. Here we show that by listening to the acoustic signal emitted by a laboratory fault, machine learning can predict the time remaining before it fails with great accuracy. These predictions are based solely on the instantaneous physical characteristics of the acoustical signal and do not make use of its history. Surprisingly, machine learning identifies a signal emitted from the fault zone previously thought to be low-amplitude noise that enables failure forecasting throughout the laboratory quake cycle. We infer that this signal originates from continuous grain motions of the fault gouge as the fault blocks displace. We posit that applying this approach to continuous seismic data may lead to significant advances in identifying currently unknown signals, in providing new insights into fault physics, and in placing bounds on fault failure times._\"\n\n\n### 2. Method\n\nRouet-Leduc et al. (2017) applied the **Random Forest algorithm** to the continuous acoustic time series data. For each time window, the authors computed a set of **c. 100 potentially relevant statistical features (incl. mean, variance, kurtosis & autocorrelation)**. Then, they used **recursive feature elimination** for feature reduction.\n\n#### a. Random Forest\n\nTo apply a Random Forest, we will need to define several hyperparameters:"},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\n\n#model_RF = RandomForestRegressor(n_estimators = n_estimators,     # number of trees in the forest\n#                                 criterion = criterion,           # quality of split measure\n#                                 max_depth = max_depth,\n#                                 max_features = max_features)     # nb. of features to consider for best split","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We learn from their Fig. 1 that \"_the RF model predicts the time remaining before the next failure by averaging the predictions of 1,000 decision trees for each time window._\" which gives us **n_estimators = 1000**. In Fig. 2, we learn that the RF was trained of c. 150 s of data (c. 10 slip events) and tested on the following c. 150 s.\n\nNo more details are given in the main text and we have to turn to the supplementary material to get more details about the engineered features and hyperparameters:\n\nAt each node, they selected a random subset of 40% of the available features, which gives us **max_features = 0.4**. We learn later on that this value had been selected via 3-fold cross-validation (from a range 30-40%).\n\nAs split criterion, they used maximum reduction of empirical variance, which is **criterion = 'mse'** (the default in RandomForestRegressor). No information was given about maximum depth.\n\n#### b. Feature engineering\n\nThe authors used moving time windows of 1.8 s each, with an offset of 0.18 s (i.e. 90% overlap between consecutive windows). Each window corresponds to one observation X_i (feature vector) and label y_i (time remaining until next failure).\n\nHere is the list of features:\n- Signal distribution (to capture evolution of signal's energy) = 7 features: **kth-order moment** of acoustic data corresponding to mean, normalized/non-normalized variance, skewness, kurtosis;\n- Precursors (bursts of acoustic emissions) = 18+10+2 = 30 features: **1st-to-9th and 91th-to-99th percentiles** (per 1% increment), the fraction of times the strain is greater than **thresholds f_0 derived from signal processing** with **f_0 = 1e-9, 5e-9, 1e-8, 5.e-8 and 1e-7** and the fraction of times it is lower than their negatives -f_0, and finally **min/max** values;\n- Time correlation = 5 features: **Fourier transforms** on frequency bands {(19.65, 20.65), (39.8, 40.8), (80.1, 81.1)} in kHz, the **autocorrelation** and the **partial autocorrelation**.\n\nFeatures were originally defined from the raw data and the **first finite difference of the data**. The authors found that the RF had a **slight performance advantage when only using the derivative signal**.\n\nWe verify that we get (7+30+5) times 2 = **84 features**. The doubling comes from a **window split**, which gave the algorithm a notion of **short-term evolution of the signal**. It means that the 42 features were estimated twice for the two subwindows of each \"observation\". We will describe in detail these different features in part II.\n\n\n### 3. Results\n\nRouet-Leduc et al. (2017) obtained R^2 = 0.89 for the RF model compared to R^2 = 0.3 for a naive model based on event periodicity, which means that the RF model explains 89% of the data variance. No mean absolute error (MAE) was given in the paper, which hampers direct comparison with the Kaggle competition leaderboard.\n\n<br><br>\n\n\n\n## II. Comparison with the LAN Kaggle competition\n\nComparison of their experiment with the Kaggle LAN competition. I once again took some inspiration from Grand Kernel Master [Andrew Lukyanenko](https://www.kaggle.com/artgor/earthquakes-fe-more-features-and-samples) for EDA.\n\n### 1. LANL EDA\n\nThe training data set is c. 10 GB"},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.options.display.precision = 15\nimport scipy\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### a. Training set"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.dtypes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.tail()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_acoustic_data_small = train['acoustic_data'].values[::100]      #skip 100 rows\ntrain_time_to_failure_small = train['time_to_failure'].values[::100]\n\ntf = train.loc[np.diff(train['time_to_failure']) > 0].index.values\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title('Smoothed training data')\nplt.plot(train_acoustic_data_small, color = 'b')\nfor ti in tf:\n    plt.axvline(x = ti/100, color = 'r', linestyle = '--')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\n#del train_acoustic_data_small\n#del train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"p = 0.01   # first p percentage of training data\nindp = round(629145480*p)\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title('First 1% of training data')\nplt.plot(train['acoustic_data'].values[:indp], color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train['time_to_failure'].values[:indp], color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### b. Test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv')\nseg_id = submission['seg_id']\nn_test = len(seg_id)\nn_test","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test0 = pd.read_csv('../input/test/' + seg_id[0] + '.csv')\nnp.shape(test0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.subplots(figsize=(16, 8))\nplt.title('First test sample')\nplt.plot(test0['acoustic_data'], color='b')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.subplots(figsize=(16, 22))\nfor i in range(48):\n    testi = pd.read_csv('../input/test/' + seg_id[i] + '.csv')\n    plt.subplot(6, 8, i + 1)\n    plt.plot(testi['acoustic_data'])\n    plt.title(seg_id[i])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. Feature engineering\n\nRouet-Leduc et al. (2017) sampled their data into time windows of 1.8 s each with an offset of 0.18 s. However we have no information about time binning in the test data. Moreover the time bin in the training set is irregular. Therefore we cannot reproduce the same sampling. The simplest option is therefore to sample the training set in samples of 150,000 successive records to match the test set samples."},{"metadata":{"trusted":true},"cell_type":"code","source":"records = 150000\nn_train = int(np.floor(train.shape[0] / records))\nn_train","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This seems reasonable, with c. 60% of observations in training set and 40% in test set. Let's now define the same features as in the 2017 study.\n\n#### Define features for training set"},{"metadata":{"trusted":true},"cell_type":"code","source":"f0pos = (1e-9, 5e-9, 1e-8, 5e-8, 1e-7)\nf0neg = (-1e-9, -5e-9, -1e-8, -5e-8, -1e-7)\n\nX_train = pd.DataFrame(index = range(n_train), columns = \n                    ['meanA', 'varA', 'varAnorm', 'skewA', 'skewAnorm', 'kurtA', 'kurtAnorm',\n                     'meanB', 'varB', 'varBnorm', 'skewB', 'skewBnorm', 'kurtB', 'kurtBnorm',\n                     'q01A', 'q02A', 'q03A', 'q04A', 'q05A', 'q06A', 'q07A', 'q08A', 'q09A',\n                     'q01B', 'q02B', 'q03B', 'q04B', 'q05B', 'q06B', 'q07B', 'q08B', 'q09B',\n                     'q91A', 'q92A', 'q93A', 'q94A', 'q95A', 'q96A', 'q97A', 'q98A', 'q99A',\n                     'q91B', 'q92B', 'q93B', 'q94B', 'q95B', 'q96B', 'q97B', 'q98B', 'q99B',\n                     'f00pA', 'f01pA', 'f02pA', 'f03pA', 'f04pA', 'f00nA', 'f01nA', 'f02nA', 'f03nA', 'f04nA',\n                     'f00pB', 'f01pB', 'f02pB', 'f03pB', 'f04pB', 'f00nB', 'f01nB', 'f02nB', 'f03nB', 'f04nB',\n                     'minA', 'maxA', 'minB', 'maxB'\n                    ])\n\ny_train = pd.DataFrame(index = range(n_train), columns = ['time_to_failure'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Not sure how the higher moments should be normalized. At the present time, the time-correlation set of features from Rouet-Leduc et al. (2017) are not included since the authors mentioned that they \"_only very marginally improve the predictions._\" Also, I only use the raw data for now since the authors obtained only slightly better results with the first finite difference of the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(n_train):\n    if i % 100 == 0:\n        print(i)\n    \n    segment = train.iloc[i*records : i*records + records]\n    y_train.loc[i, 'time_to_failure'] = segment['time_to_failure'].values[-1]\n\n    X_train.loc[i, 'meanA'] = segment['acoustic_data'][0:round(records/2)].mean()\n    X_train.loc[i, 'meanB'] = segment['acoustic_data'][round(records/2)+1:records].mean()\n    X_train.loc[i, 'varA'] = segment['acoustic_data'][0:round(records/2)].var()\n    X_train.loc[i, 'varB'] = segment['acoustic_data'][round(records/2)+1:records].var()\n    X_train.loc[i, 'skewA'] = scipy.stats.skew(segment['acoustic_data'][0:round(records/2)])\n    X_train.loc[i, 'skewB'] = scipy.stats.skew(segment['acoustic_data'][round(records/2)+1:records])\n    X_train.loc[i, 'kurtA'] = scipy.stats.kurtosis(segment['acoustic_data'][0:round(records/2)])\n    X_train.loc[i, 'kurtB'] = scipy.stats.kurtosis(segment['acoustic_data'][round(records/2)+1:records])\n    X_train.loc[i, 'varAnorm'] = X_train.loc[i, 'varA']/(X_train.loc[i, 'varA']+X_train.loc[i, 'varB'])\n    X_train.loc[i, 'varBnorm'] = X_train.loc[i, 'varB']/(X_train.loc[i, 'varA']+X_train.loc[i, 'varB'])\n    X_train.loc[i, 'skewAnorm'] = X_train.loc[i, 'skewA']/(X_train.loc[i, 'skewA']+X_train.loc[i, 'skewB'])\n    X_train.loc[i, 'skewBnorm'] = X_train.loc[i, 'skewB']/(X_train.loc[i, 'skewA']+X_train.loc[i, 'skewB'])\n    X_train.loc[i, 'kurtAnorm'] = X_train.loc[i, 'kurtA']/(X_train.loc[i, 'kurtA']+X_train.loc[i, 'kurtB'])\n    X_train.loc[i, 'kurtBnorm'] = X_train.loc[i, 'kurtB']/(X_train.loc[i, 'kurtA']+X_train.loc[i, 'kurtB'])\n    \n    X_train.loc[i, 'q01A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.01)\n    X_train.loc[i, 'q02A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.02)\n    X_train.loc[i, 'q03A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.03)\n    X_train.loc[i, 'q04A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.04)\n    X_train.loc[i, 'q05A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.05)\n    X_train.loc[i, 'q06A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.06)\n    X_train.loc[i, 'q07A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.07)\n    X_train.loc[i, 'q08A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.08)\n    X_train.loc[i, 'q09A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.09)\n    X_train.loc[i, 'q91A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.91)\n    X_train.loc[i, 'q92A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.92)\n    X_train.loc[i, 'q93A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.93)\n    X_train.loc[i, 'q94A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.94)\n    X_train.loc[i, 'q95A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.95)\n    X_train.loc[i, 'q96A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.96)\n    X_train.loc[i, 'q97A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.97)\n    X_train.loc[i, 'q98A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.98)\n    X_train.loc[i, 'q99A'] = np.quantile(segment['acoustic_data'][0:round(records/2)], 0.99)\n    X_train.loc[i, 'q01B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.01)\n    X_train.loc[i, 'q02B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.02)\n    X_train.loc[i, 'q03B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.03)\n    X_train.loc[i, 'q04B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.04)\n    X_train.loc[i, 'q05B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.05)\n    X_train.loc[i, 'q06B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.06)\n    X_train.loc[i, 'q07B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.07)\n    X_train.loc[i, 'q08B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.08)\n    X_train.loc[i, 'q09B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.09)\n    X_train.loc[i, 'q91B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.91)\n    X_train.loc[i, 'q92B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.92)\n    X_train.loc[i, 'q93B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.93)\n    X_train.loc[i, 'q94B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.94)\n    X_train.loc[i, 'q95B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.95)\n    X_train.loc[i, 'q96B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.96)\n    X_train.loc[i, 'q97B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.97)\n    X_train.loc[i, 'q98B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.98)\n    X_train.loc[i, 'q99B'] = np.quantile(segment['acoustic_data'][round(records/2)+1:records], 0.99)\n\n    X_train.loc[i, 'f00pA'] = sum(segment['acoustic_data'][0:round(records/2)] >= f0pos[0])/75000\n    X_train.loc[i, 'f01pA'] = sum(segment['acoustic_data'][0:round(records/2)] >= f0pos[1])/75000\n    X_train.loc[i, 'f02pA'] = sum(segment['acoustic_data'][0:round(records/2)] >= f0pos[2])/75000\n    X_train.loc[i, 'f03pA'] = sum(segment['acoustic_data'][0:round(records/2)] >= f0pos[3])/75000\n    X_train.loc[i, 'f04pA'] = sum(segment['acoustic_data'][0:round(records/2)] >= f0pos[4])/75000\n    X_train.loc[i, 'f00nA'] = sum(segment['acoustic_data'][0:round(records/2)] <= f0neg[0])/75000\n    X_train.loc[i, 'f01nA'] = sum(segment['acoustic_data'][0:round(records/2)] <= f0neg[1])/75000\n    X_train.loc[i, 'f02nA'] = sum(segment['acoustic_data'][0:round(records/2)] <= f0neg[2])/75000\n    X_train.loc[i, 'f03nA'] = sum(segment['acoustic_data'][0:round(records/2)] <= f0neg[3])/75000\n    X_train.loc[i, 'f04nA'] = sum(segment['acoustic_data'][0:round(records/2)] <= f0neg[4])/75000\n    X_train.loc[i, 'f00pB'] = sum(segment['acoustic_data'][round(records/2)+1:records] >= f0pos[0])/74999\n    X_train.loc[i, 'f01pB'] = sum(segment['acoustic_data'][round(records/2)+1:records] >= f0pos[1])/74999\n    X_train.loc[i, 'f02pB'] = sum(segment['acoustic_data'][round(records/2)+1:records] >= f0pos[2])/74999\n    X_train.loc[i, 'f03pB'] = sum(segment['acoustic_data'][round(records/2)+1:records] >= f0pos[3])/74999\n    X_train.loc[i, 'f04pB'] = sum(segment['acoustic_data'][round(records/2)+1:records] >= f0pos[4])/74999\n    X_train.loc[i, 'f00nB'] = sum(segment['acoustic_data'][round(records/2)+1:records] <= f0neg[0])/74999\n    X_train.loc[i, 'f01nB'] = sum(segment['acoustic_data'][round(records/2)+1:records] <= f0neg[1])/74999\n    X_train.loc[i, 'f02nB'] = sum(segment['acoustic_data'][round(records/2)+1:records] <= f0neg[2])/74999\n    X_train.loc[i, 'f03nB'] = sum(segment['acoustic_data'][round(records/2)+1:records] <= f0neg[3])/74999\n    X_train.loc[i, 'f04nB'] = sum(segment['acoustic_data'][round(records/2)+1:records] <= f0neg[4])/74999\n    \n    X_train.loc[i, 'minA'] = min(segment['acoustic_data'][0:round(records/2)])\n    X_train.loc[i, 'maxA'] = max(segment['acoustic_data'][0:round(records/2)])\n    X_train.loc[i, 'minB'] = min(segment['acoustic_data'][round(records/2)+1:records])\n    X_train.loc[i, 'maxB'] = max(segment['acoustic_data'][round(records/2)+1:records])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let us produce a figure similar to Fig. 1b of Rouet-Leduc et al. (2017):"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16,8))\nplt.subplot(411)\nplt.title('Smoothed training data')\nplt.plot(train_acoustic_data_small)\nfor ti in tf:\n    plt.axvline(x = ti/100, color = 'r', linestyle = '--')\n\nplt.subplot(412)\nplt.title('Mean')\nplt.plot(X_train['meanA'])\n\nplt.subplot(413)\nplt.title('Skewness')\nplt.plot(X_train['skewA'])\n\nplt.subplot(414)\nplt.title('Kurtosis')\nplt.plot(X_train['kurtA'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Define same features for test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test = pd.DataFrame(index = range(n_test), columns = \n                    ['meanA', 'varA', 'varAnorm', 'skewA', 'skewAnorm', 'kurtA', 'kurtAnorm',\n                     'meanB', 'varB', 'varBnorm', 'skewB', 'skewBnorm', 'kurtB', 'kurtBnorm',\n                     'q01A', 'q02A', 'q03A', 'q04A', 'q05A', 'q06A', 'q07A', 'q08A', 'q09A',\n                     'q01B', 'q02B', 'q03B', 'q04B', 'q05B', 'q06B', 'q07B', 'q08B', 'q09B',\n                     'q91A', 'q92A', 'q93A', 'q94A', 'q95A', 'q96A', 'q97A', 'q98A', 'q99A',\n                     'q91B', 'q92B', 'q93B', 'q94B', 'q95B', 'q96B', 'q97B', 'q98B', 'q99B',\n                     'f00pA', 'f01pA', 'f02pA', 'f03pA', 'f04pA', 'f00nA', 'f01nA', 'f02nA', 'f03nA', 'f04nA',\n                     'f00pB', 'f01pB', 'f02pB', 'f03pB', 'f04pB', 'f00nB', 'f01nB', 'f02nB', 'f03nB', 'f04nB',\n                     'minA', 'maxA', 'minB', 'maxB'\n                    ])\n\nfor i in range(n_test):\n    if i % 100 == 0:\n        print(i)\n\n    testi = pd.read_csv('../input/test/' + seg_id[i] + '.csv')\n    X_test.loc[i, 'meanA'] = testi['acoustic_data'][0:round(records/2)].mean()\n    X_test.loc[i, 'meanB'] = testi['acoustic_data'][round(records/2)+1:records].mean()\n    X_test.loc[i, 'varA'] = testi['acoustic_data'][0:round(records/2)].var()\n    X_test.loc[i, 'varB'] = testi['acoustic_data'][round(records/2)+1:records].var()\n    X_test.loc[i, 'skewA'] = scipy.stats.skew(testi['acoustic_data'][0:round(records/2)])\n    X_test.loc[i, 'skewB'] = scipy.stats.skew(testi['acoustic_data'][round(records/2)+1:records])\n    X_test.loc[i, 'kurtA'] = scipy.stats.kurtosis(testi['acoustic_data'][0:round(records/2)])\n    X_test.loc[i, 'kurtB'] = scipy.stats.kurtosis(testi['acoustic_data'][round(records/2)+1:records])\n    X_test.loc[i, 'varAnorm'] = X_test.loc[i, 'varA']/(X_test.loc[i, 'varA']+X_test.loc[i, 'varB'])\n    X_test.loc[i, 'varBnorm'] = X_test.loc[i, 'varB']/(X_test.loc[i, 'varA']+X_test.loc[i, 'varB'])\n    X_test.loc[i, 'skewAnorm'] = X_test.loc[i, 'skewA']/(X_test.loc[i, 'skewA']+X_test.loc[i, 'skewB'])\n    X_test.loc[i, 'skewBnorm'] = X_test.loc[i, 'skewB']/(X_test.loc[i, 'skewA']+X_test.loc[i, 'skewB'])\n    X_test.loc[i, 'kurtAnorm'] = X_test.loc[i, 'kurtA']/(X_test.loc[i, 'kurtA']+X_test.loc[i, 'kurtB'])\n    X_test.loc[i, 'kurtBnorm'] = X_test.loc[i, 'kurtB']/(X_test.loc[i, 'kurtA']+X_test.loc[i, 'kurtB'])\n\n    X_test.loc[i, 'q01A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.01)\n    X_test.loc[i, 'q02A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.02)\n    X_test.loc[i, 'q03A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.03)\n    X_test.loc[i, 'q04A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.04)\n    X_test.loc[i, 'q05A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.05)\n    X_test.loc[i, 'q06A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.06)\n    X_test.loc[i, 'q07A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.07)\n    X_test.loc[i, 'q08A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.08)\n    X_test.loc[i, 'q09A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.09)\n    X_test.loc[i, 'q91A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.91)\n    X_test.loc[i, 'q92A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.92)\n    X_test.loc[i, 'q93A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.93)\n    X_test.loc[i, 'q94A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.94)\n    X_test.loc[i, 'q95A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.95)\n    X_test.loc[i, 'q96A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.96)\n    X_test.loc[i, 'q97A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.97)\n    X_test.loc[i, 'q98A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.98)\n    X_test.loc[i, 'q99A'] = np.quantile(testi['acoustic_data'][0:round(records/2)], 0.99)\n    X_test.loc[i, 'q01B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.01)\n    X_test.loc[i, 'q02B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.02)\n    X_test.loc[i, 'q03B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.03)\n    X_test.loc[i, 'q04B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.04)\n    X_test.loc[i, 'q05B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.05)\n    X_test.loc[i, 'q06B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.06)\n    X_test.loc[i, 'q07B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.07)\n    X_test.loc[i, 'q08B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.08)\n    X_test.loc[i, 'q09B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.09)\n    X_test.loc[i, 'q91B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.91)\n    X_test.loc[i, 'q92B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.92)\n    X_test.loc[i, 'q93B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.93)\n    X_test.loc[i, 'q94B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.94)\n    X_test.loc[i, 'q95B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.95)\n    X_test.loc[i, 'q96B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.96)\n    X_test.loc[i, 'q97B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.97)\n    X_test.loc[i, 'q98B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.98)\n    X_test.loc[i, 'q99B'] = np.quantile(testi['acoustic_data'][round(records/2)+1:records], 0.99)\n\n    X_test.loc[i, 'f00pA'] = sum(testi['acoustic_data'][0:round(records/2)] >= f0pos[0])/75000\n    X_test.loc[i, 'f01pA'] = sum(testi['acoustic_data'][0:round(records/2)] >= f0pos[1])/75000\n    X_test.loc[i, 'f02pA'] = sum(testi['acoustic_data'][0:round(records/2)] >= f0pos[2])/75000\n    X_test.loc[i, 'f03pA'] = sum(testi['acoustic_data'][0:round(records/2)] >= f0pos[3])/75000\n    X_test.loc[i, 'f04pA'] = sum(testi['acoustic_data'][0:round(records/2)] >= f0pos[4])/75000\n    X_test.loc[i, 'f00nA'] = sum(testi['acoustic_data'][0:round(records/2)] <= f0neg[0])/75000\n    X_test.loc[i, 'f01nA'] = sum(testi['acoustic_data'][0:round(records/2)] <= f0neg[1])/75000\n    X_test.loc[i, 'f02nA'] = sum(testi['acoustic_data'][0:round(records/2)] <= f0neg[2])/75000\n    X_test.loc[i, 'f03nA'] = sum(testi['acoustic_data'][0:round(records/2)] <= f0neg[3])/75000\n    X_test.loc[i, 'f04nA'] = sum(testi['acoustic_data'][0:round(records/2)] <= f0neg[4])/75000\n    X_test.loc[i, 'f00pB'] = sum(testi['acoustic_data'][round(records/2)+1:records] >= f0pos[0])/74999\n    X_test.loc[i, 'f01pB'] = sum(testi['acoustic_data'][round(records/2)+1:records] >= f0pos[1])/74999\n    X_test.loc[i, 'f02pB'] = sum(testi['acoustic_data'][round(records/2)+1:records] >= f0pos[2])/74999\n    X_test.loc[i, 'f03pB'] = sum(testi['acoustic_data'][round(records/2)+1:records] >= f0pos[3])/74999\n    X_test.loc[i, 'f04pB'] = sum(testi['acoustic_data'][round(records/2)+1:records] >= f0pos[4])/74999\n    X_test.loc[i, 'f00nB'] = sum(testi['acoustic_data'][round(records/2)+1:records] <= f0neg[0])/74999\n    X_test.loc[i, 'f01nB'] = sum(testi['acoustic_data'][round(records/2)+1:records] <= f0neg[1])/74999\n    X_test.loc[i, 'f02nB'] = sum(testi['acoustic_data'][round(records/2)+1:records] <= f0neg[2])/74999\n    X_test.loc[i, 'f03nB'] = sum(testi['acoustic_data'][round(records/2)+1:records] <= f0neg[3])/74999\n    X_test.loc[i, 'f04nB'] = sum(testi['acoustic_data'][round(records/2)+1:records] <= f0neg[4])/74999\n    \n    X_test.loc[i, 'minA'] = min(testi['acoustic_data'][0:round(records/2)])\n    X_test.loc[i, 'maxA'] = max(testi['acoustic_data'][0:round(records/2)])\n    X_test.loc[i, 'minB'] = min(testi['acoustic_data'][round(records/2)+1:records])\n    X_test.loc[i, 'maxB'] = max(testi['acoustic_data'][round(records/2)+1:records])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_test.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 3. Modelling\n\nAs in the 2017 paper, we will use Random Forest (see parameterization in I.2.a)."},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nfrom sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor\nfrom sklearn.metrics import mean_absolute_error\n\nfrom catboost import CatBoostRegressor","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Xsub_train, Xsub_val, ysub_train, ysub_val = train_test_split(X_train, y_train, test_size=0.4)\n\nmodel_RF = RandomForestRegressor(n_estimators = 1000,\n                                 criterion = 'mse',\n                                 max_features = 0.4)\n\nmodel_RF.fit(Xsub_train, ysub_train)\nmodel_RF_pred_valset = model_RF.predict(Xsub_val)\nMAE = mean_absolute_error(ysub_val, model_RF_pred_valset)\nMAE","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since Boosting often outperforms Random Forests, I also try CatBoost to slightly improve my score."},{"metadata":{"trusted":true},"cell_type":"code","source":"### try CatBoost?\n#Xsub_train, Xsub_val, ysub_train, ysub_val = train_test_split(X_train, y_train, test_size=0.4)\n\nmodel_CatBoost = CatBoostRegressor(silent=True)\n\nmodel_CatBoost.fit(Xsub_train, ysub_train)\nmodel_CatBoost_pred_valset = model_CatBoost.predict(Xsub_val)\nMAE = mean_absolute_error(ysub_val, model_CatBoost_pred_valset)\nMAE","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 4. Submission\n\n#### Random Forest"},{"metadata":{"trusted":true},"cell_type":"code","source":"#rerun model on full training set and predict on test set\nmodel_RF = RandomForestRegressor(n_estimators = 1000,\n                                 criterion = 'mse',\n                                 max_features = 0.4)\n\nmodel_RF.fit(X_train, y_train)\nmodel_RF_pred_testset = model_RF.predict(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission['time_to_failure'] = model_RF_pred_testset\nsubmission.to_csv('submission_RF.csv', index = False)\nsubmission.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_CatBoost = CatBoostRegressor(silent=True)\n\nmodel_CatBoost.fit(X_train, y_train)\nmodel_CatBoost_pred_testset = model_CatBoost.predict(X_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission['time_to_failure'] = model_CatBoost_pred_testset\nsubmission.to_csv('submission_CatBoost.csv', index = False)\nsubmission.head()","execution_count":null,"outputs":[]}],"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}