{"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\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","metadata":{"id":"GF4WsCxg-lL9"}},{"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-11-17T23:45:42.755517Z","iopub.execute_input":"2023-11-17T23:45:42.756985Z","iopub.status.idle":"2023-11-17T23:45:44.289504Z","shell.execute_reply.started":"2023-11-17T23:45:42.756932Z","shell.execute_reply":"2023-11-17T23:45:44.288451Z"},"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-11-18T00:12:50.526999Z","iopub.execute_input":"2023-11-18T00:12:50.527659Z","iopub.status.idle":"2023-11-18T00:16:47.685520Z","shell.execute_reply.started":"2023-11-18T00:12:50.527606Z","shell.execute_reply":"2023-11-18T00:16:47.683852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.kernel_approximation import RBFSampler\n\ndef 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 add_fourier_features(x_data, X, segment, rbfsampler):\n    # Aplicar la Transformada de Fourier\n    fft_features = np.abs(np.fft.fft(x_data))\n\n    # Seleccionar solo la primera mitad de las características FFT (debido a la simetría)\n    fft_features = fft_features[:len(fft_features) // 2]\n\n    # Utilizar RBFSampler en las características FFT\n    fft_transformed = rbfsampler.fit_transform(fft_features.reshape(-1, 1))\n\n    # Agregar las características transformadas al DataFrame\n    for i, fft_feat in enumerate(fft_transformed.T):\n        X.loc[segment, f'fft_{i}'] = fft_feat.mean()\n\n    return X\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    # Init rbf sampler\n    rbfsampler = RBFSampler(gamma=1.0, n_components=100, random_state=1)\n    X = add_fourier_features(x_data, X, segment, rbfsampler)\n    return X\n\n\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-11-18T00:19:21.938940Z","iopub.execute_input":"2023-11-18T00:19:21.939411Z","iopub.status.idle":"2023-11-18T00:19:21.976634Z","shell.execute_reply.started":"2023-11-18T00:19:21.939370Z","shell.execute_reply":"2023-11-18T00:19:21.975123Z"},"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-11-18T00:19:22.297080Z","iopub.execute_input":"2023-11-18T00:19:22.298246Z","iopub.status.idle":"2023-11-18T00:51:09.232479Z","shell.execute_reply.started":"2023-11-18T00:19:22.298169Z","shell.execute_reply":"2023-11-18T00:51:09.230439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-11-18T00:51:09.235347Z","iopub.execute_input":"2023-11-18T00:51:09.235988Z","iopub.status.idle":"2023-11-18T00:51:09.258179Z","shell.execute_reply.started":"2023-11-18T00:51:09.235942Z","shell.execute_reply":"2023-11-18T00:51:09.255927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-11-18T00:51:21.333667Z","iopub.execute_input":"2023-11-18T00:51:21.334233Z","iopub.status.idle":"2023-11-18T00:51:23.422445Z","shell.execute_reply.started":"2023-11-18T00:51:21.334189Z","shell.execute_reply":"2023-11-18T00:51:23.420736Z"},"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-11-18T00:51:24.253898Z","iopub.execute_input":"2023-11-18T00:51:24.255988Z","iopub.status.idle":"2023-11-18T00:51:24.274853Z","shell.execute_reply.started":"2023-11-18T00:51:24.255902Z","shell.execute_reply":"2023-11-18T00:51:24.272691Z"},"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-11-18T00:51:25.681783Z","iopub.execute_input":"2023-11-18T00:51:25.682415Z","iopub.status.idle":"2023-11-18T00:51:25.734032Z","shell.execute_reply.started":"2023-11-18T00:51:25.682352Z","shell.execute_reply":"2023-11-18T00:51:25.732335Z"},"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-11-18T00:51:27.336751Z","iopub.execute_input":"2023-11-18T00:51:27.337330Z","iopub.status.idle":"2023-11-18T00:51:34.511510Z","shell.execute_reply.started":"2023-11-18T00:51:27.337282Z","shell.execute_reply":"2023-11-18T00:51:34.509576Z"},"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-11-18T00:51:34.514374Z","iopub.execute_input":"2023-11-18T00:51:34.514873Z","iopub.status.idle":"2023-11-18T00:51:35.993459Z","shell.execute_reply.started":"2023-11-18T00:51:34.514831Z","shell.execute_reply":"2023-11-18T00:51:35.991935Z"},"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-11-18T00:51:35.996389Z","iopub.execute_input":"2023-11-18T00:51:35.997064Z","iopub.status.idle":"2023-11-18T00:51:36.139657Z","shell.execute_reply.started":"2023-11-18T00:51:35.996987Z","shell.execute_reply":"2023-11-18T00:51:36.137956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from xgboost.sklearn import XGBRegressor\nfrom sklearn.model_selection import GridSearchCV\n\nxgb1 = XGBRegressor()\nparameters = {'nthread':[4], #when use hyperthread, xgboost may become slower\n              'objective':['reg:squarederror'],\n              'learning_rate': [.03, 0.05, .07], #so called `eta` value\n              'max_depth': [3, 5, 6, 7],\n              'min_child_weight': [4],\n              'silent': [1],\n              'subsample': [0.7],\n              'colsample_bytree': [0.7],\n              'n_estimators': [500]}\n\nxgb_grid = GridSearchCV(xgb1,\n                        parameters,\n                        cv = 2,\n                        n_jobs = 5,\n                        verbose=True)","metadata":{"execution":{"iopub.status.busy":"2023-11-18T01:45:02.467390Z","iopub.execute_input":"2023-11-18T01:45:02.467982Z","iopub.status.idle":"2023-11-18T01:45:02.481174Z","shell.execute_reply.started":"2023-11-18T01:45:02.467923Z","shell.execute_reply":"2023-11-18T01:45:02.479180Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb_grid.fit(X_train, y_train)\npredictions = xgb_grid.predict(X_test)\nscore=xgb_grid.score(X_test, y_test)\nprint(score)\nmse = mean_squared_error(y_test, predictions)\nprint(mse)","metadata":{"execution":{"iopub.status.busy":"2023-11-18T03:21:28.653906Z","iopub.execute_input":"2023-11-18T03:21:28.655276Z","iopub.status.idle":"2023-11-18T04:09:26.839128Z","shell.execute_reply.started":"2023-11-18T03:21:28.655175Z","shell.execute_reply":"2023-11-18T04:09:26.837146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-11-18T00:04:04.350637Z","iopub.execute_input":"2023-11-18T00:04:04.351181Z","iopub.status.idle":"2023-11-18T00:10:12.430608Z","shell.execute_reply.started":"2023-11-18T00:04:04.351136Z","shell.execute_reply":"2023-11-18T00:10:12.429451Z"},"trusted":true},"execution_count":null,"outputs":[]}]}