{"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":"code","source":"#import libraries\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom os import listdir\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import train_test_split\nfrom sklearn import preprocessing\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.decomposition import PCA\nfrom sklearn.metrics import mean_squared_error","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:48:53.723548Z","iopub.execute_input":"2023-04-27T19:48:53.723959Z","iopub.status.idle":"2023-04-27T19:48:54.619389Z","shell.execute_reply.started":"2023-04-27T19:48:53.723926Z","shell.execute_reply":"2023-04-27T19:48:54.618104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Introduction\nName: Ian van Ingen\n\\\nUsername: ianvaningen\n\\\nPrivate Score: 3.87917\n\\\nPublic score: 3.35218\n\\\nPrint functions used for confirming data values / shapes have been commented.","metadata":{}},{"cell_type":"markdown","source":"## 2. Data","metadata":{}},{"cell_type":"markdown","source":"### 2.1 Dataset\n\nBefore being able to make our models we need to load, explore and process our data. To begin we will load the training recording and a test segment to look at their formatting. We will also plot a small segment of the training recording to get a better understanding of our data.\n\\\nThe training recording is a Pandas DataFrame with 62914580 samples. Each sample has an acoustic_data value and a time_to_failure value. Our test segments have the same data type but have only 150000 samples with only acoustic_data values. From this information we can deduce that acoustic_data is our input and time_to_failure is our output.\n\\\nWe change the datatype of the training recording and the test segments to Numpy Arrays and split the training recording into 4194 segments of each 150000 elements. The output of each segment is the time_to_failure value at the end of the segment.\n\n","metadata":{}},{"cell_type":"code","source":"# Read training data\ndata = pd.read_csv(\"/kaggle/input/LANL-Earthquake-Prediction/train.csv\", dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:48:54.621451Z","iopub.execute_input":"2023-04-27T19:48:54.621781Z","iopub.status.idle":"2023-04-27T19:52:36.316678Z","shell.execute_reply.started":"2023-04-27T19:48:54.621748Z","shell.execute_reply":"2023-04-27T19:52:36.314556Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Print the training data, its shape and its data type \nprint(data)\nprint(data.shape)\nprint(type(data))","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:52:36.319376Z","iopub.execute_input":"2023-04-27T19:52:36.319932Z","iopub.status.idle":"2023-04-27T19:52:36.349453Z","shell.execute_reply.started":"2023-04-27T19:52:36.319870Z","shell.execute_reply":"2023-04-27T19:52:36.348209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# import test segments and take a look at their size\ntest_path = \"/kaggle/input/LANL-Earthquake-Prediction/test/\"\ntest_files = listdir(\"/kaggle/input/LANL-Earthquake-Prediction/test\")\n\nx_sub = np.zeros((len(test_files),150000))\nfor n in range(len(test_files)):\n    seg = pd.read_csv(test_path  + test_files[n])\n    x_sub[n] = seg.acoustic_data\n    #print(seg)\n    #print(type(seg))\n    #print(seg.shape)\nprint(x_sub.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:52:36.355304Z","iopub.execute_input":"2023-04-27T19:52:36.355758Z","iopub.status.idle":"2023-04-27T19:53:45.849897Z","shell.execute_reply.started":"2023-04-27T19:52:36.355724Z","shell.execute_reply":"2023-04-27T19:53:45.848495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot a small section of the data to have a visual representation\nn = 20 # percentage shown\nN = int((n/100)*len(data.acoustic_data))\nfig,ax=plt.subplots()\nax.plot(range(int(N/1000)+1),data.acoustic_data[0:N:1000], color = \"red\")\nax.set_ylabel(\"Acoustic data\")\nax2=ax.twinx()\nax2.plot(range(int(N/1000)+1), data.time_to_failure[0:N:1000], color = \"blue\")\nax2.set_ylabel(\"Time to Failure\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:53:45.851348Z","iopub.execute_input":"2023-04-27T19:53:45.851710Z","iopub.status.idle":"2023-04-27T19:53:46.298139Z","shell.execute_reply.started":"2023-04-27T19:53:45.851675Z","shell.execute_reply":"2023-04-27T19:53:46.297148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Change data type training data to numpy array\ndata = data.to_numpy()\nprint(type(data))\nprint(data.shape)\nprint(data)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:53:46.299904Z","iopub.execute_input":"2023-04-27T19:53:46.300278Z","iopub.status.idle":"2023-04-27T19:53:47.974156Z","shell.execute_reply.started":"2023-04-27T19:53:46.300241Z","shell.execute_reply":"2023-04-27T19:53:47.973146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make the input array consisting of 4194 inputs each consisting of 150000 samples and the output array\nsegment_size = 150000\nnum_segments = np.floor(629145480/segment_size)\nprint(num_segments)\ndata_short = data[0:int(num_segments*segment_size):1]\ndata = None\nprint(data_short.shape)\nx = np.zeros((int(num_segments), int(segment_size)))\ny = np.zeros((int(num_segments)))\nfor i in range(int(num_segments)):\n    if i == 0:\n        j = 0\n    else:\n        j += 150000\n    segment = data_short[j:j+150000:1]\n    x[i] = segment[:,0]\n    y[i] = segment[-1,1]\ndata_short = None\nsegment = None\n#print(x)\nprint(x.shape)\n#print(y)\n#print(y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:53:47.975715Z","iopub.execute_input":"2023-04-27T19:53:47.976357Z","iopub.status.idle":"2023-04-27T19:53:53.270696Z","shell.execute_reply.started":"2023-04-27T19:53:47.976320Z","shell.execute_reply":"2023-04-27T19:53:53.269323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\nWhen the training data has the preferred format, feature extraction can start. When looking at the previously plotted segment of the training recording we can see that there are large peaks before time_to_failure reaches its minimum. The amplitude and frequency of the peaks seems to be related to the time_to_failure values so we will focus on statistical features based on these characteristics. Inspiration for several features has been taken from other notebooks that are mentioned in section 5. ","metadata":{}},{"cell_type":"code","source":"# Features (based on things seen in plot and other notebooks)\nnum_features = 17\nX = np.zeros((int(num_features), int(num_segments)))\nX[0] = np.mean(x, axis=1)     # average acoustic data\nX[1] = np.max(x, axis=1)      # max value \nX[2] = np.max(abs(x), axis=1) # absolute max value\nX[3] = np.min(x, axis=1)      # min value\nX[4] = np.max(abs(x), axis=1) # absolute min value\nX[5] = np.std(x, axis=1)      # Standard deviation\nX[6] = np.quantile(x,0.01, axis=1)   # 1% quantile\nX[7] = np.quantile(x,0.05, axis=1)   # 5% quantile\nX[8] = np.quantile(x,0.25, axis=1)   # 25% quantile\nX[9] = np.quantile(x,0.5, axis=1)    # 50% quantile\nX[10] = np.quantile(x,0.75, axis=1)   # 75% quantile\nX[11] = np.quantile(x,0.95, axis=1)   # 95% quantile\nX[12] = np.quantile(x,0.99, axis=1)   # 99% quantile\nX[13] = np.quantile(x,0.999, axis=1)  # 99.9% quantile\nX[14] = X[5] - X[3]           # interquartile range\nX[15] = np.mean(np.diff(x, n=1, axis = 1), axis =1) # difference between upfollowing samples\nX[16] = np.mean(x[:,-25000:-1:1], axis =1) # mean value of the last 25000 samples\nX = X.T\n#print(X)\n#print(X.shape)\nx = None\n\n# Features test_segments\nX_sub = np.zeros((int(num_features), 2624))\nX_sub[0] = np.mean(x_sub, axis=1)     # average acoustic data\nX_sub[1] = np.max(x_sub, axis=1)      # max value \nX_sub[2] = np.max(abs(x_sub), axis=1) # absolute max value\nX_sub[3] = np.min(x_sub, axis=1)      # min value\nX_sub[4] = np.max(abs(x_sub), axis=1) # absolute min value\nX_sub[5] = np.std(x_sub, axis=1)      # Standard deviation\nX_sub[6] = np.quantile(x_sub,0.01, axis=1)   # 1% quantile\nX_sub[7] = np.quantile(x_sub,0.05, axis=1)   # 5% quantile\nX_sub[8] = np.quantile(x_sub,0.25, axis=1)   # 25% quantile\nX_sub[9] = np.quantile(x_sub,0.5, axis=1)    # 50% quantile\nX_sub[10] = np.quantile(x_sub,0.75, axis=1)   # 75% quantile\nX_sub[11] = np.quantile(x_sub,0.95, axis=1)   # 95% quantile\nX_sub[12] = np.quantile(x_sub,0.99, axis=1)   # 99% quantile\nX_sub[13] = np.quantile(x_sub,0.999, axis=1)  # 99.9% quantile\nX_sub[14] = X_sub[5] - X_sub[3]           # interquartile range\nX_sub[15] = np.mean(np.diff(x_sub, n=1, axis = 1), axis =1) # difference between upfollowing samples\nX_sub[16] = np.mean(x_sub[:,-25000:-1:1], axis =1) # mean value of the last 25000 samples\nX_sub = X_sub.T\n#print(X)\n#print(X.shape)\nx_sub = None\n","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:53:53.272793Z","iopub.execute_input":"2023-04-27T19:53:53.273444Z","iopub.status.idle":"2023-04-27T19:57:52.828612Z","shell.execute_reply.started":"2023-04-27T19:53:53.273391Z","shell.execute_reply":"2023-04-27T19:57:52.826900Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.3 Data Preparation\nTo be able to use these features we will have to split our data into a training and test set. We want a large training set but we have to make sure that our test set is big enough to test the accuracy of the model in a reliable way. Considering we have a sample size of 4194, we have a test set of around 800 samples when we use a train-test-ratio of 80 to 20. I consider this large enough for reliable testing.\n\\\nAfter splitting the data we normalise the features based on the train set so we can use these normalised features for our model training.","metadata":{}},{"cell_type":"code","source":"# Split data in test and train set\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,shuffle=True, random_state=66)\n#print(X_train)\n#print(X_test)\n#print(y_train)\n#print(y_test)\n#print(X_train.shape)\n#print(X_test.shape)\n#print(y_train.shape)\n#print(y_test.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:52.830683Z","iopub.execute_input":"2023-04-27T19:57:52.831094Z","iopub.status.idle":"2023-04-27T19:57:52.841942Z","shell.execute_reply.started":"2023-04-27T19:57:52.831055Z","shell.execute_reply":"2023-04-27T19:57:52.840521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Normalise data\nscaler = preprocessing.StandardScaler().fit(X_train)\nX_train = scaler.transform(X_train)\nX_test = scaler.transform(X_test)\nX_sub = scaler.transform(X_sub)\n#print(X_train_norm)\n#print(X_test_norm)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:52.846449Z","iopub.execute_input":"2023-04-27T19:57:52.846961Z","iopub.status.idle":"2023-04-27T19:57:52.860424Z","shell.execute_reply.started":"2023-04-27T19:57:52.846918Z","shell.execute_reply":"2023-04-27T19:57:52.858934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Training and results","metadata":{}},{"cell_type":"markdown","source":"### 3.1. Linear Regression\nThe first model we will train is the linear regression model. It is a simple, straight-forward model, but its drawback is that our predictions may have large deviations because the output is difficult to predict using a linear model.\n\\\nThe model has a class defined accuracy of 0.335 and MSE of 8.98 for the train set and values of 0.281 and 9.68 for the test set. \n\\\nWhen looking at the importance of the values we can see that some features are more relevant than others. To decrease the number of features we will apply PCA on our data. More than 96% of the variance can be explained with 7 PCA components so this is the number of components we will use.\n\\\nNow our results are an accuracy of 0.293 and MSE of 9.54 for the train set and 0.277 and 9.75 for the test set. This is a slightly worse score than the model not using PCA but this can be easily explained by the fact that we use less data to make our predictions.\n","metadata":{}},{"cell_type":"code","source":"# Make a linear regression model\n# Define model\nlin_model = LinearRegression()\n\n# Fit model\nlin_model.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:52.862325Z","iopub.execute_input":"2023-04-27T19:57:52.862731Z","iopub.status.idle":"2023-04-27T19:57:52.896258Z","shell.execute_reply.started":"2023-04-27T19:57:52.862695Z","shell.execute_reply":"2023-04-27T19:57:52.894564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Predict outputs of train and test data and get the accuracy defined by the class method and the MSE\nlin_predictions_train = lin_model.predict(X_train)\nlin_score_train = lin_model.score(X_train, y_train)\nlin_MSE_train = mean_squared_error(y_train, lin_predictions_train)\n\nlin_predictions_test = lin_model.predict(X_test)\nlin_score_test = lin_model.score(X_test, y_test)\nlin_MSE_test = mean_squared_error(y_test, lin_predictions_test)\n\nprint(lin_score_train)\nprint(lin_MSE_train)\n\nprint(lin_score_test)\nprint(lin_MSE_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:52.898923Z","iopub.execute_input":"2023-04-27T19:57:52.900060Z","iopub.status.idle":"2023-04-27T19:57:52.919972Z","shell.execute_reply.started":"2023-04-27T19:57:52.899979Z","shell.execute_reply":"2023-04-27T19:57:52.918151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Plot the importance of the features\nimportance = lin_model.coef_\nprint(importance)\nplt.bar([x for x in range(len(importance))], importance)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:52.922850Z","iopub.execute_input":"2023-04-27T19:57:52.924738Z","iopub.status.idle":"2023-04-27T19:57:53.148373Z","shell.execute_reply.started":"2023-04-27T19:57:52.924667Z","shell.execute_reply":"2023-04-27T19:57:53.147083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train pca\nfeatures = 7\npca = PCA(n_components=features)\npca = pca.fit(X_train)\n\n# perform pca on features\nX_train_pca=pca.transform(X_train);\nX_test_pca=pca.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:53.149817Z","iopub.execute_input":"2023-04-27T19:57:53.150204Z","iopub.status.idle":"2023-04-27T19:57:53.183771Z","shell.execute_reply.started":"2023-04-27T19:57:53.150169Z","shell.execute_reply":"2023-04-27T19:57:53.181964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# show variance \nplt.bar(range(0, features), pca.explained_variance_ratio_, label=\"individual var\");\nplt.step(range(0, features), np.cumsum (pca.explained_variance_ratio_), 'r', label=\"cumulative var\");\nplt.xlabel('Principal component index'); plt.ylabel('explained variance ratio %');\nplt.legend()\nplt.show()\n\n","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:53.186359Z","iopub.execute_input":"2023-04-27T19:57:53.187434Z","iopub.status.idle":"2023-04-27T19:57:53.433698Z","shell.execute_reply.started":"2023-04-27T19:57:53.187357Z","shell.execute_reply":"2023-04-27T19:57:53.432510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#training and testing on model\nlin_model_pca = LinearRegression()\n# fit the model\nlin_model_pca.fit(X_train_pca, y_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:53.435238Z","iopub.execute_input":"2023-04-27T19:57:53.435617Z","iopub.status.idle":"2023-04-27T19:57:53.447759Z","shell.execute_reply.started":"2023-04-27T19:57:53.435582Z","shell.execute_reply":"2023-04-27T19:57:53.446186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# results\nlin_pca_predictions_train = lin_model_pca.predict(X_train_pca)\nlin_pca_score_train = lin_model_pca.score(X_train_pca, y_train)\nlin_pca_MSE_train = mean_squared_error(y_train, lin_pca_predictions_train)\n\nlin_pca_predictions_test = lin_model_pca.predict(X_test_pca)\nlin_pca_score_test = lin_model_pca.score(X_test_pca, y_test)\nlin_pca_MSE_test = mean_squared_error(y_test, lin_pca_predictions_test)\n\nprint(lin_pca_score_train)\nprint(lin_pca_MSE_train)\n\nprint(lin_pca_score_test)\nprint(lin_pca_MSE_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:53.450641Z","iopub.execute_input":"2023-04-27T19:57:53.451242Z","iopub.status.idle":"2023-04-27T19:57:53.473607Z","shell.execute_reply.started":"2023-04-27T19:57:53.451190Z","shell.execute_reply":"2023-04-27T19:57:53.471612Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Plot the importance of the features\nimportance = lin_model_pca.coef_\nprint(importance)\nplt.bar([x for x in range(len(importance))], importance)\nplt.show()\n\n# Show the contribution of the original features to the PCA components\nfor i in range(features):\n  print(abs(pca.components_[i]))","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:53.475682Z","iopub.execute_input":"2023-04-27T19:57:53.476836Z","iopub.status.idle":"2023-04-27T19:57:53.688836Z","shell.execute_reply.started":"2023-04-27T19:57:53.476785Z","shell.execute_reply":"2023-04-27T19:57:53.687432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.2. Random Forest Regressor\nThe second model we will train is the Random Forest Regressor. This model shoulc be more accurate than the Linear Regression model because we are no longer limited to linear combinations of features but this model is more prone to overfitting. Overfitting can be prevented by choosing the right parameters for the model, primarily max_depth. The parameters have been chosen in a way that the train and test set have similar values for accuracy and MSE.\n\\\nThe score for the test data is the highest for a max_depth of 4. The accuracy and MSE of the train data are respectively 0.476 and 7.07 and the values for the test data are 0.416 and 7.87.\n","metadata":{}},{"cell_type":"code","source":"# Define and fit model\nRF_model = RandomForestRegressor(max_depth=4)\nRF_model.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:53.690887Z","iopub.execute_input":"2023-04-27T19:57:53.692280Z","iopub.status.idle":"2023-04-27T19:57:54.612306Z","shell.execute_reply.started":"2023-04-27T19:57:53.692227Z","shell.execute_reply":"2023-04-27T19:57:54.611055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# results\nRF_predictions_train = RF_model.predict(X_train)\nRF_score_train = RF_model.score(X_train, y_train)\nRF_MSE_train = mean_squared_error(y_train, RF_predictions_train)\n\nRF_predictions_test = RF_model.predict(X_test)\nRF_score_test = RF_model.score(X_test, y_test)\nRF_MSE_test = mean_squared_error(y_test, RF_predictions_test)\n\nprint(RF_score_train)\nprint(RF_MSE_train)\n\nprint(RF_score_test)\nprint(RF_MSE_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:54.614198Z","iopub.execute_input":"2023-04-27T19:57:54.615070Z","iopub.status.idle":"2023-04-27T19:57:54.704332Z","shell.execute_reply.started":"2023-04-27T19:57:54.615001Z","shell.execute_reply":"2023-04-27T19:57:54.703074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submission\nsubmission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\nsubmission['time_to_failure'] = RF_model.predict(X_sub) \nsubmission.to_csv('submission.csv')\nprint(submission)","metadata":{"execution":{"iopub.status.busy":"2023-04-27T19:57:54.705824Z","iopub.execute_input":"2023-04-27T19:57:54.706368Z","iopub.status.idle":"2023-04-27T19:57:54.766666Z","shell.execute_reply.started":"2023-04-27T19:57:54.706318Z","shell.execute_reply":"2023-04-27T19:57:54.764800Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion\n\nThe training recording was split into segments similar to the test segments. Feature extraction was done on these segments after which they were split into a training and test set and the features were normalised.\\\nLinear Regression was used first on the data after which PCA was implemented to decrease the number of features. However, this did not result in a higher accuracy. Random Forest Regressor gave a higher accuracy than both the Linear Regression Models and was used to predict the output of the test segments.\nThe accuracy could possibly be improved by adding additional features or using other, more relevant, regression models. The features that were chosen were rather simple and it is possible that more advanced features are more reliable. The same is true for the used models. The models were ones I already knew and there are models that are much more reliable and suitable for this assignment. In addition to this, it is very likely that the Random Forest Regressor would have had higher accuracy when different parameters had been chosen.\n","metadata":{}},{"cell_type":"markdown","source":"## 5. References\n\nThe following notebooks were used to get insight into the assignment and the data exploration and to get inspiration for feature extraction.\\\nhttps://www.kaggle.com/code/allunia/shaking-earth \\\nhttps://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples \\\nhttps://www.kaggle.com/code/andrekos/basic-feature-benchmark-with-quantiles/notebook","metadata":{}}]}