{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":11000,"databundleVersionId":875412,"sourceType":"competition"}],"dockerImageVersionId":30458,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import KFold, StratifiedKFold, RepeatedKFold, train_test_split, cross_val_score\nfrom sklearn.linear_model import LinearRegression","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:22:34.694933Z","iopub.execute_input":"2024-04-21T16:22:34.695481Z","iopub.status.idle":"2024-04-21T16:22:34.705051Z","shell.execute_reply.started":"2024-04-21T16:22:34.695427Z","shell.execute_reply":"2024-04-21T16:22:34.703005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Introduction\n","metadata":{"id":"3NnAgB-y0knj"}},{"cell_type":"markdown","source":"<font color=#6698FF>\n    Olivierkd <br>\n    Wesley Heijenberg <br>\n    Score: 2.68631 \n    \n    \n    \n","metadata":{"id":"DPJV5e8sDVnT"}},{"cell_type":"markdown","source":"## 2. Data\n","metadata":{"id":"PB7FLdQn-dmK"}},{"cell_type":"markdown","source":"### 2.1 Dataset\n<font color=#6698FF>In this section, we load and explore the dataset. <br>\n    We split the data into 70% training and 30% testing data. These are typical values to make sure there is enough training data, but not too little test data to check performance with.","metadata":{"id":"vfGikOxUAxJB"}},{"cell_type":"code","source":"data = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure':np.float64})\n\ntrain, test = train_test_split(data, test_size=0.3, shuffle=False) # Split data into 70% training and 30% test\ntrain.head(10)\ntest.head(10)\n\ndel data","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:22:37.642244Z","iopub.execute_input":"2024-04-21T16:22:37.643590Z","iopub.status.idle":"2024-04-21T16:26:32.558688Z","shell.execute_reply.started":"2024-04-21T16:22:37.643493Z","shell.execute_reply":"2024-04-21T16:26:32.556683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(train.shape)\nprint(test.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:27:54.876425Z","iopub.execute_input":"2024-04-21T16:27:54.878114Z","iopub.status.idle":"2024-04-21T16:27:54.886002Z","shell.execute_reply.started":"2024-04-21T16:27:54.878054Z","shell.execute_reply":"2024-04-21T16:27:54.884167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\n\n<font color=#6698FF> Explore the features and target variables of the dataset. Think about making some scatter plots, box plots, histograms or printing the data, but feel free to choose any method that suits you.\nWhat do you think is the right performance\nmetric to use for this dataset? Clearly explain which performance metric you\nchoose and why.\nAlgorithmic bias can be a real problem in Machine Learning.  Explain what you believe. <br>\n\nThe data has a lot of spikes in it which would indicate a failure that has happened or is about to happen. Using metrics such as minimum, maximum, mean but also other statistical metrics such as kurtosis and skew may be useful metrics.\n</font>","metadata":{"id":"QN5jS2_p_Y9N"}},{"cell_type":"code","source":"#https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples\nacoustic_data_ = train['acoustic_data'].values[::75]  \ntime_to_failure_ = train['time_to_failure'].values[::75]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"Trends of acoustic_data and time_to_failure. 2% of data (sampled)\")\nplt.plot(acoustic_data_, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(time_to_failure_, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:27:56.509011Z","iopub.execute_input":"2024-04-21T16:27:56.509545Z","iopub.status.idle":"2024-04-21T16:28:04.698158Z","shell.execute_reply.started":"2024-04-21T16:27:56.509482Z","shell.execute_reply":"2024-04-21T16:28:04.696862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The green line shows the the time to failure and the acoustic data is in blue. The acoustic data shows several spikes after that is some time after failure (negative time to failure)","metadata":{}},{"cell_type":"code","source":"acoustic_data_ = train['acoustic_data'].values[:6291455]  \ntime_to_failure_ = train['time_to_failure'].values[:6291455]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title(\"More detailed acoustic data (1% of data)\")\nplt.plot(acoustic_data_, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(time_to_failure_, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel acoustic_data_\ndel time_to_failure_","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:28:10.674280Z","iopub.execute_input":"2024-04-21T16:28:10.674740Z","iopub.status.idle":"2024-04-21T16:28:18.242416Z","shell.execute_reply.started":"2024-04-21T16:28:10.674691Z","shell.execute_reply":"2024-04-21T16:28:18.240976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is a more detailed zoomed in example of the dataset, which shows a large spike.","metadata":{}},{"cell_type":"markdown","source":"### 2.3 Data Preparation\n\n<font color=#6698FF>  This dataset hasn’t been cleaned yet. Meaning that some attributes (features) are in numerical format and some are in categorial format. Moreover, there are missing values as well. However, all Scikit-learn’s implementations of these algorithms expect numerical features. Check for all features if they are in categorial and use a method to transform them to numerical values. For the numerical data, handle the missing data and normalize the data. \nNote that you are only allowed to use training data for preprocessing but you then need to perform similar changes on test data too.\nYou can use [pipelining](https://scikit-learn.org/stable/modules/generated/sklearn.pipeline.Pipeline.html) to help with the preprocessing.</font> <br>\nThe features used are: mean, std, min, max, skew, kurtosis, but also a few fft features such as imaginary and real mean, min and max.","metadata":{"id":"4A5GsVzgATRu"}},{"cell_type":"code","source":"rows = 150000 # Amount of rows per segment\ntrain_segments = int(np.floor(train.shape[0] / rows)) # Amount of segments in dataset\nprint(\"Amount of segments: \", train_segments)\n\nx_train = pd.DataFrame(index=range(train_segments), dtype=np.float64, columns=['mean', 'std', 'min', 'max', 'skew', 'kurtosis', 'Imean', 'Rmean', 'Imin', 'Rmin', 'Imax', 'Rmax', 'max_to_min'])\ny_train = pd.DataFrame(index=range(train_segments), dtype=np.float64, columns=['time_to_failure'])","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:28:21.136122Z","iopub.execute_input":"2024-04-21T16:28:21.136694Z","iopub.status.idle":"2024-04-21T16:28:21.153448Z","shell.execute_reply.started":"2024-04-21T16:28:21.136636Z","shell.execute_reply":"2024-04-21T16:28:21.151739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def extract_features(seg_id, seg, X):\n    values = pd.Series(seg['acoustic_data'].values)\n    values_fft = np.fft.fft(values)\n    values_real = np.real(values_fft)\n    values_imag = np.imag(values_fft)\n    \n    X.loc[seg_id, 'mean'] = values.mean()\n    X.loc[seg_id, 'std'] = values.std()\n    X.loc[seg_id, 'min'] = values.min()\n    X.loc[seg_id, 'max'] = values.max()\n    X.loc[seg_id, 'skew'] = values.skew()\n    X.loc[seg_id, 'kurtosis'] = values.kurt()\n    \n    X.loc[seg_id, 'Imean'] = values_imag.mean()\n    X.loc[seg_id, 'Rmean'] = values_real.mean()\n    \n    X.loc[seg_id, 'Imin'] = values_imag.min()\n    X.loc[seg_id, 'Rmin'] = values_real.min()\n    \n    X.loc[seg_id, 'Imax'] = values_imag.max()\n    X.loc[seg_id, 'Rmax'] = values_real.max()\n    \n    X.loc[seg_id, 'max_to_min'] = values.max() - values.min()","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:28:22.659117Z","iopub.execute_input":"2024-04-21T16:28:22.659683Z","iopub.status.idle":"2024-04-21T16:28:22.672153Z","shell.execute_reply.started":"2024-04-21T16:28:22.659628Z","shell.execute_reply":"2024-04-21T16:28:22.670340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for seg_id in range(train_segments):\n    seg = train.iloc[seg_id * rows : seg_id * rows + rows] # Select segment data points\n    extract_features(seg_id, seg, x_train) # Extract the features for this segment\n    y_train.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1] # Copy corresponding output (time_to_failure)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:28:32.070948Z","iopub.execute_input":"2024-04-21T16:28:32.071438Z","iopub.status.idle":"2024-04-21T16:29:03.678597Z","shell.execute_reply.started":"2024-04-21T16:28:32.071398Z","shell.execute_reply":"2024-04-21T16:29:03.677150Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Output shape: \", x_train.shape)\nx_train.head(10) # Show first 10 rows","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:29:11.562559Z","iopub.execute_input":"2024-04-21T16:29:11.563713Z","iopub.status.idle":"2024-04-21T16:29:11.605601Z","shell.execute_reply.started":"2024-04-21T16:29:11.563649Z","shell.execute_reply":"2024-04-21T16:29:11.604064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_segments = int(np.floor(test.shape[0] / rows))\nprint(\"Amount of test segments: \", test_segments)\n\nx_test = pd.DataFrame(index=range(test_segments), dtype=np.float64, columns=['mean', 'std', 'min', 'max', 'skew', 'kurtosis', 'Imean', 'Rmean', 'Imin', 'Rmin', 'Imax', 'Rmax', 'max_to_min'])\ny_test = pd.DataFrame(index=range(test_segments), dtype=np.float64, columns=['time_to_failure'])\n\nfor seg_id in range(test_segments):\n    seg = test.iloc[seg_id * rows : seg_id * rows + rows]\n    extract_features(seg_id, seg, x_test)\n    y_test.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1]\n\nprint(\"Test output shape: \", x_test.shape)\nx_test.head(10)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:29:14.987621Z","iopub.execute_input":"2024-04-21T16:29:14.989235Z","iopub.status.idle":"2024-04-21T16:29:28.457865Z","shell.execute_reply.started":"2024-04-21T16:29:14.989159Z","shell.execute_reply":"2024-04-21T16:29:28.456452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Scale data\ntrain_scaler = StandardScaler().fit(x_train)\ntest_scaler = StandardScaler().fit(x_test)\n\nx_train_scaled = train_scaler.transform(x_train)\nx_test_scaled = test_scaler.transform(x_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:29:50.130777Z","iopub.execute_input":"2024-04-21T16:29:50.131471Z","iopub.status.idle":"2024-04-21T16:29:50.151608Z","shell.execute_reply.started":"2024-04-21T16:29:50.131382Z","shell.execute_reply":"2024-04-21T16:29:50.149816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del train\ndel test","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:29:56.339314Z","iopub.execute_input":"2024-04-21T16:29:56.340623Z","iopub.status.idle":"2024-04-21T16:29:56.464097Z","shell.execute_reply.started":"2024-04-21T16:29:56.340538Z","shell.execute_reply":"2024-04-21T16:29:56.461282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 3. Training and Results\n<font color=#6698FF> Briefly introduce algorithms you choose. \nPresent your balanced accuracies for both test and training data for all classifiers. Analyse the performance on test and training in terms of bias and variance. Give one advantage and one drawback of the method you use. <br>\n    \nFor the linear model we used linear regression and for the non-linear model we chose random forest regression. Linear regression is simple and fast, however its linearity means it is not very useful for non-linear datasets. Random forest regression is non-linear which is better for non-linear datasets, however this model is more likely to overfit.\n","metadata":{"id":"IUcc2ue7yrOv"}},{"cell_type":"code","source":"def train_model(x_train, y_train, x_test, y_test, model):\n    model.fit(x_train, y_train.values.flatten())\n    \n    train_pred = model.predict(x_train)\n    test_pred = model.predict(x_test)\n    \n    train_score = mean_squared_error(y_train, train_pred) # Calculate MSE for training data\n    test_score = mean_squared_error(y_test, test_pred) # Calculate MSE for test data\n    \n    print(\"Train MSE: \", train_score)\n    print(\"Test MSE: \", test_score)\n    \n    plt.scatter(y_train.values.flatten(), train_pred, label=\"Train\")\n    plt.scatter(y_test.values.flatten(), test_pred, label=\"Test\")\n    plt.plot([(0, 0), (16, 16)], [(0, 0), (16, 16)], color='g')\n    plt.xlim([0, 16])\n    plt.ylim([0, 16])\n    plt.xlabel(\"Expected\")\n    plt.ylabel(\"Actual\")\n    plt.legend()","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:30:02.220286Z","iopub.execute_input":"2024-04-21T16:30:02.221956Z","iopub.status.idle":"2024-04-21T16:30:02.234250Z","shell.execute_reply.started":"2024-04-21T16:30:02.221856Z","shell.execute_reply":"2024-04-21T16:30:02.232672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_model(x_train_scaled, y_train, x_test_scaled, y_test, LinearRegression())","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:30:05.478989Z","iopub.execute_input":"2024-04-21T16:30:05.479484Z","iopub.status.idle":"2024-04-21T16:30:05.936383Z","shell.execute_reply.started":"2024-04-21T16:30:05.479441Z","shell.execute_reply":"2024-04-21T16:30:05.934113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The linear model performs extremely poorly, the results on training data is very flat being only accurate between 4 and 8. While the test data is so far off that it is not in the limits of the plot.","metadata":{}},{"cell_type":"code","source":"train_model(x_train_scaled, y_train, x_test_scaled, y_test, RandomForestRegressor())","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:30:09.286145Z","iopub.execute_input":"2024-04-21T16:30:09.286877Z","iopub.status.idle":"2024-04-21T16:30:12.495367Z","shell.execute_reply.started":"2024-04-21T16:30:09.286823Z","shell.execute_reply":"2024-04-21T16:30:12.493634Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The random forest regressor performs better, but does show overfitting. MSE on training data is decent at around 1.16, while the test data is rather poor at 8.96. Both are however still an improvement over the linear regression. Changing the random forest max depth in order to reduce unfortunately did not work.","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\nsubmission_x_test = pd.DataFrame(columns=x_train.columns, dtype=np.float64, index=submission.index)\n\nfor seg_id in submission_x_test.index:\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')\n    extract_features(seg_id, seg, submission_x_test)\n\nsubmission_x_test_scaled = train_scaler.transform(submission_x_test)\n\nmodel = RandomForestRegressor()\nmodel.fit(x_train_scaled, y_train.values.flatten())\nsubmission['time_to_failure'] = model.predict(submission_x_test_scaled)\nsubmission.to_csv('submission.csv', index=True)\nprint(submission)","metadata":{"execution":{"iopub.status.busy":"2024-04-21T16:30:18.597985Z","iopub.execute_input":"2024-04-21T16:30:18.598507Z","iopub.status.idle":"2024-04-21T16:31:49.672367Z","shell.execute_reply.started":"2024-04-21T16:30:18.598463Z","shell.execute_reply":"2024-04-21T16:31:49.670970Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion\n\n<font color=#6698FF>Discuss all the choices you made during the process and your final conclusions.  Highlight the strong points of your approach, discuss its shortcomings  and suggest some future approaches that may improve it. Please be self critical here. The assignment is not about achieving a state of the art performance, but about showing what you have learned the concepts during the course.</font> <br>\n\nThe chosen features are rather simple and could be the reason why the model did not perform very well. In the future it may be a good idea to add more features and then use PCA to reduce features and select the best ones. Using a different kind of model could also have made a difference. This assignment did teach us more about data exploration and preparation such as the graphing and segmentation.\n<br><br>\n\nIn the end our best attempt had an MSE on training data of 1.1679571234609787, while the MSE on test data was 8.967878453126216. This does show some overfitting which is further supported by the scatterplot","metadata":{"id":"yekgnxu7y0gE"}},{"cell_type":"markdown","source":"## 5. References","metadata":{"id":"wn7PBEmWESIi"}},{"cell_type":"markdown","source":"Gabriel Preda: https://www.kaggle.com/code/gpreda/lanl-earthquake-data-exploration-and-baseline (feature extraction and data exploration graphing) <br>\nInversion: https://www.kaggle.com/code/inversion/basic-feature-benchmark/notebook (results graphing) <br>\nJessteijn: https://www.kaggle.com/code/jessteijn/ml23-week8 (submission file code) <br>\nCheatsheet: https://colab.research.google.com/drive/12h-QBlsaWXkjGIRoJXfF4yqi1elnX9qn?usp=sharing#scrollTo=3p70JoCSyFK6 (Some functions such as train test split and models)","metadata":{}}]}