{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 1. Introduction","metadata":{"_uuid":"59da319c-8a3d-453f-ac00-362ab21a5789","_cell_guid":"e8e0abec-b767-4a9b-a19b-a436644b328f","trusted":true}},{"cell_type":"markdown","source":"<font color=#6698FF>Name: Tomasz Abels and Jack Chen\n\n<font color=#6698FF>Username: JackChenXJ\n\n<font color=#6698FF>Score: 2.59039","metadata":{"_uuid":"f0673542-b16f-4ea3-9053-12bb852ba20e","_cell_guid":"3a176642-ab7f-4473-8d9b-834e8e0025aa","trusted":true}},{"cell_type":"markdown","source":"# 2. Data","metadata":{"_uuid":"23a8c03a-a960-44de-84a0-1972160f3e8e","_cell_guid":"bc110290-d051-410a-bdf8-246514dfd454","trusted":true}},{"cell_type":"markdown","source":"## 2.1 Dataset\nIn this section, the dataset and the test set will be loaded and explored.  ","metadata":{"_uuid":"69dd81c8-8529-4d1b-9b29-469f21c5bbb1","_cell_guid":"ec14d3cd-6a85-42b6-ae82-cc1260073b0a","trusted":true}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt \n\nimport os\nprint(os.listdir(\"../input/LANL-Earthquake-Prediction\"))","metadata":{"_uuid":"cd2d6ef0-f99f-4024-b2e2-83da6c143176","_cell_guid":"2d3ea610-6dc3-4997-bb3d-2e06dfed17b5","execution":{"iopub.status.busy":"2023-04-23T19:01:51.610430Z","iopub.execute_input":"2023-04-23T19:01:51.610816Z","iopub.status.idle":"2023-04-23T19:01:51.639982Z","shell.execute_reply.started":"2023-04-23T19:01:51.610773Z","shell.execute_reply":"2023-04-23T19:01:51.639079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64}) # read the csv data file\nprint(train.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:01:51.641709Z","iopub.execute_input":"2023-04-23T19:01:51.642340Z","iopub.status.idle":"2023-04-23T19:05:28.874186Z","shell.execute_reply.started":"2023-04-23T19:01:51.642299Z","shell.execute_reply":"2023-04-23T19:05:28.872203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acoustic_data = train['acoustic_data'].values[::100] # sample at every 100th data point\ntime_data = train['time_to_failure'].values[::100] # sample at every 100th data point\nprint(acoustic_data.shape)\nprint(time_data.shape)","metadata":{"_uuid":"edab80fb-df72-4d0c-91ba-fdf9cab68dda","_cell_guid":"e45a5cc2-1b26-416a-bf79-b78d4a62033a","execution":{"iopub.status.busy":"2023-04-23T19:05:28.876539Z","iopub.execute_input":"2023-04-23T19:05:28.876923Z","iopub.status.idle":"2023-04-23T19:05:28.892179Z","shell.execute_reply.started":"2023-04-23T19:05:28.876882Z","shell.execute_reply":"2023-04-23T19:05:28.890948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Dataset","metadata":{}},{"cell_type":"code","source":"# plot the data set\nfig, ax1 = plt.subplots(figsize=(16, 4))\nplt.title(\"Acoustic data and the time to failure. Sampled every 100th data point\")\nplt.plot(acoustic_data, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(time_data, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.8))","metadata":{"_uuid":"89f611e5-d9bf-41f5-8eaa-f6659a44525e","_cell_guid":"f9340c17-4547-49a7-aa51-adc2f73fb852","execution":{"iopub.status.busy":"2023-04-23T19:05:28.895155Z","iopub.execute_input":"2023-04-23T19:05:28.895605Z","iopub.status.idle":"2023-04-23T19:05:36.339356Z","shell.execute_reply.started":"2023-04-23T19:05:28.895563Z","shell.execute_reply":"2023-04-23T19:05:36.338105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test data","metadata":{}},{"cell_type":"code","source":"# First three test data\ntest1 = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/seg_00030f.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\ntest2 = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/seg_0012b5.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\ntest3 = pd.read_csv('/kaggle/input/LANL-Earthquake-Prediction/test/seg_00184e.csv' ,dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\n\n# Initialize the subplots\nfig, ax = plt.subplots(3, 1, figsize=(16,6))\n\n# Plot the time domain, label the graph and limit the axis\nax[0].plot(test1, color='b')\nax[0].legend(['acoustic_data'])\nax[0].set_ylabel('acoustic_data')\n\nax[1].plot(test2, color='b')\nax[1].legend(['acoustic_data'])\nax[1].set_ylabel('acoustic_data')\n\nax[2].plot(test3, color='b')\nax[2].legend(['acoustic_data'])\nax[2].set_ylabel('acoustic_data')","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:05:36.340609Z","iopub.execute_input":"2023-04-23T19:05:36.341008Z","iopub.status.idle":"2023-04-23T19:05:37.744111Z","shell.execute_reply.started":"2023-04-23T19:05:36.340962Z","shell.execute_reply":"2023-04-23T19:05:37.743069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Feature extraction\nThe features that are used for this project mainly come from other people, which are referenced in the last section. Basic statistical measures are used such as, mean, standard deviation, minimum and maximum. These are also applied on smaller intervals. Also the quantiles of various intervals are calulated. The rolling window is used to determine the statistical measures on a certain moving interval over the whole data set. \n\nAfter the features are extracted, the data will be normalized and then PCA will be performed. The PCA histogram shows that a certain feature is the best feature. This feature will be plotted with the time to failure data. It can be seen in the plot that the value of the feature peaks when the time to failure data also peaks, which indicates that the feature is a good feature. \n\nFeatures such as race and sex can be used for the machine learning algorithm, but it depends on whether the features are good and on the situation. If these features are much better than the rest, it can be considered to use them, but if there are not, these features can be left out of the feature pool.  ","metadata":{"_uuid":"7a960c75-32e6-4385-bd6e-a46fada82908","_cell_guid":"44058ade-c4bc-4a5b-a951-37ba9a403166","trusted":true}},{"cell_type":"markdown","source":"### Feature extraction","metadata":{}},{"cell_type":"code","source":"# Code from Basic Feature Benchmark by INVERSION, https://www.kaggle.com/code/inversion/basic-feature-benchmark\n# Code from Earthquakes FE. More features and samples by ANDREW LUKYANENKO, https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\n# Code from Rolling + Quantiles by PATRICK YAM, https://www.kaggle.com/code/wimwim/rolling-quantiles/notebook\n\nfrom tqdm import tqdm\nfrom scipy import stats\n\nrows = 150000 \nsegments = int(np.floor(train.shape[0] / rows)) # divide the data into smaller segments\n\nX_train = pd.DataFrame(index=range(segments), dtype=np.float64)\ny_train = pd.DataFrame(index=range(segments), dtype=np.float64,columns=['time_to_failure'])\n\nfor segment in tqdm(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n    x = pd.Series(seg['acoustic_data'].values)\n    y = seg['time_to_failure'].values[-1]\n    \n    y_train.loc[segment, 'time_to_failure'] = y\n    \n    # Basic statistical measures \n    X_train.loc[segment, 'ave'] = x.mean()\n    X_train.loc[segment, 'std'] = x.std()\n    X_train.loc[segment, 'max'] = x.max()\n    X_train.loc[segment, 'min'] = x.min()\n    \n    X_train.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x))\n    X_train.loc[segment, 'abs_max'] = np.abs(x).max()\n    X_train.loc[segment, 'abs_min'] = np.abs(x).min()\n    \n    # Basic statistical measures with intervals\n    X_train.loc[segment, 'std_first_50000'] = x[:50000].std()\n    X_train.loc[segment, 'std_last_50000'] = x[-50000:].std()\n    X_train.loc[segment, 'std_first_10000'] = x[:10000].std()\n    X_train.loc[segment, 'std_last_10000'] = x[-10000:].std()\n    \n    X_train.loc[segment, 'avg_first_50000'] = x[:50000].mean()\n    X_train.loc[segment, 'avg_last_50000'] = x[-50000:].mean()\n    X_train.loc[segment, 'avg_first_10000'] = x[:10000].mean()\n    X_train.loc[segment, 'avg_last_10000'] = x[-10000:].mean()\n    \n    X_train.loc[segment, 'min_first_50000'] = x[:50000].min()\n    X_train.loc[segment, 'min_last_50000'] = x[-50000:].min()\n    X_train.loc[segment, 'min_first_10000'] = x[:10000].min()\n    X_train.loc[segment, 'min_last_10000'] = x[-10000:].min()\n    \n    X_train.loc[segment, 'max_first_50000'] = x[:50000].max()\n    X_train.loc[segment, 'max_last_50000'] = x[-50000:].max()\n    X_train.loc[segment, 'max_first_10000'] = x[:10000].max()\n    X_train.loc[segment, 'max_last_10000'] = x[-10000:].max()\n    \n    X_train.loc[segment, 'max_to_min'] = x.max() / np.abs(x.min())\n    X_train.loc[segment, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n    X_train.loc[segment, 'count_big'] = len(x[np.abs(x) > 500])\n    X_train.loc[segment, 'sum'] = x.sum()\n    \n    # Quantiles\n    X_train.loc[segment, 'q95'] = np.quantile(x, 0.95)\n    X_train.loc[segment, 'q99'] = np.quantile(x, 0.99)\n    X_train.loc[segment, 'q05'] = np.quantile(x, 0.05)\n    X_train.loc[segment, 'q01'] = np.quantile(x, 0.01)\n    \n    X_train.loc[segment, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n    X_train.loc[segment, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n    X_train.loc[segment, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n    X_train.loc[segment, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    # Rolling statistical measures\n    for windows in [10, 100, 1000]:\n        x_roll_std = x.rolling(windows).std().dropna().values\n        x_roll_mean = x.rolling(windows).mean().dropna().values\n        \n        X_train.loc[segment, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X_train.loc[segment, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X_train.loc[segment, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X_train.loc[segment, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X_train.loc[segment, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X_train.loc[segment, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X_train.loc[segment, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X_train.loc[segment, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X_train.loc[segment, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X_train.loc[segment, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X_train.loc[segment, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X_train.loc[segment, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X_train.loc[segment, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X_train.loc[segment, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X_train.loc[segment, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X_train.loc[segment, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X_train.loc[segment, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X_train.loc[segment, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X_train.loc[segment, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X_train.loc[segment, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","metadata":{"_uuid":"d0521e2f-3209-4203-bbdd-114cde9ad927","_cell_guid":"557ca7b5-7b5a-41b6-97f2-a0540a6860d8","execution":{"iopub.status.busy":"2023-04-23T19:05:37.745752Z","iopub.execute_input":"2023-04-23T19:05:37.746092Z","iopub.status.idle":"2023-04-23T19:07:36.757044Z","shell.execute_reply.started":"2023-04-23T19:05:37.746060Z","shell.execute_reply":"2023-04-23T19:07:36.755717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Printing shapes and X_train\nprint(X_train.shape)\nprint(y_train.shape)\nprint(X_train)","metadata":{"_uuid":"fa64cb94-35e0-4e0c-841a-04205c9891b0","_cell_guid":"85b882a8-8ab8-44ca-9569-dfedb4e949b6","execution":{"iopub.status.busy":"2023-04-23T19:07:36.759004Z","iopub.execute_input":"2023-04-23T19:07:36.759485Z","iopub.status.idle":"2023-04-23T19:07:36.792823Z","shell.execute_reply.started":"2023-04-23T19:07:36.759437Z","shell.execute_reply":"2023-04-23T19:07:36.791875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Perform PCA","metadata":{}},{"cell_type":"code","source":"#standard normalize the data\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X_train)\nX_train = scaler.transform(X_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:36.794216Z","iopub.execute_input":"2023-04-23T19:07:36.794517Z","iopub.status.idle":"2023-04-23T19:07:36.902941Z","shell.execute_reply.started":"2023-04-23T19:07:36.794489Z","shell.execute_reply":"2023-04-23T19:07:36.901967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train pca\nfrom sklearn.decomposition import PCA\npca = PCA(n_components=20) # doing pca and keeping only n_components, shows the first 20 components\npca = pca.fit(X_train) # the correct dimension of X for sklearn is P*N (samples*features)\nX_pca_skl = pca.transform(X_train)\n\n# perform pca on features\nplt.bar(range(0,20), pca.explained_variance_ratio_, label=\"individual var\");\nplt.step(range(0,20), np.cumsum(pca.explained_variance_ratio_),'r', label=\"cumulative var\");\nplt.xlabel('Principal component index'); plt.ylabel('explained variance ratio %');\nplt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:36.906768Z","iopub.execute_input":"2023-04-23T19:07:36.907103Z","iopub.status.idle":"2023-04-23T19:07:37.421290Z","shell.execute_reply.started":"2023-04-23T19:07:36.907072Z","shell.execute_reply":"2023-04-23T19:07:37.420065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plots the best feature and the time to failure data\nfig, ax1 = plt.subplots(figsize=(16, 4))\nplt.title(\"Trends of best feature and time_to_failure. 2% of data (sampled)\")\nplt.plot(X_pca_skl[:,0], color='b')\nax1.set_ylabel('best feature', color='b')\nplt.legend(['best feature'])\nax2 = ax1.twinx()\nplt.plot(y_train, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:58:34.161394Z","iopub.execute_input":"2023-04-23T19:58:34.161877Z","iopub.status.idle":"2023-04-23T19:58:34.560916Z","shell.execute_reply.started":"2023-04-23T19:58:34.161836Z","shell.execute_reply":"2023-04-23T19:58:34.559749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.3 Data Preparation\n\nThe features which are extracted from the data are already numerical and there aren't any missing values. So there is no need to transfrom the data. First the data will be splitted and then normalized. ","metadata":{"_uuid":"2069db4f-6a52-4d88-a8f8-73e77a86a551","_cell_guid":"1a04ef55-9530-454e-a9d0-b6187ca8d79f","trusted":true}},{"cell_type":"markdown","source":"### 2.3.1 Train-test split\nThe data will be splitted into a training and test set. For the X_train, the three best features of the feature extraction are used. The size of the test set is 25 % of the total data set. The random state is set to 102. This is done to get the same random state each time that the code is run to get consistent results for testing. ","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX_train = X_pca_skl[:,0:3] # Three best features\n\n# Randomize the set and split it into a training and test set\nx_train, x_test, Y_train, Y_test = train_test_split(X_train, y_train, test_size=0.25,shuffle=False, random_state=102)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:37.863829Z","iopub.execute_input":"2023-04-23T19:07:37.864170Z","iopub.status.idle":"2023-04-23T19:07:37.871342Z","shell.execute_reply.started":"2023-04-23T19:07:37.864138Z","shell.execute_reply":"2023-04-23T19:07:37.870074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#standard normalize the data\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(x_train)\nx_train_norm = scaler.transform(x_train)\nx_test_norm= scaler.transform(x_test)","metadata":{"_uuid":"1b2f9833-88d2-4fce-9dea-804b31259ff9","_cell_guid":"b241c027-94b7-4edd-8eca-378ee927901e","execution":{"iopub.status.busy":"2023-04-23T19:07:37.872731Z","iopub.execute_input":"2023-04-23T19:07:37.873588Z","iopub.status.idle":"2023-04-23T19:07:37.880282Z","shell.execute_reply.started":"2023-04-23T19:07:37.873555Z","shell.execute_reply":"2023-04-23T19:07:37.879319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Training and Results","metadata":{}},{"cell_type":"markdown","source":"Two different models will be used for the training and testing of the data. One linear model and one non-linear model. The linear model that is used is SVM and the non-linear model is the random forest. The scores and the MSE of the models are also shown to show how well the model resembles the actual data. The linear model works better on the test data and the non-linear model works better on the training data.\n    \nThe linear model is a fast method to get a reasonable score, but it has a lower accuracy than non-linear models. \nThe non-linear model is accurate and has a higher score, but it is slower and it overfits the training data, which results in a lower test data score. ","metadata":{"_uuid":"f4f38413-6a4c-4c15-8eaa-ba351449306d","_cell_guid":"ed5e8ac3-008a-4384-be59-7eefdf4c2776","trusted":true}},{"cell_type":"markdown","source":"## 3.1 Training and testing","metadata":{}},{"cell_type":"markdown","source":"### Linear model","metadata":{}},{"cell_type":"code","source":"# Linear model\n# Support Vector Regression model\nfrom sklearn.svm import SVR\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_absolute_error\n\nmodel = SVR()\nmodel.fit(x_train_norm, Y_train.values.flatten())\n\n# Predict\nY_train_pred = model.predict(x_train_norm)\nY_test_pred = model.predict(x_test_norm)\n\n# accuracy\nscore1=model.score(x_train_norm, Y_train)\nscore2=model.score(x_test_norm, Y_test)\nprint(\"Training\\nR2:\",score1)\nprint(\"MSE:\",mean_squared_error(Y_train, Y_train_pred))\nprint(\"MAE:\",mean_absolute_error(Y_train, Y_train_pred))\n\nprint(\"\\nTest\\nR2:\",score2)\nprint(\"MSE:\",mean_squared_error(Y_test, Y_test_pred))\nprint(\"MAE:\",mean_absolute_error(Y_test, Y_test_pred))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T20:01:12.951275Z","iopub.execute_input":"2023-04-23T20:01:12.952335Z","iopub.status.idle":"2023-04-23T20:01:14.695314Z","shell.execute_reply.started":"2023-04-23T20:01:12.952295Z","shell.execute_reply":"2023-04-23T20:01:14.694155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction train data\nfig, ax1 = plt.subplots(figsize=(20, 4))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(Y_train_pred, color='b')\nax1.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure predicted'],)\nax2 = ax1.twinx()\nplt.plot(Y_train, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:39.613582Z","iopub.execute_input":"2023-04-23T19:07:39.614533Z","iopub.status.idle":"2023-04-23T19:07:40.067199Z","shell.execute_reply.started":"2023-04-23T19:07:39.614485Z","shell.execute_reply":"2023-04-23T19:07:40.065987Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction test data\nfig, ax1 = plt.subplots(figsize=(20, 4))\nx = np.linspace(0, 1049, 1049)\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(Y_test_pred, color='b')\nax1.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure predicted'],)\nax2 = ax1.twinx()\nplt.plot(x,Y_test, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.8))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:40.068786Z","iopub.execute_input":"2023-04-23T19:07:40.069236Z","iopub.status.idle":"2023-04-23T19:07:40.456548Z","shell.execute_reply.started":"2023-04-23T19:07:40.069191Z","shell.execute_reply":"2023-04-23T19:07:40.455316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Non-Linear model","metadata":{}},{"cell_type":"code","source":"# Non-Linear model\n# Random forest regressor model\nfrom sklearn.ensemble import RandomForestRegressor\n\nmodel = RandomForestRegressor(max_depth=5, random_state=102)\n\nmodel.fit(x_train_norm, Y_train.values.flatten())\n\n# Predict\nY_train_pred = model.predict(x_train_norm)\nY_test_pred = model.predict(x_test_norm)\n\n# accuracy\nscore1=model.score(x_train_norm, Y_train)\nscore2=model.score(x_test_norm, Y_test)\nprint(\"Training\\nR2:\",score1)\nprint(\"MSE:\",mean_squared_error(Y_train, Y_train_pred))\nprint(\"MAE:\",mean_absolute_error(Y_train, Y_train_pred))\n\nprint(\"\\nTest\\nR2:\",score2)\nprint(\"MSE:\",mean_squared_error(Y_test, Y_test_pred))\nprint(\"MAE:\",mean_absolute_error(Y_test, Y_test_pred))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T20:01:36.905628Z","iopub.execute_input":"2023-04-23T20:01:36.906091Z","iopub.status.idle":"2023-04-23T20:01:37.402066Z","shell.execute_reply.started":"2023-04-23T20:01:36.906051Z","shell.execute_reply":"2023-04-23T20:01:37.401175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction train data\nfig, ax1 = plt.subplots(figsize=(20, 4))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(Y_train_pred, color='b')\nax1.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure predicted'],)\nax2 = ax1.twinx()\nplt.plot(Y_train, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:41.277453Z","iopub.execute_input":"2023-04-23T19:07:41.278455Z","iopub.status.idle":"2023-04-23T19:07:41.720955Z","shell.execute_reply.started":"2023-04-23T19:07:41.278420Z","shell.execute_reply":"2023-04-23T19:07:41.719752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Prediction test data\nfig, ax1 = plt.subplots(figsize=(20, 4))\nx = np.linspace(0, 1049, 1049)\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(Y_test_pred, color='b')\nax1.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure predicted'],)\nax2 = ax1.twinx()\nplt.plot(x,Y_test, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.8))","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:41.722728Z","iopub.execute_input":"2023-04-23T19:07:41.723454Z","iopub.status.idle":"2023-04-23T19:07:42.246702Z","shell.execute_reply.started":"2023-04-23T19:07:41.723409Z","shell.execute_reply":"2023-04-23T19:07:42.245536Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2 Submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id') # submission csv\nX_test = pd.DataFrame(dtype=np.float64, index=submission.index)\n    \nfor seg_id in X_test.index:\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv') # read every test data and extract features\n    \n    x = pd.Series(seg['acoustic_data'].values)\n    \n    X_test.loc[seg_id, 'ave'] = x.mean()\n    X_test.loc[seg_id, 'std'] = x.std()\n    X_test.loc[seg_id, 'max'] = x.max()\n    X_test.loc[seg_id, 'min'] = x.min()\n    \n    X_test.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(x))\n    X_test.loc[seg_id, 'abs_max'] = np.abs(x).max()\n    X_test.loc[seg_id, 'abs_min'] = np.abs(x).min()\n    \n    X_test.loc[seg_id, 'std_first_50000'] = x[:50000].std()\n    X_test.loc[seg_id, 'std_last_50000'] = x[-50000:].std()\n    X_test.loc[seg_id, 'std_first_10000'] = x[:10000].std()\n    X_test.loc[seg_id, 'std_last_10000'] = x[-10000:].std()\n    \n    X_test.loc[seg_id, 'avg_first_50000'] = x[:50000].mean()\n    X_test.loc[seg_id, 'avg_last_50000'] = x[-50000:].mean()\n    X_test.loc[seg_id, 'avg_first_10000'] = x[:10000].mean()\n    X_test.loc[seg_id, 'avg_last_10000'] = x[-10000:].mean()\n\n    X_test.loc[seg_id, 'min_first_50000'] = x[:50000].min()\n    X_test.loc[seg_id, 'min_last_50000'] = x[-50000:].min()\n    X_test.loc[seg_id, 'min_first_10000'] = x[:10000].min()\n    X_test.loc[seg_id, 'min_last_10000'] = x[-10000:].min()\n    \n    X_test.loc[seg_id, 'max_first_50000'] = x[:50000].max()\n    X_test.loc[seg_id, 'max_last_50000'] = x[-50000:].max()\n    X_test.loc[seg_id, 'max_first_10000'] = x[:10000].max()\n    X_test.loc[seg_id, 'max_last_10000'] = x[-10000:].max()\n    \n    X_test.loc[seg_id, 'max_to_min'] = x.max() / np.abs(x.min())\n    X_test.loc[seg_id, 'max_to_min_diff'] = x.max() - np.abs(x.min())\n    X_test.loc[seg_id, 'count_big'] = len(x[np.abs(x) > 500])\n    X_test.loc[seg_id, 'sum'] = x.sum()\n\n    X_test.loc[seg_id, 'q95'] = np.quantile(x, 0.95)\n    X_test.loc[seg_id, 'q99'] = np.quantile(x, 0.99)\n    X_test.loc[seg_id, 'q05'] = np.quantile(x, 0.05)\n    X_test.loc[seg_id, 'q01'] = np.quantile(x, 0.01)\n    \n    X_test.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(x), 0.95)\n    X_test.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(x), 0.99)\n    X_test.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(x), 0.05)\n    X_test.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    for windows in [10, 100, 1000]:\n        x_roll_std = x.rolling(windows).std().dropna().values\n        x_roll_mean = x.rolling(windows).mean().dropna().values\n        \n        X_test.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n        X_test.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n        X_test.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n        X_test.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n        X_test.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n        X_test.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n        X_test.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n        X_test.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n        X_test.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n        X_test.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n        \n        X_test.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n        X_test.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n        X_test.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n        X_test.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n        X_test.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n        X_test.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n        X_test.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n        X_test.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n        X_test.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n        X_test.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max()","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:07:42.248530Z","iopub.execute_input":"2023-04-23T19:07:42.249240Z","iopub.status.idle":"2023-04-23T19:09:58.830261Z","shell.execute_reply.started":"2023-04-23T19:07:42.249185Z","shell.execute_reply":"2023-04-23T19:09:58.829123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(X_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:09:58.831941Z","iopub.execute_input":"2023-04-23T19:09:58.832410Z","iopub.status.idle":"2023-04-23T19:09:58.859922Z","shell.execute_reply.started":"2023-04-23T19:09:58.832363Z","shell.execute_reply":"2023-04-23T19:09:58.858750Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#standard normalize the data\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X_test)\nX_test_norm = scaler.transform(X_test)\n\n# train pca\nfrom sklearn.decomposition import PCA\npca = PCA(n_components=20) # doing pca and keeping only n_components, shows the first 20 components\npca = pca.fit(X_test_norm) # the correct dimension of X for sklearn is P*N (samples*features)\nX_pca_skl = pca.transform(X_test_norm)\nX_test_pca = X_pca_skl[:,0:3]","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:09:58.861407Z","iopub.execute_input":"2023-04-23T19:09:58.861736Z","iopub.status.idle":"2023-04-23T19:09:58.890938Z","shell.execute_reply.started":"2023-04-23T19:09:58.861703Z","shell.execute_reply":"2023-04-23T19:09:58.889387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submit to the submission file\nscaler = preprocessing.StandardScaler().fit(X_test_pca)\nX_test_pca_norm = scaler.transform(X_test_pca)\nsubmission['time_to_failure'] = model.predict(X_test_pca_norm) \nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:09:58.892740Z","iopub.execute_input":"2023-04-23T19:09:58.893307Z","iopub.status.idle":"2023-04-23T19:09:58.955616Z","shell.execute_reply.started":"2023-04-23T19:09:58.893259Z","shell.execute_reply":"2023-04-23T19:09:58.954237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(submission)","metadata":{"execution":{"iopub.status.busy":"2023-04-23T19:09:58.957349Z","iopub.execute_input":"2023-04-23T19:09:58.958070Z","iopub.status.idle":"2023-04-23T19:09:58.967665Z","shell.execute_reply.started":"2023-04-23T19:09:58.958002Z","shell.execute_reply":"2023-04-23T19:09:58.966198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Discussion and Conclusion\nThe best feature from the PCA was a good feature, but it was still not the best, because it still had random unwanted peaks. So more features could lead to a better scores. The training-test split could also be optimized, so choosing a different split and look if the score improves. Or even using no split and use the whole data for training. \n\nThe models worked quite well, but it could be better to explore different models. The random tree model is a good model, but it overfits the training data, which will lead to a worse test score. Next time we could use K-Folds Cross Validation to improve the accuracy of the scores. In our case, the random tree with depth 5 worked the best. We tested it with the test data and this model gave the best score. \n\nIn our case there were no catogorial features, so cleaning the data was not needed. But if it was needed then encoding the data could be a solution. So that means that each category will be assigned a certain integer, which makes the feature numerical. Missing data can be solved by interpolating the data. So for example the mean or the previous value of the feature can be filled into the missing rows. \n\nThe Mean Squared Error and Mean Absolute Error are used for the cost function. The mean absolute error should give a more accurate result of the performace than the mean squared error, because mean absolute error works better when the data has noise. Also the score is calculated, which is the R^2 score in this case. This indicates how well the predicted data represents the real data. ","metadata":{"_uuid":"e69c2c27-58df-4ddc-ad71-951affad4d06","_cell_guid":"e407ef1b-1c60-48f8-b1cd-ba51280b2a93","trusted":true}},{"cell_type":"markdown","source":"# 5. References\nEarthquakes FE. More features and samples by ANDREW LUKYANENKO, https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\n\nBasic Feature Benchmark by INVERSION, https://www.kaggle.com/code/inversion/basic-feature-benchmark/notebook\n\nRolling + Quantiles by PATRICK YAM, https://www.kaggle.com/code/wimwim/rolling-quantiles/notebook","metadata":{}}]}