{"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":"# **LANL Earthquake Prediction**\n### Data used: LANL Earthquake Prediction ","metadata":{}},{"cell_type":"markdown","source":"# Nội dung\n\n- Mô tả bài toán\n- Chuẩn bị dữ liệu \n- Trích lọc đặc trưng \n- Mô tả các đặc trưng\n- Huấn luyện mô hình\n- Bagging Regression \n- XGBoost\n- Submission","metadata":{}},{"cell_type":"markdown","source":"#### Back to top\n#### Next\n# Mô tả bài toán\n

Mục tiêu của bài toán là sử dụng sóng địa chấn nhằm dự đoán thời gian xảy ra động đất trong phòng thí nghiệm. Dữ liệu thu được từ những thí nghiệm nổi tiếng nhằm nghiên cứu những tính chất vật lý của động đất. Các nhà nghiên cứu mô phỏng lại các trận động đất trong phòng thí nghiệm với mục đích thu lại dữ liệu về các sóng địa chấn tại mỗi thời điểm và thời gian còn lại đến khi xảy ra động đất tại mỗi thời điểm đó. acoustic_data là dữ liệu về sóng địa chấn sử dụng để dự đoán thời gian còn lại trước khi trận động đất tiếp theo xảy ra (time_to_failure).

\n\n

Dữ liệu training là những phần nhỏ (segment), liên tục của dữ liệu thí nghiệm. Dữ liệu test gồm một thư mục chứa nhiều file mỗi file cũng là những segment nhỏ chứa dữ liệu liên tục nhưng không phải là một phần của dữ liệu thí nghiệm. Vì vậy, việc dự đoán không thể giả sử là tuân theo những quy luật thông thường như có thể thấy trong file dữ liệu train.

\n\n

Với mỗi seg_id trong thư mục test, ta cần dự đoán được time_to_failure duy nhất tương ứng với thời gian từ thời điểm cùng trong segment cho đến khi động đất xảy ra.

","metadata":{}},{"cell_type":"markdown","source":"### Mô tả các file dữ liệu\n*

train.csv - Một đoạn liên tục duy nhất lấy từ dữ liệu thí nghiệm.

\n*

test - Một thư mục chứa nhiều segment của dữ liệu test với cùng độ dài 150000 mỗi segment.

\n*

sample_submission.csv - Mẫu file submission.

","metadata":{}},{"cell_type":"markdown","source":"### Các trường dữ liệu\n*

acoustic_data - Cường độ sóng địa chấn [int16]

\n*

time_to_failure - Thời gian (giây) cho đến khi xảy ra động đất [float64]

\n*

seg_id - id của các segment trong tập dữ liệu test

","metadata":{}},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n#### Next\n# Chuẩn bị dữ liệu\n### Tải các packages cần thiết\n

Nạp các package cần thiết để sử dụng, trích lọc đặc trưng từ dữ liệu và train model.

","metadata":{}},{"cell_type":"code","source":"import gc\nimport os\nimport time\nimport logging\nimport librosa\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport xgboost as xgb\nfrom scipy import stats\nfrom tqdm import tqdm_notebook\nimport matplotlib.pyplot as plt\nfrom sklearn.ensemble import BaggingRegressor\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.model_selection import GridSearchCV, KFold","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:11:33.541762Z","iopub.execute_input":"2021-06-11T04:11:33.542245Z","iopub.status.idle":"2021-06-11T04:11:36.15778Z","shell.execute_reply.started":"2021-06-11T04:11:33.542137Z","shell.execute_reply":"2021-06-11T04:11:36.1566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Không hiển thị các warning messages","metadata":{}},{"cell_type":"code","source":"warnings.filterwarnings(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:11:36.159288Z","iopub.execute_input":"2021-06-11T04:11:36.159623Z","iopub.status.idle":"2021-06-11T04:11:36.163934Z","shell.execute_reply.started":"2021-06-11T04:11:36.159591Z","shell.execute_reply":"2021-06-11T04:11:36.16303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Nạp dữ liệu","metadata":{}},{"cell_type":"code","source":"test_path = '../input/LANL-Earthquake-Prediction/test'\nsubmission_path = '../input/LANL-Earthquake-Prediction/sample_submission.csv'\ntrain_path = '../input/LANL-Earthquake-Prediction/train.csv'","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:11:38.240903Z","iopub.execute_input":"2021-06-11T04:11:38.241255Z","iopub.status.idle":"2021-06-11T04:11:38.245412Z","shell.execute_reply.started":"2021-06-11T04:11:38.241227Z","shell.execute_reply":"2021-06-11T04:11:38.24431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Đọc dữ liệu từ file train.csv

","metadata":{}},{"cell_type":"code","source":"%%time\nraw_df = pd.read_csv(train_path,\n dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:11:39.851888Z","iopub.execute_input":"2021-06-11T04:11:39.852199Z","iopub.status.idle":"2021-06-11T04:15:04.050139Z","shell.execute_reply.started":"2021-06-11T04:11:39.852173Z","shell.execute_reply":"2021-06-11T04:15:04.048475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"raw_df.shape","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:16:57.668344Z","iopub.execute_input":"2021-06-11T04:16:57.668986Z","iopub.status.idle":"2021-06-11T04:16:57.67978Z","shell.execute_reply.started":"2021-06-11T04:16:57.668941Z","shell.execute_reply":"2021-06-11T04:16:57.678737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Lấy mẫu dữ liệu test

","metadata":{}},{"cell_type":"code","source":"sample_test = pd.read_csv('../input/LANL-Earthquake-Prediction/test/seg_017314.csv')\nsample_test.shape","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:17:01.593618Z","iopub.execute_input":"2021-06-11T04:17:01.594224Z","iopub.status.idle":"2021-06-11T04:17:01.642095Z","shell.execute_reply.started":"2021-06-11T04:17:01.594178Z","shell.execute_reply":"2021-06-11T04:17:01.641037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n#### Next\n# Trích lọc đặc trưng\n

Vì các đoạn trong tập dữ liệu test đều có độ dài 150000 nên dữ liệu train cũng nên được chia thành các phần nhỏ với cùng độ dài 150000

","metadata":{}},{"cell_type":"code","source":"seg_size = 150000\nsegments = int(np.floor(raw_df.shape[0] / seg_size))\nsegments","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:17:04.169923Z","iopub.execute_input":"2021-06-11T04:17:04.170584Z","iopub.status.idle":"2021-06-11T04:17:04.178972Z","shell.execute_reply.started":"2021-06-11T04:17:04.170532Z","shell.execute_reply":"2021-06-11T04:17:04.177695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Plot segment đầu tiên và time_to_failure tương ứng

","metadata":{}},{"cell_type":"code","source":"sample = raw_df['acoustic_data'][:150000]\nttf = raw_df['time_to_failure'][:150000]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title('first segment and time_to_failure')\nplt.plot(sample, color='b')\nax1.set_xlabel('training samples')\nax1.set_ylabel('acoustic data')\nplt.legend(['acoustic data'], loc=(0.01, 0.95))\nax2 = ax1.twinx()\nplt.plot(ttf, color='y')\nax2.set_ylabel('time to failure')\nplt.legend(['time to failure'], loc=(0.01, 0.9))\nplt.grid(True)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:17:06.733687Z","iopub.execute_input":"2021-06-11T04:17:06.734084Z","iopub.status.idle":"2021-06-11T04:17:07.296254Z","shell.execute_reply.started":"2021-06-11T04:17:06.734052Z","shell.execute_reply":"2021-06-11T04:17:07.295463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.DataFrame(index=range(segments), dtype=np.float64)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:17:10.641751Z","iopub.execute_input":"2021-06-11T04:17:10.642375Z","iopub.status.idle":"2021-06-11T04:17:10.647788Z","shell.execute_reply.started":"2021-06-11T04:17:10.642317Z","shell.execute_reply":"2021-06-11T04:17:10.646762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Vì dữ liệu có dạng sóng, nên ta sẽ trích lọc các đặc trưng từ phổ biên độ và trong miền số phức bằng cách sử dụng Fast Fourier Transform (np.fft). Với mỗi segment dài 150000, hàm dưới đây sẽ trích lọc ra các đặc trưng thống kê như sau:

\n

mean - Giá trị trung bình

\n

std - Độ lệch chuẩn

\n

max - Giá trị max

\n

min - Giá trị min

\n

sum - Tổng các phần tử trong segment

\n

mad - Mean absolute deviation - Trung bình độ lệch tuyệt đối

\n

kurt - Kurtosis (mô tả hình dạng đuôi của phân phối xác suất)

\n

skew - Skewness (mô tả tính bất đối xứng của phân phối xác suất)

\n

med - Phần tử trung vị

\n

abs_mean - Trung bình giá trị tuyệt đối

\n

q95, q99, q05, q01 - Phân vị 0.95, 0.99, 0.05 và 0.01

\n\n

Ngoài ra, ứng với mỗi segment, chúng ta sẽ trích lọc đặc trưng MFCC sử dụng thư viện librosa. Sau đó, ứng với mỗi frame trong MFCC, ta sẽ lấy giá trị trung bình để có được vector MFCC dùng cho việc dự đoán. Việc trích lọc đặc trưng MFCC sử dụng hàm biến đổi Fourier cho thời gian ngắn Short Time Fourier Transform hay STFT với các thông số n_fft = 2048, hop_length = 512","metadata":{}},{"cell_type":"code","source":"def create_features(segment_id, segment, X, targets=True):\n xc = pd.Series(segment['acoustic_data'].values)\n zc = np.fft.fft(xc)\n realFFT = np.real(zc)\n imagFFT = np.imag(zc)\n \n X.loc[seg_id, 'mean'] = xc.mean()\n X.loc[seg_id, 'std'] = xc.std()\n X.loc[seg_id, 'max'] = xc.max()\n X.loc[seg_id, 'min'] = xc.min()\n X.loc[seg_id, 'sum'] = xc.sum()\n X.loc[seg_id, 'mad'] = xc.mad()\n X.loc[seg_id, 'kurt'] = xc.kurtosis()\n X.loc[seg_id, 'skew'] = xc.skew()\n X.loc[seg_id, 'med'] = xc.median()\n \n X.loc[seg_id, 'q95'] = np.quantile(xc, 0.95)\n X.loc[seg_id, 'q99'] = np.quantile(xc, 0.99)\n X.loc[seg_id, 'q05'] = np.quantile(xc, 0.05)\n X.loc[seg_id, 'q01'] = np.quantile(xc, 0.01)\n \n X.loc[seg_id, 'abs_mean'] = np.abs(xc).mean()\n X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n X.loc[seg_id, 'abs_min'] = np.abs(xc).min()\n X.loc[seg_id, 'abs_q95'] = np.quantile(np.abs(xc), 0.95)\n X.loc[seg_id, 'abs_q99'] = np.quantile(np.abs(xc), 0.99)\n X.loc[seg_id, 'abs_q05'] = np.quantile(np.abs(xc), 0.05)\n X.loc[seg_id, 'abs_q01'] = np.quantile(np.abs(xc), 0.01)\n \n X.loc[seg_id, 'Rmean'] = realFFT.mean()\n X.loc[seg_id, 'Rstd'] = realFFT.std()\n X.loc[seg_id, 'Rmax'] = realFFT.max()\n X.loc[seg_id, 'Rmin'] = realFFT.min()\n \n X.loc[seg_id, 'Imean'] = imagFFT.mean()\n X.loc[seg_id, 'Istd'] = imagFFT.std()\n X.loc[seg_id, 'Imax'] = imagFFT.max()\n X.loc[seg_id, 'Imin'] = imagFFT.min()\n \n X.loc[seg_id, 'std_first_50000'] = xc[:50000].std()\n X.loc[seg_id, 'std_last_50000'] = xc[-50000:].std()\n X.loc[seg_id, 'std_first_25000'] = xc[:25000].std()\n X.loc[seg_id, 'std_last_25000'] = xc[-25000:].std()\n X.loc[seg_id, 'std_first_10000'] = xc[:10000].std()\n X.loc[seg_id, 'std_last_10000'] = xc[-10000:].std()\n \n X.loc[seg_id, 'Rstd_first_50000'] = realFFT[:50000].std()\n X.loc[seg_id, 'Rstd_last_50000'] = realFFT[-50000:].std()\n X.loc[seg_id, 'Rstd_first_25000'] = realFFT[:25000].std()\n X.loc[seg_id, 'Rstd_last_25000'] = realFFT[-25000:].std()\n X.loc[seg_id, 'Rstd_first_10000'] = realFFT[:10000].std()\n X.loc[seg_id, 'Rstd_last_10000'] = realFFT[-10000:].std()\n \n X.loc[seg_id, 'mean_first_50000'] = xc[:50000].mean()\n X.loc[seg_id, 'mean_last_50000'] = xc[-50000:].mean()\n X.loc[seg_id, 'mean_first_25000'] = xc[:25000].mean()\n X.loc[seg_id, 'mean_last_25000'] = xc[-25000:].mean()\n X.loc[seg_id, 'mean_first_10000'] = xc[:10000].mean()\n X.loc[seg_id, 'mean_last_10000'] = xc[-10000:].mean()\n \n X.loc[seg_id, 'Rmean_first_50000'] = realFFT[:50000].mean()\n X.loc[seg_id, 'Rmean_last_50000'] = realFFT[-50000:].mean()\n X.loc[seg_id, 'Rmean_first_25000'] = realFFT[:25000].mean()\n X.loc[seg_id, 'Rmean_last_25000'] = realFFT[-25000:].mean()\n X.loc[seg_id, 'Rmean_first_10000'] = realFFT[:10000].mean()\n X.loc[seg_id, 'Rmean_last_10000'] = realFFT[-10000:].mean()\n \n X.loc[seg_id, 'max_first_50000'] = xc[:50000].max()\n X.loc[seg_id, 'max_last_50000'] = xc[-50000:].max()\n X.loc[seg_id, 'max_first_25000'] = xc[:25000].max()\n X.loc[seg_id, 'max_last_25000'] = xc[-25000:].max()\n X.loc[seg_id, 'max_first_10000'] = xc[:10000].max()\n X.loc[seg_id, 'max_last_10000'] = xc[-10000:].max()\n \n X.loc[seg_id, 'Rmax_first_50000'] = realFFT[:50000].max()\n X.loc[seg_id, 'Rmax_last_50000'] = realFFT[-50000:].max()\n X.loc[seg_id, 'Rmax_first_25000'] = realFFT[:25000].max()\n X.loc[seg_id, 'Rmax_last_25000'] = realFFT[-25000:].max()\n X.loc[seg_id, 'Rmax_first_10000'] = realFFT[:10000].max()\n X.loc[seg_id, 'Rmax_last_10000'] = realFFT[-10000:].max()\n \n X.loc[seg_id, 'min_first_50000'] = xc[:50000].min()\n X.loc[seg_id, 'min_last_50000'] = xc[-50000:].min()\n X.loc[seg_id, 'min_first_25000'] = xc[:25000].min()\n X.loc[seg_id, 'min_last_25000'] = xc[-25000:].min()\n X.loc[seg_id, 'min_first_10000'] = xc[:10000].min()\n X.loc[seg_id, 'min_last_10000'] = xc[-10000:].min()\n \n X.loc[seg_id, 'Rmin_first_50000'] = realFFT[:50000].min()\n X.loc[seg_id, 'Rmin_last_50000'] = realFFT[-50000:].min()\n X.loc[seg_id, 'Rmin_first_25000'] = realFFT[:25000].min()\n X.loc[seg_id, 'Rmin_last_25000'] = realFFT[-25000:].min()\n X.loc[seg_id, 'Rmin_first_10000'] = realFFT[:10000].min()\n X.loc[seg_id, 'Rmin_last_10000'] = realFFT[-10000:].min()\n \n for windows in [5, 10, 50, 100, 500, 1000, 5000, 10000]:\n x_roll_std = xc.rolling(windows).std().dropna().values\n x_roll_mean = xc.rolling(windows).mean().dropna().values\n X.loc[seg_id, 'ave_roll_std_' + str(windows)] = x_roll_std.mean()\n X.loc[seg_id, 'std_roll_std_' + str(windows)] = x_roll_std.std()\n X.loc[seg_id, 'max_roll_std_' + str(windows)] = x_roll_std.max()\n X.loc[seg_id, 'min_roll_std_' + str(windows)] = x_roll_std.min()\n X.loc[seg_id, 'q01_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.01)\n X.loc[seg_id, 'q05_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.05)\n X.loc[seg_id, 'q95_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.95)\n X.loc[seg_id, 'q99_roll_std_' + str(windows)] = np.quantile(x_roll_std, 0.99)\n X.loc[seg_id, 'av_change_abs_roll_std_' + str(windows)] = np.mean(np.diff(x_roll_std))\n X.loc[seg_id, 'av_change_rate_roll_std_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n X.loc[seg_id, 'abs_max_roll_std_' + str(windows)] = np.abs(x_roll_std).max()\n \n X.loc[seg_id, 'ave_roll_mean_' + str(windows)] = x_roll_mean.mean()\n X.loc[seg_id, 'std_roll_mean_' + str(windows)] = x_roll_mean.std()\n X.loc[seg_id, 'max_roll_mean_' + str(windows)] = x_roll_mean.max()\n X.loc[seg_id, 'min_roll_mean_' + str(windows)] = x_roll_mean.min()\n X.loc[seg_id, 'q01_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.01)\n X.loc[seg_id, 'q05_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.05)\n X.loc[seg_id, 'q95_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.95)\n X.loc[seg_id, 'q99_roll_mean_' + str(windows)] = np.quantile(x_roll_mean, 0.99)\n \n X.loc[seg_id, 'max_to_min_diff'] = xc.max() - np.abs(xc.min())\n \n X.loc[seg_id, 'count_big'] = len(xc[np.abs(xc) > xc.mean()])\n \n mfcc = librosa.feature.mfcc(np.asarray(xc.astype(np.float16)))\n mfcc_mean = mfcc.mean(axis=1)\n for i in range(len(mfcc_mean)):\n X.loc[seg_id, 'mfcc_' + str(i)] = mfcc_mean[i]\n \n if targets:\n X.loc[seg_id, 'time_to_failure'] = segment['time_to_failure'].values[-1]","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:17:14.336507Z","iopub.execute_input":"2021-06-11T04:17:14.336908Z","iopub.status.idle":"2021-06-11T04:17:14.386186Z","shell.execute_reply.started":"2021-06-11T04:17:14.336876Z","shell.execute_reply":"2021-06-11T04:17:14.385215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Tính các feature và thêm vào dataframe train_df

","metadata":{}},{"cell_type":"code","source":"for seg_id in tqdm_notebook(range(segments)):\n segment = raw_df.iloc[seg_id*seg_size:seg_id*seg_size+seg_size]\n create_features(seg_id, segment, train_df)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:17:16.165704Z","iopub.execute_input":"2021-06-11T04:17:16.166349Z","iopub.status.idle":"2021-06-11T04:39:50.63645Z","shell.execute_reply.started":"2021-06-11T04:17:16.166297Z","shell.execute_reply":"2021-06-11T04:39:50.635464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Lấy các giá trị (seg_id) từ file submission mẫu và khởi tạo dataframe test_df

","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv(submission_path, index_col='seg_id')\ntest_df = pd.DataFrame(columns=train_df.columns, dtype=np.float64, index=submission.index)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:48:45.71427Z","iopub.execute_input":"2021-06-11T04:48:45.714772Z","iopub.status.idle":"2021-06-11T04:48:45.763901Z","shell.execute_reply.started":"2021-06-11T04:48:45.714731Z","shell.execute_reply":"2021-06-11T04:48:45.762656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Do lượng feature lớn nên ta sẽ sử dụng correlation matrix để lựa chọn ra những feature có độ tương quan đối với time_to_failure cao

","metadata":{}},{"cell_type":"code","source":"corrMat = train_df.corr()\ncorrTarget = abs(corrMat['time_to_failure'])\ncorrTarget = corrTarget[corrTarget > 0.5]\nhigh_corr_features = corrTarget.index.drop(['time_to_failure'])\nhigh_corr_features","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:48:47.964823Z","iopub.execute_input":"2021-06-11T04:48:47.96526Z","iopub.status.idle":"2021-06-11T04:48:48.689993Z","shell.execute_reply.started":"2021-06-11T04:48:47.965227Z","shell.execute_reply":"2021-06-11T04:48:48.688907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Khởi tạo giá trị cho X_train và y_train","metadata":{}},{"cell_type":"code","source":"X_train = train_df[high_corr_features]\ny_train = train_df['time_to_failure']\nX_train","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:48:51.982887Z","iopub.execute_input":"2021-06-11T04:48:51.983354Z","iopub.status.idle":"2021-06-11T04:48:52.041604Z","shell.execute_reply.started":"2021-06-11T04:48:51.983314Z","shell.execute_reply":"2021-06-11T04:48:52.040583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Tính toán các đặc trưng từ dữ liệu test

","metadata":{}},{"cell_type":"code","source":"for seg_id in tqdm_notebook(test_df.index):\n seg_file = seg_id + '.csv'\n seg = pd.read_csv(os.path.join(test_path, seg_file))\n create_features(seg_id, seg, test_df, targets=False)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T04:48:56.649927Z","iopub.execute_input":"2021-06-11T04:48:56.650342Z","iopub.status.idle":"2021-06-11T05:02:26.781088Z","shell.execute_reply.started":"2021-06-11T04:48:56.650309Z","shell.execute_reply":"2021-06-11T05:02:26.780048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Loại bỏ các feature có độ tương quan thấp khỏi test_df","metadata":{}},{"cell_type":"code","source":"final_X_test = test_df[high_corr_features]\nfinal_X_test","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:08:54.614912Z","iopub.execute_input":"2021-06-11T05:08:54.615611Z","iopub.status.idle":"2021-06-11T05:08:54.65842Z","shell.execute_reply.started":"2021-06-11T05:08:54.615549Z","shell.execute_reply":"2021-06-11T05:08:54.657419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n#### Next\n# Mô tả các đặc trưng","metadata":{}},{"cell_type":"markdown","source":"Hàm bên dưới scatter các điểm dữ liệu theo từng feature và time_to_failure tương ứng","metadata":{}},{"cell_type":"code","source":"def plot_feature_scatter(features, X=train_df):\n i = 0\n plt.figure()\n nlines = int(len(features)/2) \n fig, ax = plt.subplots(nlines, 2, figsize=(20, 5*nlines ))\n for feature in features:\n i+=1\n plt.subplot(nlines,2,i)\n ax[int((i-1)/2)][(i+1) % 2].set_xlabel(feature)\n ax[int((i-1)/2)][(i+1) % 2].set_ylabel('time_to_failure')\n plt.title('{} - time_to_falure correlation)'.format(feature), color='r')\n plt.scatter(x = X[feature], y = X['time_to_failure'])\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:08:57.890443Z","iopub.execute_input":"2021-06-11T05:08:57.89089Z","iopub.status.idle":"2021-06-11T05:08:57.89969Z","shell.execute_reply.started":"2021-06-11T05:08:57.890848Z","shell.execute_reply":"2021-06-11T05:08:57.89853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_feature_scatter(high_corr_features)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:08:58.706247Z","iopub.execute_input":"2021-06-11T05:08:58.706663Z","iopub.status.idle":"2021-06-11T05:09:04.251157Z","shell.execute_reply.started":"2021-06-11T05:08:58.706627Z","shell.execute_reply":"2021-06-11T05:09:04.249907Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Đoạn code dưới sử dụng hàm sns.histplot mô tả histogram và density (mật độ) của các đặc trưng, màu đỏ ứng với dữ liệu train và màu xanh ứng với dữ liệu test

","metadata":{}},{"cell_type":"code","source":"def plot_histplot_features(features, colors=['red', 'blue'], df1=X_train, df2=final_X_test):\n i = 0\n plt.figure()\n nlines = int(len(features)/2)\n fig, ax = plt.subplots(nlines,2,figsize=(20,5*nlines))\n for feature in features:\n i += 1\n plt.subplot(nlines,2,i)\n sns.histplot(df1[feature],color=colors[0], kde=True, bins=40, label='train')\n sns.histplot(df2[feature],color=colors[1], kde=True, bins=40, label='test')\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:07.81868Z","iopub.execute_input":"2021-06-11T05:09:07.819086Z","iopub.status.idle":"2021-06-11T05:09:07.827654Z","shell.execute_reply.started":"2021-06-11T05:09:07.819051Z","shell.execute_reply":"2021-06-11T05:09:07.826201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_histplot_features(high_corr_features)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:08.926702Z","iopub.execute_input":"2021-06-11T05:09:08.927128Z","iopub.status.idle":"2021-06-11T05:09:20.111802Z","shell.execute_reply.started":"2021-06-11T05:09:08.927082Z","shell.execute_reply":"2021-06-11T05:09:20.111167Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Features scaling\n

Đoạn code dưới sẽ scale lại các feature để phù hợp hơn với phân phối chuẩn và sẽ visualize lại bằng các đồ thị mới. Chúng ta sẽ scale với cả dữ liệu train và test. Xanh lá ứng với dữ liệu train và xanh dương ứng với dữ liệu test.

","metadata":{}},{"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(pd.concat([X_train, final_X_test]))\nscaled_X_train = pd.DataFrame(scaler.transform(X_train), columns=X_train.columns)\nscaled_X_test = pd.DataFrame(scaler.transform(final_X_test), columns=final_X_test.columns, index=final_X_test.index)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:23.713876Z","iopub.execute_input":"2021-06-11T05:09:23.714316Z","iopub.status.idle":"2021-06-11T05:09:23.73584Z","shell.execute_reply.started":"2021-06-11T05:09:23.714264Z","shell.execute_reply":"2021-06-11T05:09:23.734716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_histplot_features(high_corr_features, colors=['green', 'blue'], df1=scaled_X_train, df2=scaled_X_test)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:28.078053Z","iopub.execute_input":"2021-06-11T05:09:28.078497Z","iopub.status.idle":"2021-06-11T05:09:39.169639Z","shell.execute_reply.started":"2021-06-11T05:09:28.078453Z","shell.execute_reply":"2021-06-11T05:09:39.168689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Mô tả các features và time_to_failure\n

Đoạn code dưới mô tả các features và time_to_failure tương ứng trên cùng một đồ thị.

","metadata":{}},{"cell_type":"code","source":"def plot_feature_ttf(features):\n i = 0\n plt.figure()\n nlines = int(len(features)/2)\n fig, ax = plt.subplots(nlines,2,figsize=(32,8*nlines))\n for feature in features:\n i += 1\n plt.subplot(nlines,2,i)\n plt.title('({}) and time to failure'.format(feature))\n plt.plot(X_train[feature], color='r')\n ax[int((i-1)/2)][(i+1) % 2].set_xlabel('training samples')\n ax[int((i-1)/2)][(i+1) % 2].set_ylabel('acoustic data ({})'.format(feature), color='r')\n plt.legend(['acoustic data ({})'.format(feature)], loc=(0.01, 0.95))\n ax2 = ax[int((i-1)/2)][(i+1) % 2].twinx()\n plt.plot(y_train, color='b')\n ax2.set_ylabel('time to failure', color='b')\n plt.legend(['time to failure'], loc=(0.01, 0.9))\n plt.grid(True)\n plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:43.126565Z","iopub.execute_input":"2021-06-11T05:09:43.126954Z","iopub.status.idle":"2021-06-11T05:09:43.137044Z","shell.execute_reply.started":"2021-06-11T05:09:43.126919Z","shell.execute_reply":"2021-06-11T05:09:43.136214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_feature_ttf(high_corr_features)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:43.934158Z","iopub.execute_input":"2021-06-11T05:09:43.934797Z","iopub.status.idle":"2021-06-11T05:09:53.616712Z","shell.execute_reply.started":"2021-06-11T05:09:43.934758Z","shell.execute_reply":"2021-06-11T05:09:53.615951Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n#### Next\n# Huấn luyện mô hình","metadata":{}},{"cell_type":"markdown","source":"

Dữ liệu sau khi đã được scale sẽ được đưa vào huấn luyện mô hình SVR để dự đoán. Mô hình được đánh giá dựa vào RMSE và MAE.

","metadata":{}},{"cell_type":"code","source":"from sklearn.svm import SVR\nfrom sklearn.metrics import mean_squared_error, mean_absolute_error","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:09:58.654148Z","iopub.execute_input":"2021-06-11T05:09:58.654733Z","iopub.status.idle":"2021-06-11T05:09:58.658117Z","shell.execute_reply.started":"2021-06-11T05:09:58.654697Z","shell.execute_reply":"2021-06-11T05:09:58.657499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Chia tập dữ liệu train thành dữ liệu train và test

","metadata":{}},{"cell_type":"code","source":"X_train, X_test, y_train, y_test = train_test_split(scaled_X_train, y_train, test_size=0.33, random_state=42)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:00.018446Z","iopub.execute_input":"2021-06-11T05:10:00.019001Z","iopub.status.idle":"2021-06-11T05:10:00.029561Z","shell.execute_reply.started":"2021-06-11T05:10:00.018968Z","shell.execute_reply":"2021-06-11T05:10:00.028643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Sử dụng *GridSearchCV* để lựa chọn tham số cho mô hình","metadata":{}},{"cell_type":"code","source":"# parameters = {'kernel': ('linear', 'rbf','poly'), 'C':[0.01, 0.1, 1, 10, 100],'epsilon':[0.1,0.2,0.3,0.4,0.5,0.6]}\n# svr = SVR(gamma = 'scale')\n# clf = GridSearchCV(svr, parameters)\n# clf.fit(X_train,y_train)\n# clf.best_params_","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:02.778831Z","iopub.execute_input":"2021-06-11T05:10:02.779399Z","iopub.status.idle":"2021-06-11T05:10:02.782656Z","shell.execute_reply.started":"2021-06-11T05:10:02.779365Z","shell.execute_reply":"2021-06-11T05:10:02.781704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Sau khi dùng GridSearchCV thu được bộ tham số: kernel='linear', C=100, epsilon=0.6","metadata":{}},{"cell_type":"code","source":"model = SVR(kernel='linear', C=100, epsilon=0.6, gamma='scale')\nmodel.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:03.89452Z","iopub.execute_input":"2021-06-11T05:10:03.895126Z","iopub.status.idle":"2021-06-11T05:10:23.419627Z","shell.execute_reply.started":"2021-06-11T05:10:03.895091Z","shell.execute_reply":"2021-06-11T05:10:23.418891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Dự đoán trên tập X_test","metadata":{}},{"cell_type":"code","source":"y_pred = model.predict(X_test)\nrmse = mean_squared_error(y_test, y_pred, squared=False)\nmae = mean_absolute_error(y_pred, y_test)\nprint('rmse = {} \\n mae = {}'.format(rmse,mae))","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:23.420754Z","iopub.execute_input":"2021-06-11T05:10:23.421146Z","iopub.status.idle":"2021-06-11T05:10:23.509953Z","shell.execute_reply.started":"2021-06-11T05:10:23.421118Z","shell.execute_reply":"2021-06-11T05:10:23.509077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n#### Next\n# Bagging Regression","metadata":{}},{"cell_type":"markdown","source":"Mô tả sai số của các dự đoán","metadata":{}},{"cell_type":"code","source":"chart = pd.DataFrame()\nchart['target'] = y_test\nchart['pred'] = y_pred\nchart['off-target'] = abs(chart['target'] - chart['pred'])","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:28.903051Z","iopub.execute_input":"2021-06-11T05:10:28.903479Z","iopub.status.idle":"2021-06-11T05:10:28.912695Z","shell.execute_reply.started":"2021-06-11T05:10:28.903439Z","shell.execute_reply":"2021-06-11T05:10:28.911407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"chart.describe()","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:30.134268Z","iopub.execute_input":"2021-06-11T05:10:30.134716Z","iopub.status.idle":"2021-06-11T05:10:30.160739Z","shell.execute_reply.started":"2021-06-11T05:10:30.134676Z","shell.execute_reply":"2021-06-11T05:10:30.159638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Có thể thấy bias của kết quả dự đoán tương đối chuẩn tuy nhiên sai số vẫn khá cao. Ta sẽ sử dụng kỹ thuật Bagging cho bài toán regression để làm giảm sai số của mô hình.\nBase model được chọn là mô hình SVR(kernel='rbf',gamma='scale',C=1.0,epsilon=0.1)","metadata":{}},{"cell_type":"code","source":"model = BaggingRegressor(base_estimator=SVR(), n_estimators=10, random_state=42).fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:32.902335Z","iopub.execute_input":"2021-06-11T05:10:32.902747Z","iopub.status.idle":"2021-06-11T05:10:35.070514Z","shell.execute_reply.started":"2021-06-11T05:10:32.902713Z","shell.execute_reply":"2021-06-11T05:10:35.069474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Có thể thấy kết quả được cải thiện đôi chút, tuy nhiên không đáng kể","metadata":{}},{"cell_type":"code","source":"y_pred = model.predict(X_test)\nrmse = mean_squared_error(y_test, y_pred, squared=False)\nmae = mean_absolute_error(y_pred, y_test)\nprint('rmse = {} \\n mae = {}'.format(rmse,mae))","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:40.554453Z","iopub.execute_input":"2021-06-11T05:10:40.554882Z","iopub.status.idle":"2021-06-11T05:10:42.476841Z","shell.execute_reply.started":"2021-06-11T05:10:40.554848Z","shell.execute_reply":"2021-06-11T05:10:42.476048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n#### Next\n# XGBoost","metadata":{}},{"cell_type":"markdown","source":"

SVR với Bagging cho kết quả tương đối tốt trên tập train. Tuy nhiên khi dùng để dự đoán với dữ liệu test của cuộc thi thì MAE tăng lên đáng kể. Nguyên nhân có thể là do đã xảy ra hiện tượng overfitting.

\n

Chúng ta sẽ sử dụng mô hình huấn luyện XGBoost để có thể giảm thiểu hiện tượng overfitting đối với tập dữ liệu train.

","metadata":{}},{"cell_type":"code","source":"X_train = scaled_X_train\ny_train = train_df['time_to_failure']","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:10:47.27822Z","iopub.execute_input":"2021-06-11T05:10:47.278621Z","iopub.status.idle":"2021-06-11T05:10:47.28363Z","shell.execute_reply.started":"2021-06-11T05:10:47.278586Z","shell.execute_reply":"2021-06-11T05:10:47.282853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Sử dụng K-fold cross validation với số fold là 5 cho việc huấn luyện và đánh giá mô hình

","metadata":{}},{"cell_type":"code","source":"%%time\nn_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=42)\n \npredictions = np.zeros(len(scaled_X_test))\nxgb_params = {'eta': 0.05,\n 'max_depth': 10,\n 'subsample': 0.8,\n 'objective': 'reg:linear',\n 'eval_metric': 'mae',\n 'nthread': 4}\n#run model\nfor fold_, (trn_idx, val_idx) in enumerate(folds.split(X_train,y_train.values)):\n X_tr, X_val = X_train.iloc[trn_idx], X_train.iloc[val_idx]\n y_tr, y_val = y_train.iloc[trn_idx], y_train.iloc[val_idx]\n\n train_data = xgb.DMatrix(data=X_tr, label=y_tr, feature_names=X_train.columns)\n valid_data = xgb.DMatrix(data=X_val, label=y_val, feature_names=X_train.columns)\n\n watchlist = [(valid_data, 'valid_data')]\n xgb_model = xgb.train(params=xgb_params, dtrain=train_data, num_boost_round=1000, evals=watchlist, early_stopping_rounds=40, verbose_eval=False)\n\n y_pred = xgb_model.predict(xgb.DMatrix(data=X_val, feature_names=X_tr.columns))\n mae = mean_absolute_error(y_pred, y_val.values)\n print('fold {}: MAE = {}'.format(fold_, mae))\n \n predictions += xgb_model.predict(xgb.DMatrix(scaled_X_test, feature_names=scaled_X_test.columns))/folds.n_splits","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:51:58.862488Z","iopub.execute_input":"2021-06-11T05:51:58.862839Z","iopub.status.idle":"2021-06-11T05:52:08.46131Z","shell.execute_reply.started":"2021-06-11T05:51:58.862799Z","shell.execute_reply":"2021-06-11T05:52:08.460574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Back to top\n#### Previous\n# Submission","metadata":{}},{"cell_type":"markdown","source":"XGBoost Submission","metadata":{}},{"cell_type":"code","source":"submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\nsubmission['time_to_failure'] = predictions\nsubmission.to_csv('submission_xgb.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-06-11T05:54:05.010448Z","iopub.execute_input":"2021-06-11T05:54:05.010832Z","iopub.status.idle":"2021-06-11T05:54:05.038783Z","shell.execute_reply.started":"2021-06-11T05:54:05.010801Z","shell.execute_reply":"2021-06-11T05:54:05.038107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"SVM Bagging Submission","metadata":{}},{"cell_type":"code","source":"# submission = pd.read_csv('../input/LANL-Earthquake-Prediction/sample_submission.csv')\n# predictions = model.predict(scaled_X_test)\n# submission['time_to_failure'] = predictions\n# submission.to_csv('submission_svm.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}