{"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":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Final Project\n\nAs in all machine learning problems you should complete the following steps:\n1. Load and explore the data (using plots and histograms )\n2. Clean/preprocess/transform the data if necessary (first performed on training data and next on test data)\n3. Train the machine learning model (in this assignment two models, one linear and one non-linear)\n4. Evaluate and optimise the model\n\nHowever, while performing the steps you should consider our main goals and research questions for doing this project and try to adress them. This can be either in each step or afterwards in discussion and conclusion section (your call!). Below are the research question we are interested in:\n\n* What are the necessary step to clean and prepare the dataset as we have categorial features and missing data.\n* What model/classifier provides the best result for this application.\n* What is the best choice for our cost function and performance metrics for this problem.\n\nPlease also note the following points during the assignment:\n\n* Use functions from open source libraries like sci-kit learn and keras and avoid using your own hand-written functions from previous assignments. You can also use our [cheatsheet](https://colab.research.google.com/drive/12h-QBlsaWXkjGIRoJXfF4yqi1elnX9qn?usp=sharing).\n* Feel free to contact us on Teams if you need more description.\n* This notebook is structured like a scientific paper. The text should provide a high-level overview of your approach. Please don't include any details about your code in the text but add them as comments in the code itself. Your code should be cleane and readable with enough comments.\n* There are some instructions and questions in each section, remove the highlighted text in blue and replace it with your explanations and answers.\n\nYou can delete this section before submission.","metadata":{"id":"GF4WsCxg-lL9"}},{"cell_type":"markdown","source":"## 1. Introduction\n","metadata":{"id":"3NnAgB-y0knj"}},{"cell_type":"markdown","source":"<font color=#6698FF> Bram van de Ruit, bramvanderuit <br>\n                       Teun Mertens, teunmertens <br>\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.\n    ","metadata":{"id":"vfGikOxUAxJB"}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.svm import NuSVR\nfrom sklearn import svm\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.metrics import mean_squared_error\nfrom scipy import stats\nimport numpy as np\nimport pandas as pd\nimport os","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:00:12.769265Z","iopub.execute_input":"2024-04-20T10:00:12.770711Z","iopub.status.idle":"2024-04-20T10:00:12.778729Z","shell.execute_reply.started":"2024-04-20T10:00:12.770655Z","shell.execute_reply":"2024-04-20T10:00:12.777301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = pd.read_csv('../input/LANL-Earthquake-Prediction/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:00:17.830403Z","iopub.execute_input":"2024-04-20T10:00:17.830893Z","iopub.status.idle":"2024-04-20T10:03:32.004641Z","shell.execute_reply.started":"2024-04-20T10:00:17.830849Z","shell.execute_reply":"2024-04-20T10:03:32.002976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_change_rate(x):\n    change = (np.diff(x) / x[:-1]).values\n    change = change[np.nonzero(change)[0]]\n    change = change[~np.isnan(change)]\n    change = change[change != -np.inf]\n    change = change[change != np.inf]\n    return np.mean(change)\n\ndef add_trend_feature(arr, abs_values=False):\n    idx = np.array(range(len(arr)))\n    if abs_values:\n        arr = np.abs(arr)\n    lr = LinearRegression()\n    lr.fit(idx.reshape(-1, 1), arr)\n    return lr.coef_[0]\n\ndef featuresx(x_data, X, segment):\n    X.loc[segment, 'ave'] = x_data.mean()\n    X.loc[segment, 'std'] = x_data.std()\n    X.loc[segment, 'max'] = x_data.max()\n    X.loc[segment, 'min'] = x_data.min()\n    \n    X.loc[segment, 'q01'] = np.quantile(x_data,0.01)\n    X.loc[segment, 'q05'] = np.quantile(x_data,0.05)\n    X.loc[segment, 'q95'] = np.quantile(x_data,0.95)\n    X.loc[segment, 'q99'] = np.quantile(x_data,0.99)\n    \n    X.loc[segment, 'mean_change_abs'] = np.mean(np.diff(x_data))\n    X.loc[segment, 'mean_change_rate'] = calc_change_rate(x_data)\n    X.loc[segment, 'abs_max'] = np.abs(x_data).max()\n    X.loc[segment, 'abs_min'] = np.abs(x_data).min()\n\n    X.loc[segment, 'max_to_min'] = x_data.max() / np.abs(x_data.min())\n    X.loc[segment, 'max_to_min_diff'] = x_data.max() - np.abs(x_data.min())\n    X.loc[segment, 'count_big'] = len(x_data[np.abs(x_data) > 500])\n    X.loc[segment, 'sum'] = x_data.sum()\n\n    X.loc[segment, 'abs_trend'] = add_trend_feature(x_data, abs_values=True)\n    X.loc[segment, 'abs_mean'] = np.abs(x_data).mean()\n    X.loc[segment, 'abs_std'] = np.abs(x_data).std()\n    X.loc[segment, 'abs_median'] = np.median(np.abs(x_data))\n\n    X.loc[segment, 'trend'] = add_trend_feature(x_data)\n    X.loc[segment, 'mad'] = x_data.mad()\n    X.loc[segment, 'kurt'] = x_data.kurtosis()\n    X.loc[segment, 'skew'] = x_data.skew()\n    X.loc[segment, 'med'] = x_data.median()\n    \n    X.loc[segment, 'abs_q95'] = np.quantile(np.abs(x_data),0.95)\n    X.loc[segment, 'abs_q99'] = np.quantile(np.abs(x_data),0.99)\n    X.loc[segment, 'F_test'], X.loc[segment, 'p_test'] = stats.f_oneway(x_data[:30000],x_data[30000:60000],x_data[60000:90000],x_data[90000:120000],x_data[120000:])\n    X.loc[segment, 'av_change_abs'] = np.mean(np.diff(x_data))\n        \n    return X\n\ndef feature_extraction(data,segments):\n    X = pd.DataFrame(index=range(segments), dtype=np.float64)\n    y = pd.DataFrame(index=range(segments), dtype=np.float64,columns=['time_to_failure'])\n    for segment in tqdm(range(segments)):\n        seg = data.iloc[segment*rows:segment*rows+rows]\n        x_data = pd.Series(seg['acoustic_data'].values)\n        y_data = seg['time_to_failure'].values[-1]\n\n        y.loc[segment,'time_to_failure'] = y_data\n        \n        X = featuresx(x_data, X, segment)\n        \n    return X, y\n    ","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:03:42.268667Z","iopub.execute_input":"2024-04-20T10:03:42.269114Z","iopub.status.idle":"2024-04-20T10:03:42.300831Z","shell.execute_reply.started":"2024-04-20T10:03:42.269075Z","shell.execute_reply":"2024-04-20T10:03:42.299090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rows = 150_000\nsegments = int(np.floor(data.shape[0] / rows))\n\nX,y = feature_extraction(data, segments)\n\ny=y.values.flatten()\nprint(X.shape)\nprint(y.shape)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:04:19.495682Z","iopub.execute_input":"2024-04-20T10:04:19.496184Z","iopub.status.idle":"2024-04-20T10:14:08.352468Z","shell.execute_reply.started":"2024-04-20T10:04:19.496118Z","shell.execute_reply":"2024-04-20T10:14:08.350847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"###### 2.1.1 Train-test split\n<font color=#6698FF>In the below, we split the train data into a test and a train set. We played around with multiple values for test_size and ended up going for a value of 0.25 after looking at the end results for all of the values. \n","metadata":{"id":"udqAFg1C3JHe"}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nX_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25,shuffle=True, random_state=102)\nx_data = data['acoustic_data'].values[:300000000] \ny_data = data['time_to_failure'].values[:300000000]","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:17:06.181033Z","iopub.execute_input":"2024-04-20T10:17:06.181561Z","iopub.status.idle":"2024-04-20T10:17:06.193976Z","shell.execute_reply.started":"2024-04-20T10:17:06.181515Z","shell.execute_reply":"2024-04-20T10:17:06.192435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\n\n<font color=#6698FF> We chose to plot the data in order to get a visual image of what we are working with, after analysing this image we found out that there is a relationship between the peaks of the acoustic data and time to failure. We came to the conclusion that the acoustic data peaks right before the time to failure peaks. We also concluded that the time to failure is 0 when the acoustic data peaks.</font>\n","metadata":{"id":"QN5jS2_p_Y9N"}},{"cell_type":"code","source":"acoustic_data = x_data[::100]\ntime_to_failure_data = y_data[::100]\n\nfig, ax1 = plt.subplots(figsize=(12, 8))\nplt.title(\"Time to failure and Acoustic data from 500 000 000 samples\")\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_data, color='r')\nax2.set_ylabel('time to failure', color='r')\nplt.legend(['time to failure'])\nplt.grid()","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:31:03.370187Z","iopub.execute_input":"2024-04-20T10:31:03.370919Z","iopub.status.idle":"2024-04-20T10:31:10.503167Z","shell.execute_reply.started":"2024-04-20T10:31:03.370854Z","shell.execute_reply":"2024-04-20T10:31:10.501578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.3 Data Preparation\n\n<font color=#6698FF> We first used preprocessing.standardScaker().fit() to check for all the features if they are in categorial or not and we used .transform() to transform the categorial ones to numerical.\n</font>","metadata":{"id":"4A5GsVzgATRu"}},{"cell_type":"code","source":"from sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X_train)\nX_train = scaler.transform(X_train)\nX_test = scaler.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:18:11.096939Z","iopub.execute_input":"2024-04-20T10:18:11.098538Z","iopub.status.idle":"2024-04-20T10:18:11.127298Z","shell.execute_reply.started":"2024-04-20T10:18:11.098447Z","shell.execute_reply":"2024-04-20T10:18:11.125587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 3. Training and Results\n<font color=#6698FF> \nWe used multiple methods as shown below; they consist of the following methods:\n   - Linear regression\n   - Neural support vector machine\n   - Decision tree regressor\n   - Random forest regressor\n\nThe reason why we chose these algorithms is because we already used some of them in earlier ML assignments and some of the others we used because we saw others using them during the LANL earthquake prediction challenge.\n    \nThe score and mse for each algorithm is printed. The algorithm with the highest score and lowest mse is Random forest regressor so in this case that would be the best algorithm with the highest accuracy.","metadata":{"id":"IUcc2ue7yrOv"}},{"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\n#linear regression\nmodel = LinearRegression()\nmodel.fit(X_train, y_train)\npred = model.predict(X_test)\nscore = model.score(X_test, y_test)\nmse = mean_squared_error(y_test, pred) #print the mse for accuracy\nprint(score)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:23:43.142953Z","iopub.execute_input":"2024-04-20T10:23:43.144263Z","iopub.status.idle":"2024-04-20T10:23:43.177295Z","shell.execute_reply.started":"2024-04-20T10:23:43.144189Z","shell.execute_reply":"2024-04-20T10:23:43.175100Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#neural support vector machine\nnsvm = NuSVR(kernel='linear')\nnsvm.fit(X_train, y_train)\npred = nsvm.predict(X_test)\nscore = nsvm.score(X_test, y_test)\nmse = mean_squared_error(y_test, pred) #print the mse for accuracy\nprint(score)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:24:12.826080Z","iopub.execute_input":"2024-04-20T10:24:12.826659Z","iopub.status.idle":"2024-04-20T10:24:13.945029Z","shell.execute_reply.started":"2024-04-20T10:24:12.826609Z","shell.execute_reply":"2024-04-20T10:24:13.943484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#support vector machine\nsvm = NuSVR(kernel='rbf')\nsvm.fit(X_train, y_train)\npred = svm.predict(X_test)\nscore = svm.score(X_test, y_test)\nmse = mean_squared_error(y_test, pred) #print the mse for accuracy\nprint(score)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:24:16.086656Z","iopub.execute_input":"2024-04-20T10:24:16.087312Z","iopub.status.idle":"2024-04-20T10:24:16.977421Z","shell.execute_reply.started":"2024-04-20T10:24:16.087251Z","shell.execute_reply":"2024-04-20T10:24:16.975428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.tree import DecisionTreeRegressor\n#decision tree regressor\ntree = DecisionTreeRegressor(max_depth=3)\ntree.fit(X_train, y_train)\npred = tree.predict(X_test)\nscore = tree.score(X_test, y_test)\nmse = mean_squared_error(y_test, pred) #print the mse for accuracy\nprint(score)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:24:18.313554Z","iopub.execute_input":"2024-04-20T10:24:18.314030Z","iopub.status.idle":"2024-04-20T10:24:18.603914Z","shell.execute_reply.started":"2024-04-20T10:24:18.313988Z","shell.execute_reply":"2024-04-20T10:24:18.602435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\n#random forest regressor\nRFC = RandomForestRegressor(max_depth=5, random_state=102)\nRFC.fit(X_train, y_train)\npred = RFC.predict(X_test)\nscore = RFC.score(X_test, y_test)\nmse = mean_squared_error(y_test, pred)\nprint(score)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2024-04-20T10:43:23.604274Z","iopub.execute_input":"2024-04-20T10:43:23.604730Z","iopub.status.idle":"2024-04-20T10:43:25.668370Z","shell.execute_reply.started":"2024-04-20T10:43:23.604690Z","shell.execute_reply":"2024-04-20T10:43:25.666982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion\n\n<font color=#6698FF>\nWe used a pretty safe and common approach. We pretty much followed the steps we learned during the ML cours. One of the strengths of this approach is the amount of features we extract from the data. We've also tried multiple regression models which gives us a better insight. However the regression models used are pretty basic so incorporating more advanced methods would probably be better. We also did not use a neural network which might have improved the accuracy immensely.\n</font>","metadata":{"id":"yekgnxu7y0gE"}},{"cell_type":"markdown","source":"## 5. References","metadata":{"id":"wn7PBEmWESIi"}},{"cell_type":"markdown","source":"* https://www.kaggle.com/competitions/LANL-Earthquake-Prediction\n* https://www.kaggle.com/code/jessteijn/ml23-week8\n* https://www.kaggle.com/code/seppedijkstra/machine-learning-week-8","metadata":{}}]}