{"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","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport collections\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\n\n# Any results you write to the current directory are saved as output.\n\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nimport scipy.stats as stats\n\nfrom tsfresh.feature_extraction import feature_calculators\n\n\n# train.csv is huge, so I implement csv_fragments() function\n# which yields DataFrame of the specified length while scaning a csv file from start to end.\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-28T21:55:33.732116Z","iopub.execute_input":"2023-04-28T21:55:33.732336Z","iopub.status.idle":"2023-04-28T21:55:34.313600Z","shell.execute_reply.started":"2023-04-28T21:55:33.732291Z","shell.execute_reply":"2023-04-28T21:55:34.312794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 1. Introduction\nNames: Fynn van der Wal, Davi Spiller Beltrao \\\nUsernames: Fynn van der Wal, Davi Spiller \\\nScore: 1.716 \\\nPublic Score: 2.64457\n","metadata":{}},{"cell_type":"markdown","source":"## 2. Data\n","metadata":{}},{"cell_type":"markdown","source":"### 2.1 Loading of the data","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float64})\npd.options.display.precision = 15\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:55:34.315025Z","iopub.execute_input":"2023-04-28T21:55:34.315477Z","iopub.status.idle":"2023-04-28T21:59:27.152401Z","shell.execute_reply.started":"2023-04-28T21:55:34.315279Z","shell.execute_reply":"2023-04-28T21:59:27.151654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.2 Data Exploration\n\nBelow, the first 10 milion datapoints are plotted to get a sense of the data. The input (X) seems to be some kind of acoustic data. The output (Y) seems to be a countdown until the earthquake event. In the first 10 milion rows, there is only one instance where the time_to_failue reaches 0, indicating an earthquake in the laboratory. This time point is preceded by some minor oscillations until a large peak in the signal is reached. After this peak, there is a period of smaller oscillations before the actual earthquake takes place.</font>","metadata":{}},{"cell_type":"code","source":"# plot first milion datapoints\nfig, ax = plt.subplots(2,1, figsize=(20,12))\nn = 10000000\nax[0].plot(train.index.values[:n], train.time_to_failure.values[:n], c=\"darkred\")\nax[0].set_title(\"time to failure\")\nax[0].set_xlabel(\"Index\")\nax[0].set_ylabel(\"time to failure in ms\");\nax[1].plot(train.index.values[:n], train.acoustic_data.values[:n], c=\"mediumseagreen\")\nax[1].set_title(\"Acoustic data\")\nax[1].set_xlabel(\"Index\")\nax[1].set_ylabel(\"Acoustic Signal\");","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:59:27.154727Z","iopub.execute_input":"2023-04-28T21:59:27.155239Z","iopub.status.idle":"2023-04-28T21:59:37.155488Z","shell.execute_reply.started":"2023-04-28T21:59:27.155187Z","shell.execute_reply":"2023-04-28T21:59:37.154816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Lets now plot the entire dataset in a nice graph:**","metadata":{}},{"cell_type":"code","source":"acoustic_data = train['acoustic_data'].values[::100] # sample at every 100th data point\ntime_data = train['time_to_failure'].values[::100] # sample at every 100th data point","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:59:37.156852Z","iopub.execute_input":"2023-04-28T21:59:37.157142Z","iopub.status.idle":"2023-04-28T21:59:37.161818Z","shell.execute_reply.started":"2023-04-28T21:59:37.157075Z","shell.execute_reply":"2023-04-28T21:59:37.160755Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax1 = plt.subplots(figsize=(16, 4))\nplt.title(\"Acoustic data and the time to earthquake\")\nplt.plot(acoustic_data, color='g')\nax1.set_ylabel('acoustic_data', color='g')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(time_data, color='r')\nax2.set_ylabel('time_to_failure', color='r')\nplt.legend(['time_to_failure'], loc=(0.875, 0.8))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:59:37.163248Z","iopub.execute_input":"2023-04-28T21:59:37.163847Z","iopub.status.idle":"2023-04-28T21:59:46.230749Z","shell.execute_reply.started":"2023-04-28T21:59:37.163799Z","shell.execute_reply":"2023-04-28T21:59:46.229925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2.3 Data Preperation","metadata":{}},{"cell_type":"markdown","source":"#### 2.3.1 Feature extraction\nFor the feature extraction [this](https://www.kaggle.com/code/artgor/even-more-features) notebook was used as a guide. We used a number of features ranging from very basic (average, std, max, min) to more complex ones like rolling features using different sizes of windows. We also included the number of peaks and different quantiles becease we saw these being used by the [price winning team](https://www.kaggle.com/competitions/LANL-Earthquake-Prediction/discussion/94390) for this challenge. Our strategy regarding the features is to try as much features as possible and then use PCA as feature selection to see what sticks. Because the dataset is large and it consists of multiple earthquakes that are uncorrelated to our understanding, the data was split up in segments of 150,000 and the features were calculated on these segments. Also, the avg, std, min and max where calculated on different parts of these segments. For example for the first 50000 and the last 50000 datapoints of the segment.","metadata":{}},{"cell_type":"markdown","source":"**Features used by winning team:**\n1. number of peaks of at least support 2 on the denoised signal, \n2. 20% percentile on std of rolling window of size 50, \n3. 4th Mel-frequency cepstral coefficients mean\n4. 18th Mel-frequency cepstral coefficients mean","metadata":{}},{"cell_type":"code","source":"import librosa, librosa.display","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:59:46.232118Z","iopub.execute_input":"2023-04-28T21:59:46.232808Z","iopub.status.idle":"2023-04-28T21:59:47.653203Z","shell.execute_reply.started":"2023-04-28T21:59:46.232531Z","shell.execute_reply":"2023-04-28T21:59:47.652001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#https://www.kaggle.com/code/artgor/even-more-features was used as a guide on which features to use.\n\nrows = 150_000\nsegments = int(np.floor(train.shape[0] / rows))\n\nX = pd.DataFrame(index=range(segments), dtype=np.float64)\nY = pd.DataFrame(index=range(segments), dtype=np.float64, columns=['time_to_failure'])\n\nfor segment in tqdm(range(segments)):\n    seg = train.iloc[segment*rows:segment*rows+rows]\n    x = pd.Series(seg['acoustic_data'].values)\n    y = seg['time_to_failure'].values[-1]\n    \n    Y.loc[segment, 'time_to_failure'] = y\n    X.loc[segment, 'ave'] = x.mean()\n    X.loc[segment, 'std'] = x.std()\n    X.loc[segment, 'max'] = x.max()\n    X.loc[segment, 'min'] = x.min()\n    \n    autocorr_lag = 5\n    X.loc[segment, 'autocorrelation_5'] = feature_calculators.autocorrelation(x, autocorr_lag)\n    X.loc[segment, 'num_peaks_10'] = feature_calculators.number_peaks(x, 10)\n    \n    X.loc[segment, 'q99'] = np.quantile(x,0.99)\n    X.loc[segment, 'q95'] = np.quantile(x, 0.95)\n    X.loc[segment, 'q05'] = np.quantile(x, 0.05)\n    X.loc[segment, 'q1'] = np.quantile(x,0.01)\n    \n    X.loc[segment, 'std_down50000'] = x[:50000].std()\n    X.loc[segment, 'std_up50000'] = x[-50000:].std()\n    X.loc[segment, 'std_down10000'] = x[:10000].std()\n    X.loc[segment, 'std_up10000'] = x[-10000:].std()\n    \n    X.loc[segment, 'avg_down50000'] = x[:50000].mean()\n    X.loc[segment, 'avg_up50000'] = x[-50000:].mean()\n    X.loc[segment, 'avg_down10000'] = x[:10000].mean()\n    X.loc[segment, 'avg_up10000'] = x[-10000:].mean()\n    \n    X.loc[segment, 'min_down50000'] = x[:50000].min()\n    X.loc[segment, 'min_up50000'] = x[-50000:].min()\n    X.loc[segment, 'min_down10000'] = x[:10000].min()\n    X.loc[segment, 'min_up10000'] = x[-10000:].min()\n    \n    X.loc[segment, 'max_down50000'] = x[:50000].max()\n    X.loc[segment, 'max_up50000'] = x[-50000:].max()\n    X.loc[segment, 'max_down10000'] = x[:10000].max()\n    X.loc[segment, 'max_up10000'] = x[-10000:].max()\n    \n    for window in [10, 20, 50, 100, 500, 1000]:\n        x_window_std = x.rolling(window).std().dropna().values\n        x_window_mean = x.rolling(window).mean().dropna().values\n        \n        X.loc[segment, 'ave_window_std' + str(window)] = x_window_std.mean()\n        X.loc[segment, 'std_window_std' + str(window)] = x_window_std.std()\n        X.loc[segment, 'max_window_std' + str(window)] = x_window_std.max()\n        X.loc[segment, 'min_window_std' + str(window)] = x_window_std.min()\n        X.loc[segment, 'q01_window_std' + str(window)] = np.quantile(x_window_std, 0.01)\n        X.loc[segment, 'q05_window_std' + str(window)] = np.quantile(x_window_std, 0.05)\n        X.loc[segment, 'q20_window_std' + str(window)] = np.quantile(x_window_std, 0.20)\n        X.loc[segment, 'q95_window_std' + str(window)] = np.quantile(x_window_std, 0.95)\n        X.loc[segment, 'q99_window_std' + str(window)] = np.quantile(x_window_std, 0.99)\n        X.loc[segment, 'diff_window_std' + str(window)] = np.mean(np.diff(x_window_std))\n        \n        X.loc[segment, 'ave_window_std' + str(window)] = x_window_mean.mean()\n        X.loc[segment, 'std_window_std' + str(window)] = x_window_mean.std()\n        X.loc[segment, 'max_window_std' + str(window)] = x_window_mean.max()\n        X.loc[segment, 'min_window_std' + str(window)] = x_window_mean.min()\n        X.loc[segment, 'q01_window_std' + str(window)] = np.quantile(x_window_mean, 0.01)\n        X.loc[segment, 'q05_window_std' + str(window)] = np.quantile(x_window_mean, 0.05)\n        X.loc[segment, 'q20_window_std' + str(window)] = np.quantile(x_window_mean, 0.20)\n        X.loc[segment, 'q95_window_std' + str(window)] = np.quantile(x_window_mean, 0.95)\n        X.loc[segment, 'q99_window_std' + str(window)] = np.quantile(x_window_mean, 0.99)\n        X.loc[segment, 'diff_window_std' + str(window)] = np.mean(np.diff(x_window_mean))\n\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T21:59:47.654928Z","iopub.execute_input":"2023-04-28T21:59:47.655240Z","iopub.status.idle":"2023-04-28T22:11:47.007465Z","shell.execute_reply.started":"2023-04-28T21:59:47.655191Z","shell.execute_reply":"2023-04-28T22:11:47.006629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 2.3.2 Normalization\nThe data is normalized.","metadata":{}},{"cell_type":"code","source":"# Normalization\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(X)\nX = scaler.transform(X)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.008891Z","iopub.execute_input":"2023-04-28T22:11:47.009345Z","iopub.status.idle":"2023-04-28T22:11:47.023513Z","shell.execute_reply.started":"2023-04-28T22:11:47.009164Z","shell.execute_reply":"2023-04-28T22:11:47.022542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 2.3.3  PCA and Train/Test split\n The test size is set to 0.3. PCA is performed on the feature matrix. A bar plot is generated and a step plot to display the explained variance ratio of each principal component (PC).","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nn = 10\npca = PCA(n_components=n) # doing pca and keeping only n_components\npca = pca.fit(X) # the correct dimension of X for sklearn is P*N (samples*features)\nX_pca = pca.transform(X)\nprint(X_pca.shape)\n\n# Splitting the data into train and test data\nfrom sklearn.model_selection import train_test_split\nx_train, x_test, y_train, y_test = train_test_split(X_pca, Y, test_size=0.3,shuffle=False, random_state=102)\n\n# Standard normalize the data\nfrom sklearn import preprocessing\nscaler = preprocessing.StandardScaler().fit(x_train)\nx_train = scaler.transform(x_train)\nx_test = scaler.transform(x_test)\n\nplt.figure(figsize=(7,7))\n# perform pca on features\nplt.bar(range(0,n), pca.explained_variance_ratio_, label=\"individual var\");\nplt.step(range(0,n), np.cumsum(pca.explained_variance_ratio_),'r', label=\"cumulative var\");\nplt.xlabel('Principal component index'); plt.ylabel('explained variance ratio');\nplt.legend()\n\n\n\n#convert y values to categorical values\n#y_train = np.ravel(y_train)\n#y_test = np.ravel(y_test)\n#label = preprocessing.LabelEncoder()\n#y_train = label.fit_transform(y_train)\n#y_test = label.fit_transform(y_test)\n","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.024927Z","iopub.execute_input":"2023-04-28T22:11:47.025384Z","iopub.status.idle":"2023-04-28T22:11:47.917153Z","shell.execute_reply.started":"2023-04-28T22:11:47.025192Z","shell.execute_reply":"2023-04-28T22:11:47.834781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n## 3. Training and Results\nFor the linear model, Support Vector Regression is used. For the non-linear model, a Random Forest is used. Both models are tested on the train data, as well as the test data.","metadata":{}},{"cell_type":"code","source":"#Linear training model\n\nfrom sklearn import datasets, svm\nfrom sklearn.svm import SVR\n\nmodel1 = SVR()\n\n#Non-linear training model\n\nfrom sklearn.ensemble import RandomForestRegressor\n\nmodel2 = RandomForestRegressor(max_depth=20)\n\n# fit the model\nmodel1.fit(x_train, y_train.values.flatten())\nmodel2.fit(x_train, y_train.values.flatten())\n\ny_pred_te1 = model1.predict(x_test)\ny_pred_tr1 = model1.predict(x_train)\ny_pred_te2 = model2.predict(x_test)\ny_pred_tr2 = model2.predict(x_train)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.835904Z","iopub.status.idle":"2023-04-28T22:11:47.839124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.1 Linear Model","metadata":{}},{"cell_type":"code","source":"#Scores of the linear model\nfrom sklearn import metrics\nscore1 = model1.score(x_test, y_test)\nprint(\"Test score of linear model:\", score1)\n\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time vs. actual time to failure\")\nplt.plot(y_pred_te1, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1200)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.844388Z","iopub.status.idle":"2023-04-28T22:11:47.847581Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"score2 = model1.score(x_train, y_train)\nprint(\"Train score of linear model:\", score2)\n\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time vs. actual time to failure\")\nplt.plot(y_pred_tr1, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1200)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.850918Z","iopub.status.idle":"2023-04-28T22:11:47.851705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### 3.2 Non-Linear Model","metadata":{}},{"cell_type":"code","source":"#Scores of the non-linear model\nscore3 = model2.score(x_test, y_test) \nprint(\"Test score of non-linear model:\", score3)\n\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time vs. actual time to failure\")\nplt.plot(y_pred_te2, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1200)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.857130Z","iopub.status.idle":"2023-04-28T22:11:47.857888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"score4 = model2.score(x_train, y_train)\nprint(\"Train score of non-linear model:\", score4)\n\nfig, ax1 = plt.subplots(figsize=(20, 4))\n\nplt.title(\"Comparison of the predicted time vs. actual time to failure\")\nplt.plot(y_pred_tr2, color='r')\nax1.set_ylabel('Time to failure', color='r')\nplt.legend(['Predicted time to failure'], loc=(0.08, 0.9))\nax2 = ax1.twinx()\nplt.plot(Y, color='b')\nax2.set_ylabel('Time to failure', color='b')\nplt.legend(['Actual time to failure'], loc=(0.08, 0.85))\nplt.grid(True)\nplt.xlim(0, 1200)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.858726Z","iopub.status.idle":"2023-04-28T22:11:47.859411Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It can clearly be seen that the results from the train data are much better than the test data, this suggests an overtrained model on both the linear and non-linear models. It can also be observed that the linear model performs better overall. The linear model is also much more computationally efficient.","metadata":{}},{"cell_type":"markdown","source":"## 4. Discussion and Conclusion","metadata":{}},{"cell_type":"markdown","source":"To get the best earthquake predictions, the following steps were made:\n\n1. Data exploration\n* The first 10 milion samples were plotted to get a sense of thet data.\n2. Data preprocessing\n* The features were trained on segments of 150,000 samples .\n* The feature matrix was normalized\n* The data was split into train and test with a testsize of 0.3.\n* PCA was done on the features to reduce the size\n3. Training models\n* SVM was used as a linear model, Random Forest was used as a non-linear model. SVM performed slightly better than Random Forest.\n\n**Strong points of our approach are:**\n* testing many features \\\nImprovements could be:\n* using k-fold validation\n* using different models, maybe neural networks\n* diving deeper into earthquake prediction to learn new important features\n* make more plots with different features to deepen our understanding of them\n","metadata":{}},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id') # submission csv\nX_test = pd.DataFrame(dtype=np.float64, index=submission.index)\n    \nfor seg_id in X_test.index:\n    seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n    \n    x = pd.Series(seg['acoustic_data'].values)\n    \n    X_test.loc[seg_id, 'ave'] = x.mean()\n    X_test.loc[seg_id, 'std'] = x.std()\n    X_test.loc[seg_id, 'max'] = x.max()\n    X_test.loc[seg_id, 'min'] = x.min()\n    \n    #mfcc = librosa.feature.mfcc(x.astype('float32'), sr=40000, n_mfcc = 20) \n    #X.loc[segment, 'mfcc4'] = mfcc[4, :].mean()\n    #X.loc[segment, 'mfcc18'] = mfcc[18, :].mean()\n    \n    autocorr_lag = 5\n    X_test.loc[seg_id, 'autocorrelation_5'] = feature_calculators.autocorrelation(x, autocorr_lag)\n    X_test.loc[seg_id, 'num_peaks_10'] = feature_calculators.number_peaks(x, 10)\n    \n    X_test.loc[seg_id, 'q99'] = np.quantile(x,0.99)\n    X_test.loc[seg_id, 'q95'] = np.quantile(x, 0.95)\n    X_test.loc[seg_id, 'q05'] = np.quantile(x, 0.05)\n    X_test.loc[seg_id, 'q1'] = np.quantile(x,0.01)\n    \n    X_test.loc[seg_id, 'std_down50000'] = x[:50000].std()\n    X_test.loc[seg_id, 'std_up50000'] = x[-50000:].std()\n    X_test.loc[seg_id, 'std_down10000'] = x[:10000].std()\n    X_test.loc[seg_id, 'std_up10000'] = x[-10000:].std()\n    \n    X_test.loc[seg_id, 'avg_down50000'] = x[:50000].mean()\n    X_test.loc[seg_id, 'avg_up50000'] = x[-50000:].mean()\n    X_test.loc[seg_id, 'avg_down10000'] = x[:10000].mean()\n    X_test.loc[seg_id, 'avg_up10000'] = x[-10000:].mean()\n    \n    X_test.loc[seg_id, 'min_down50000'] = x[:50000].min()\n    X_test.loc[seg_id, 'min_up50000'] = x[-50000:].min()\n    X_test.loc[seg_id, 'min_down10000'] = x[:10000].min()\n    X_test.loc[seg_id, 'min_up10000'] = x[-10000:].min()\n    \n    X_test.loc[seg_id, 'max_down50000'] = x[:50000].max()\n    X_test.loc[seg_id, 'max_up50000'] = x[-50000:].max()\n    X_test.loc[seg_id, 'max_down10000'] = x[:10000].max()\n    X_test.loc[seg_id, 'max_up10000'] = x[-10000:].max()\n    \n    for window in [10, 20, 50, 100, 500, 1000]:\n        x_window_std = x.rolling(window).std().dropna().values\n        x_window_mean = x.rolling(window).mean().dropna().values\n        \n        X_test.loc[seg_id, 'ave_window_std' + str(window)] = x_window_std.mean()\n        X_test.loc[seg_id, 'std_window_std' + str(window)] = x_window_std.std()\n        X_test.loc[seg_id, 'max_window_std' + str(window)] = x_window_std.max()\n        X_test.loc[seg_id, 'min_window_std' + str(window)] = x_window_std.min()\n        X_test.loc[seg_id, 'q01_window_std' + str(window)] = np.quantile(x_window_std, 0.01)\n        X_test.loc[seg_id, 'q05_window_std' + str(window)] = np.quantile(x_window_std, 0.05)\n        X_test.loc[seg_id, 'q20_window_std' + str(window)] = np.quantile(x_window_std, 0.20)\n        X_test.loc[seg_id, 'q95_window_std' + str(window)] = np.quantile(x_window_std, 0.95)\n        X_test.loc[seg_id, 'q99_window_std' + str(window)] = np.quantile(x_window_std, 0.99)\n        X_test.loc[seg_id, 'diff_window_std' + str(window)] = np.mean(np.diff(x_window_std))\n        \n        X_test.loc[seg_id, 'ave_window_std' + str(window)] = x_window_mean.mean()\n        X_test.loc[seg_id, 'std_window_std' + str(window)] = x_window_mean.std()\n        X_test.loc[seg_id, 'max_window_std' + str(window)] = x_window_mean.max()\n        X_test.loc[seg_id, 'min_window_std' + str(window)] = x_window_mean.min()\n        X_test.loc[seg_id, 'q01_window_std' + str(window)] = np.quantile(x_window_mean, 0.01)\n        X_test.loc[seg_id, 'q05_window_std' + str(window)] = np.quantile(x_window_mean, 0.05)\n        X_test.loc[seg_id, 'q20_window_std' + str(window)] = np.quantile(x_window_mean, 0.20)\n        X_test.loc[seg_id, 'q95_window_std' + str(window)] = np.quantile(x_window_mean, 0.95)\n        X_test.loc[seg_id, 'q99_window_std' + str(window)] = np.quantile(x_window_mean, 0.99)\n        X_test.loc[seg_id, 'diff_window_std' + str(window)] = np.mean(np.diff(x_window_mean))","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.860317Z","iopub.status.idle":"2023-04-28T22:11:47.861000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler = preprocessing.StandardScaler().fit(X_test)\nX_test = scaler.transform(X_test)\n\npca = PCA(n_components=10)\npca = pca.fit(X_test)\nX_pca = pca.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.861837Z","iopub.status.idle":"2023-04-28T22:11:47.862543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# submit to the submission file\nscaler = preprocessing.StandardScaler().fit(X_pca)\nX_pca = scaler.transform(X_pca)\nsubmission['time_to_failure'] = model1.predict(X_pca) \nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2023-04-28T22:11:47.863394Z","iopub.status.idle":"2023-04-28T22:11:47.864078Z"},"trusted":true},"execution_count":null,"outputs":[]}]}