{"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":"# 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=#bcaeff>Jesse Wessteijn, Jessteijn <br>\n                    Evert-Jan Beiboer, evertjanbeiboer <br>\n                    score: , Leaderboard rank:","metadata":{}},{"cell_type":"markdown","source":"## 2. Data\n","metadata":{"id":"PB7FLdQn-dmK"}},{"cell_type":"markdown","source":"### 2.1 Dataset","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":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-30T23:25:31.583798Z","iopub.execute_input":"2023-04-30T23:25:31.584497Z","iopub.status.idle":"2023-04-30T23:25:33.211934Z","shell.execute_reply.started":"2023-04-30T23:25:31.584440Z","shell.execute_reply":"2023-04-30T23:25:33.210563Z"},"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":"2023-04-30T23:25:35.353348Z","iopub.execute_input":"2023-04-30T23:25:35.353810Z","iopub.status.idle":"2023-04-30T23:29:27.871057Z","shell.execute_reply.started":"2023-04-30T23:25:35.353766Z","shell.execute_reply":"2023-04-30T23:29:27.869737Z"},"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":"2023-04-30T23:29:27.872986Z","iopub.execute_input":"2023-04-30T23:29:27.873357Z","iopub.status.idle":"2023-04-30T23:29:27.895790Z","shell.execute_reply.started":"2023-04-30T23:29:27.873297Z","shell.execute_reply":"2023-04-30T23:29:27.894783Z"},"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":"2023-04-30T23:29:27.897132Z","iopub.execute_input":"2023-04-30T23:29:27.897576Z","iopub.status.idle":"2023-04-30T23:36:20.811799Z","shell.execute_reply.started":"2023-04-30T23:29:27.897541Z","shell.execute_reply":"2023-04-30T23:36:20.810050Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#bcaeff>We used to [this notebook](https://www.kaggle.com/code/artgor/earthquakes-fe-more-features-and-samples) for feature extration. I created a function for it, so I can use it later on the test data files.","metadata":{}},{"cell_type":"markdown","source":"#### 2.1.1 Train-test split","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.2,shuffle=True, random_state=102)\nprint(X_train.shape)\nprint(X_test.shape)\nprint(y_train.shape)\nprint(y_test.shape)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:20.816944Z","iopub.execute_input":"2023-04-30T23:36:20.817801Z","iopub.status.idle":"2023-04-30T23:36:20.834222Z","shell.execute_reply.started":"2023-04-30T23:36:20.817739Z","shell.execute_reply":"2023-04-30T23:36:20.833019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#bcaeff>First the test value can not be too big or too small. When the test value is too big you have nothing left for training. When the test value is too small the given accuracy will not be reliable.<br>\n<font color=#bcaeff> For the second question 'random_state' is used so that output will not change every restart.","metadata":{}},{"cell_type":"markdown","source":"### 2.2 Data Exploration","metadata":{"id":"QN5jS2_p_Y9N"}},{"cell_type":"code","source":"n=1000\ny1_acoustic_data=data['acoustic_data'][::n]\ny2_time_to_failure=data['time_to_failure'][::n]\nprint(data[0:10*n:n])\n\n\n#plotting the graphs of all the deta\nfig, ax = plt.subplots(2, 1, figsize=(27,9))\nax[0].plot(y2_time_to_failure, color='b')\nax[0].legend(['Real time to fail'])\nax[0].set_ylabel('Time left')\nax[0].set_xlabel('Sample times 1000)')\nax[1].plot(y1_acoustic_data, color='r')\nax[1].legend(['Acoustic data'])\nax[1].set_ylabel('Acoustic output' )\nax[1].set_xlabel('Sample times 1000)')","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:20.835133Z","iopub.execute_input":"2023-04-30T23:36:20.835434Z","iopub.status.idle":"2023-04-30T23:36:22.840394Z","shell.execute_reply.started":"2023-04-30T23:36:20.835405Z","shell.execute_reply":"2023-04-30T23:36:22.839407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#bcaeff>if you plot the data you can see that the acoustic value has a peak when time of failure is zero. This is the measured earthquake","metadata":{}},{"cell_type":"markdown","source":"### 2.3 Data Preparation","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":"2023-04-30T23:36:22.841806Z","iopub.execute_input":"2023-04-30T23:36:22.842568Z","iopub.status.idle":"2023-04-30T23:36:22.857010Z","shell.execute_reply.started":"2023-04-30T23:36:22.842530Z","shell.execute_reply":"2023-04-30T23:36:22.855798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 3. Training and Results","metadata":{"id":"IUcc2ue7yrOv"}},{"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\n# define the model\nmodel = LinearRegression()\n# fit the model\nmodel.fit(X_train, y_train)\n# predict\npredictions=model.predict(X_test)\n#accuracy\nscore=model.score(X_test, y_test)\nprint(score)\nmse = mean_squared_error(y_test, predictions)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:22.858847Z","iopub.execute_input":"2023-04-30T23:36:22.859767Z","iopub.status.idle":"2023-04-30T23:36:22.887074Z","shell.execute_reply.started":"2023-04-30T23:36:22.859715Z","shell.execute_reply":"2023-04-30T23:36:22.885434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#define and fit model\nnsvm = NuSVR(kernel='linear')\nnsvm.fit(X_train, y_train)\n\n# predict\npredictions=nsvm.predict(X_test)\n# accuracy\nscore=nsvm.score(X_test, y_test)\nprint(score)\nmse = mean_squared_error(y_test, predictions)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:22.888951Z","iopub.execute_input":"2023-04-30T23:36:22.890404Z","iopub.status.idle":"2023-04-30T23:36:24.152362Z","shell.execute_reply.started":"2023-04-30T23:36:22.890351Z","shell.execute_reply":"2023-04-30T23:36:24.151175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#define and fit model\nsvm = NuSVR(kernel='rbf')\nsvm.fit(X_train, y_train)\n\n# predict\npredictions=svm.predict(X_test)\n# accuracy\nscore=svm.score(X_test, y_test)\nprint(score)\nmse = mean_squared_error(y_test, predictions)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:24.153776Z","iopub.execute_input":"2023-04-30T23:36:24.154248Z","iopub.status.idle":"2023-04-30T23:36:25.044434Z","shell.execute_reply.started":"2023-04-30T23:36:24.154213Z","shell.execute_reply":"2023-04-30T23:36:25.043189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.tree import DecisionTreeRegressor\n\ntree = DecisionTreeRegressor(max_depth=3)\ntree.fit(X_train, y_train)\npredictions = tree.predict(X_test)\nscore=tree.score(X_test, y_test)\nprint(score)\nmse = mean_squared_error(y_test, predictions)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:25.048536Z","iopub.execute_input":"2023-04-30T23:36:25.048902Z","iopub.status.idle":"2023-04-30T23:36:25.227587Z","shell.execute_reply.started":"2023-04-30T23:36:25.048866Z","shell.execute_reply":"2023-04-30T23:36:25.226384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#bcaeff>I check a lot of models, those 4 are the best in my opinion. They give the best results. Two linear models and two non-linear models were used. The non-linear svm model is the best one.","metadata":{}},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion","metadata":{"id":"yekgnxu7y0gE"}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv', index_col='seg_id')\nX_test = pd.DataFrame(dtype=np.float64, index=submission.index)\n\nfor seg_id in tqdm(X_test.index):\n    seg = pd.read_csv('../input/LANL-Earthquake-Prediction/test/' + seg_id + '.csv')\n    \n    x = pd.Series(seg['acoustic_data'].values)\n    X_test = featuresx(x, X_test, seg_id)\n    \nprint(x.shape)\nprint(X_test.shape)\n\nX_test = scaler.transform(X_test)\nsubmission['time_to_failure'] = svm.predict(X_test)\nsubmission.to_csv('submission.csv')\n\nprint(submission)","metadata":{"execution":{"iopub.status.busy":"2023-04-30T23:36:25.230272Z","iopub.execute_input":"2023-04-30T23:36:25.231548Z","iopub.status.idle":"2023-04-30T23:41:22.944024Z","shell.execute_reply.started":"2023-04-30T23:36:25.231498Z","shell.execute_reply":"2023-04-30T23:41:22.943021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color=#bcaeff>We also learned to work with big data, which leads into not being able to plot everything or make a useful table to see what you doing. This had made the understanding of the deta much harder. Furthermore our results are all under 50% which is not good, but we did not do the project to achieve a certain performance but to learn with the concept of machine learning.</font>","metadata":{}}]}