{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 導入函式庫"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"import gc\nimport os\nimport time\nimport sys\nimport logging\nimport datetime\nimport warnings\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport xgboost as xgb\nimport lightgbm as lgb\nfrom scipy import stats\nfrom scipy.signal import hann\nfrom tqdm import tqdm_notebook\nimport matplotlib.pyplot as plt\nfrom scipy.signal import hilbert\nfrom scipy.signal import convolve\nfrom sklearn.svm import NuSVR, SVR\nfrom catboost import CatBoostRegressor\nfrom sklearn.kernel_ridge import KernelRidge\nfrom sklearn.linear_model import ElasticNet, Lasso, BayesianRidge, LassoLarsIC\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor\nfrom sklearn.model_selection import KFold,StratifiedKFold, RepeatedKFold\nfrom sklearn.model_selection import cross_val_score, train_test_split, cross_val_predict\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.preprocessing import RobustScaler\nfrom sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone\n\nwarnings.filterwarnings(\"ignore\")","execution_count":1,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 讀取資料\n我們先看一下input目錄中的文件"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"PATH=\"../input/\"\nos.listdir(PATH)","execution_count":2,"outputs":[{"output_type":"execute_result","execution_count":2,"data":{"text/plain":"['train.csv', 'sample_submission.csv', 'test']"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"input目錄中有train.csv、sample_submission.csv，以及另外一個test資料夾，再來看一下test資料夾內:"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"test有{}個檔案\".format(len(os.listdir(os.path.join(PATH, 'test' )))))","execution_count":3,"outputs":[{"output_type":"stream","text":"test有2624個檔案\n","name":"stdout"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\ntrain_df = pd.read_csv(os.path.join(PATH,'train.csv'), \n dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})","execution_count":4,"outputs":[{"output_type":"stream","text":"CPU times: user 2min 22s, sys: 10.6 s, total: 2min 33s\nWall time: 2min 34s\n","name":"stdout"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.options.display.precision = 15\ntrain_df.head()","execution_count":5,"outputs":[{"output_type":"execute_result","execution_count":5,"data":{"text/plain":" acoustic_data time_to_failure\n0 12 1.469099998474121\n1 6 1.469099998474121\n2 8 1.469099998474121\n3 5 1.469099998474121\n4 8 1.469099998474121","text/html":"
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
acoustic_datatime_to_failure
0121.469099998474121
161.469099998474121
281.469099998474121
351.469099998474121
481.469099998474121
\n
"},"metadata":{}}]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"資料型態轉換用\n```python=\nd = {'a': [1, 2], 'b': [3, 4]}\ndf = pd.DataFrame(data=d, columns=['a', 'b'], dtype=np.int16)\nprint(sys.getsizeof(df))\ndf = df.astype({\"a\":np.int32})\nprint(sys.getsizeof(df))\nprint(df.a.dtype)\nprint(df.b.dtype)\ndf.head()\n```"},{"metadata":{},"cell_type":"markdown","source":"看一下資料的最大與最小值"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(train_df['acoustic_data'].max())\nprint(train_df['acoustic_data'].min())\nprint(train_df['time_to_failure'].max())\nprint(train_df['time_to_failure'].min())","execution_count":6,"outputs":[{"output_type":"stream","text":"5444\n-5515\n16.1074\n9.5503965e-05\n","name":"stdout"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"# 先備份原始資料\ntrain_df_save = train_df.copy","execution_count":7,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"將所有資料轉為正數，只保留震幅大小，並且取log"},{"metadata":{"trusted":true},"cell_type":"code","source":"# np.log1p(x) : 計算 log(1 + x)\n# 對acoustic_data取絕對值並做平滑處理\ntrain_df.acoustic_data = np.log1p(abs(train_df.acoustic_data))","execution_count":8,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"訓練資料的大小: 629145480筆，每筆兩個屬性"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df.shape","execution_count":9,"outputs":[{"output_type":"execute_result","execution_count":9,"data":{"text/plain":"(629145480, 2)"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"觀察一下前幾筆長相:"},{"metadata":{},"cell_type":"markdown","source":"* accoustic_data: 實驗中測量的聲學信號\n* time_to_failure: 離發生failure還有幾秒"},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.options.display.precision = 15\ntrain_df.head(5)","execution_count":10,"outputs":[{"output_type":"execute_result","execution_count":10,"data":{"text/plain":" acoustic_data time_to_failure\n0 2.564949274063110 1.469099998474121\n1 1.945910215377808 1.469099998474121\n2 2.197224617004395 1.469099998474121\n3 1.791759490966797 1.469099998474121\n4 2.197224617004395 1.469099998474121","text/html":"
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
acoustic_datatime_to_failure
02.5649492740631101.469099998474121
11.9459102153778081.469099998474121
22.1972246170043951.469099998474121
31.7917594909667971.469099998474121
42.1972246170043951.469099998474121
\n
"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"# EDA (Exploratory Data Analysis)\n由於資料量太大，總共有6億多筆，故我們先每50筆採樣出1筆，plot出資料來看看"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''# 每50筆採樣1筆\ntrain_ad_sample_df = train_df['acoustic_data'].values[::50]\ntrain_ttf_sample_df = train_df['time_to_failure'].values[::50]\n\ndef plot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1% sampled data\"):\n fig, ax1 = plt.subplots(figsize=(12, 8))\n plt.title(title)\n plt.plot(train_ad_sample_df, color='r')\n ax1.set_ylabel('acoustic data', color='r')\n plt.legend(['acoustic data'], loc=(0.01, 0.95))\n ax2 = ax1.twinx()\n plt.plot(train_ttf_sample_df, 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\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df)\ndel train_ad_sample_df\ndel train_ttf_sample_df'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"gc.collect()","execution_count":11,"outputs":[{"output_type":"execute_result","execution_count":11,"data":{"text/plain":"22"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"上圖show出了2%的資料，可發現，在failure發生處，accoustic data通常有明顯的尖峰，且即使在兩次failure之間，也會有一些高峰存在。\n\n接著我們plot出頭兩次的地震 (time_to_failure為0的地方)"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''train_ad_sample_df = train_df['acoustic_data'].values[:50580000]\ntrain_ttf_sample_df = (train_df['time_to_failure'].values[:50580000])\nplot_acc_ttf_data(train_ad_sample_df, train_ttf_sample_df, title=\"Acoustic data and time to failure: 1st 2 quakes\")\ndel train_ad_sample_df\ndel train_ttf_sample_df'''","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"在這個放大的圖中，可以發現大震盪不只會出現在failure附近，在兩次failure中間也有出現較大的震盪。"},{"metadata":{},"cell_type":"markdown","source":"# 特徵工程"},{"metadata":{"trusted":true},"cell_type":"code","source":"# 查看test資料\ntest = pd.read_csv('../input/test/seg_00030f.csv')\ntest.shape","execution_count":12,"outputs":[{"output_type":"execute_result","execution_count":12,"data":{"text/plain":"(150000, 1)"},"metadata":{}}]},{"metadata":{"trusted":true},"cell_type":"code","source":"test.head()","execution_count":13,"outputs":[{"output_type":"execute_result","execution_count":13,"data":{"text/plain":" acoustic_data\n0 4\n1 0\n2 -2\n3 0\n4 2","text/html":"
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n
acoustic_data
04
10
2-2
30
42
\n
"},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"上方結果顯示一份測試資料是15萬維的，事實上，所有測試資料都是一樣維度的。\n\n我們將訓練資料切割成與測試資料相同維度。"},{"metadata":{"trusted":true},"cell_type":"code","source":"rows = 150000\nsegments = int(np.floor(train_df.shape[0] * 8 / rows))\nprint(\"Number of segments: \", segments)","execution_count":14,"outputs":[{"output_type":"stream","text":"Number of segments: 33554\n","name":"stdout"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"def 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 classic_sta_lta(x, length_sta, length_lta):\n sta = np.cumsum(x ** 2)\n # Convert to float\n sta = np.require(sta, dtype=np.float)\n # Copy for LTA\n lta = sta.copy()\n # Compute the STA and the LTA\n sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n sta /= length_sta\n lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n lta /= length_lta\n # Pad zeros\n sta[:length_lta - 1] = 0\n # Avoid division by zero by setting zero values to tiny float\n dtiny = np.finfo(0.0).tiny\n idx = lta < dtiny\n lta[idx] = dtiny\n return sta / lta","execution_count":15,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 處理訓練資料"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X = pd.DataFrame(index=range(segments), dtype=np.float32)\ntrain_y = pd.DataFrame(index=range(segments), dtype=np.float32, columns=['time_to_failure'])","execution_count":16,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# These may be needed later\n'''need_aggregated_features = True\nif need_aggregated_features:\n total_mean = train_df['acoustic_data'].mean()\n total_std = train_df['acoustic_data'].std()\n total_max = train_df['acoustic_data'].max()\n total_min = train_df['acoustic_data'].min()\n total_sum = train_df['acoustic_data'].sum()\n total_abs_sum = np.abs(train_df['acoustic_data']).sum()'''","execution_count":17,"outputs":[{"output_type":"execute_result","execution_count":17,"data":{"text/plain":"\"need_aggregated_features = True\\nif need_aggregated_features:\\n total_mean = train_df['acoustic_data'].mean()\\n total_std = train_df['acoustic_data'].std()\\n total_max = train_df['acoustic_data'].max()\\n total_min = train_df['acoustic_data'].min()\\n total_sum = train_df['acoustic_data'].sum()\\n total_abs_sum = np.abs(train_df['acoustic_data']).sum()\""},"metadata":{}}]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X.shape, train_y.shape","execution_count":29,"outputs":[{"output_type":"execute_result","execution_count":29,"data":{"text/plain":"((33554, 3), (33554, 1))"},"metadata":{}}]},{"metadata":{"trusted":true},"cell_type":"code","source":"def create_features(seg_id, seg, X):\n xc = pd.Series(seg['acoustic_data'].values)\n zc = np.fft.fft(xc) #做傅立葉轉換\n \n # 整個segment的平均值、標準差、最大值\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 \n # 傅立葉轉換後，有分為實部和虛部\n realFFT = np.real(zc)\n imagFFT = np.imag(zc)\n \n # 實部平方加虛部平方\n X.loc[seg_id, 'RIsquare'] = np.mean(np.square(realFFT) + np.square(imagFFT))\n X.loc[seg_id, 'RImsquare'] = np.square(realFFT.mean()) + np.square(imagFFT.mean())\n \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 # 虛部的平均、標準差、最大最小值\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 # 實部的倒數5000個點的平均、標準差、最大最小值\n X.loc[seg_id, 'Rmean_last_5000'] = realFFT[-5000:].mean()\n X.loc[seg_id, 'Rstd__last_5000'] = realFFT[-5000:].std()\n X.loc[seg_id, 'Rmax_last_5000'] = realFFT[-5000:].max()\n X.loc[seg_id, 'Rmin_last_5000'] = realFFT[-5000:].min()\n \n # 實部的倒數15000個點的平均、標準差、最大最小值\n X.loc[seg_id, 'Rmean_last_15000'] = realFFT[-15000:].mean()\n X.loc[seg_id, 'Rstd_last_15000'] = realFFT[-15000:].std()\n X.loc[seg_id, 'Rmax_last_15000'] = realFFT[-15000:].max()\n X.loc[seg_id, 'Rmin_last_15000'] = realFFT[-15000:].min()\n \n # 整個segment的平均變化(差分)、變化率\n X.loc[seg_id, 'mean_change_abs'] = np.mean(np.diff(xc))\n X.loc[seg_id, 'mean_change_rate'] = np.mean(np.nonzero((np.diff(xc) / xc[:-1]))[0])\n \n # 整個segment取絕對值後的最大最小值\n X.loc[seg_id, 'abs_max'] = np.abs(xc).max()\n X.loc[seg_id, 'abs_min'] = np.abs(xc).min()\n \n # 原始segmant的前50000、後50000、前10000、後10000的標準差\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_10000'] = xc[:10000].std()\n X.loc[seg_id, 'std_last_10000'] = xc[-10000:].std()\n # 原始segmant的前50000、後50000、前10000、後10000的平均值\n X.loc[seg_id, 'avg_first_50000'] = xc[:50000].mean()\n X.loc[seg_id, 'avg_last_50000'] = xc[-50000:].mean()\n X.loc[seg_id, 'avg_first_10000'] = xc[:10000].mean()\n X.loc[seg_id, 'avg_last_10000'] = xc[-10000:].mean()\n # 原始segmant的前50000、後50000、前10000、後10000的最小值\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_10000'] = xc[:10000].min()\n X.loc[seg_id, 'min_last_10000'] = xc[-10000:].min()\n # 原始segmant的前50000、後50000、前10000、後10000的最大值\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_10000'] = xc[:10000].max()\n X.loc[seg_id, 'max_last_10000'] = xc[-10000:].max()\n\n #X.loc[seg_id, 'max_to_min'] = xc.max() / np.abs(xc.min())\n \n # 原始segment取絕對值後的全距\n X.loc[seg_id, 'max_to_min_diff'] = xc.max() - np.abs(xc.min())\n \n # 計算絕對值大於500的個數\n X.loc[seg_id, 'count_big'] = len(xc[np.abs(xc) > 500])\n \n # 原始segment的加總\n X.loc[seg_id, 'sum'] = xc.sum()\n \n # 前50000、後50000、前10000、後100000筆的平均變化率\n X.loc[seg_id, 'mean_change_rate_first_50000'] = np.mean(np.nonzero((np.diff(xc[:50000]) / xc[:50000][:-1]))[0])\n X.loc[seg_id, 'mean_change_rate_last_50000'] = np.mean(np.nonzero((np.diff(xc[-50000:]) / xc[-50000:][:-1]))[0])\n X.loc[seg_id, 'mean_change_rate_first_10000'] = np.mean(np.nonzero((np.diff(xc[:10000]) / xc[:10000][:-1]))[0])\n X.loc[seg_id, 'mean_change_rate_last_10000'] = np.mean(np.nonzero((np.diff(xc[-10000:]) / xc[-10000:][:-1]))[0])\n \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 X.loc[seg_id, 'q999'] = np.quantile(xc,0.999)\n X.loc[seg_id, 'q001'] = np.quantile(xc,0.001)\n \n # 絕對值的分位數\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 \n X.loc[seg_id, 'trend'] = add_trend_feature(xc)\n X.loc[seg_id, 'abs_trend'] = add_trend_feature(xc, abs_values=True)\n \n # 絕對值的平均值、標準差\n X.loc[seg_id, 'abs_mean'] = np.abs(xc).mean()\n X.loc[seg_id, 'abs_std'] = np.abs(xc).std()\n \n # 原始segmant的 Mean Absolute Deviation (平均絕對偏差)\n X.loc[seg_id, 'mad'] = xc.mad() \n \n # 計算不同window大小的移動平均值\n X.loc[seg_id, 'Moving_average_700_mean'] = xc.rolling(window=700).mean().mean(skipna=True)\n X.loc[seg_id, 'Moving_average_1500_mean'] = xc.rolling(window=1500).mean().mean(skipna=True)\n X.loc[seg_id, 'Moving_average_3000_mean'] = xc.rolling(window=3000).mean().mean(skipna=True)\n X.loc[seg_id, 'Moving_average_6000_mean'] = xc.rolling(window=6000).mean().mean(skipna=True)\n \n # EWM(exponential weighted functions)\n ewma = pd.Series.ewm\n X.loc[seg_id, 'exp_Moving_average_300_mean'] = (ewma(xc, span=300).mean()).mean(skipna=True)\n X.loc[seg_id, 'exp_Moving_average_3000_mean'] = ewma(xc, span=3000).mean().mean(skipna=True)\n X.loc[seg_id, 'exp_Moving_average_30000_mean'] = ewma(xc, span=6000).mean().mean(skipna=True)\n \n no_of_std = 2\n X.loc[seg_id, 'MA_700MA_std_mean'] = xc.rolling(window=700).std().mean()\n X.loc[seg_id,'MA_700MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n X.loc[seg_id,'MA_700MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_700MA_std_mean']).mean()\n X.loc[seg_id, 'MA_400MA_std_mean'] = xc.rolling(window=400).std().mean()\n X.loc[seg_id,'MA_400MA_BB_high_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] + no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n X.loc[seg_id,'MA_400MA_BB_low_mean'] = (X.loc[seg_id, 'Moving_average_700_mean'] - no_of_std * X.loc[seg_id, 'MA_400MA_std_mean']).mean()\n X.loc[seg_id, 'MA_1000MA_std_mean'] = xc.rolling(window=1000).std().mean()\n \n # 四分位距 IQR(interquartile range) \n X.loc[seg_id, 'iqr'] = np.subtract(*np.percentile(xc, [75, 25]))\n \n # 減去分布兩端10%的值後，計算平均數\n X.loc[seg_id, 'ave10'] = stats.trim_mean(xc, 0.1)\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 X.loc[seg_id, 'av_change_abs_roll_mean_' + str(windows)] = np.mean(np.diff(x_roll_mean))\n X.loc[seg_id, 'av_change_rate_roll_mean_' + str(windows)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n X.loc[seg_id, 'abs_max_roll_mean_' + str(windows)] = np.abs(x_roll_mean).max() ","execution_count":38,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# iterate over all segments\nfor seg_id in tqdm_notebook(range(segments)):\n seg = train_df.iloc[seg_id*int(rows/8):seg_id*int(rows/8)+rows]\n create_features(seg_id, seg, train_X)\n # the y value is the last entry in the time to failure in the segment\n train_y.loc[seg_id, 'time_to_failure'] = seg['time_to_failure'].values[-1]","execution_count":39,"outputs":[{"output_type":"display_data","data":{"text/plain":"HBox(children=(IntProgress(value=0, max=33554), HTML(value='')))","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"678814fe3a084716ac464119108c335d"}},"metadata":{}},{"output_type":"error","ename":"KeyboardInterrupt","evalue":"","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)","\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mseg_id\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtqdm_notebook\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msegments\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mseg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtrain_df\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0miloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrows\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrows\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mrows\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mcreate_features\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mseg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrain_X\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0;31m# the y value is the last entry in the time to failure in the segment\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mtrain_y\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'time_to_failure'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mseg\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'time_to_failure'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m\u001b[0m in \u001b[0;36mcreate_features\u001b[0;34m(seg_id, seg, X)\u001b[0m\n\u001b[1;32m 119\u001b[0m \u001b[0mewma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSeries\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mewm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 120\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'exp_Moving_average_300_mean'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mewma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mspan\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m300\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mskipna\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 121\u001b[0;31m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'exp_Moving_average_3000_mean'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mewma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mspan\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m3000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mskipna\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 122\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mseg_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'exp_Moving_average_30000_mean'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mewma\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mspan\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m6000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mskipna\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 123\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.6/site-packages/pandas/core/window.py\u001b[0m in \u001b[0;36mmean\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 2176\u001b[0m \u001b[0;34m\"\"\"exponential weighted moving average\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2177\u001b[0m \u001b[0mnv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalidate_window_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'mean'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2178\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_apply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ewma'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2179\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2180\u001b[0m \u001b[0;34m@\u001b[0m\u001b[0mSubstitution\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'ewm'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.6/site-packages/pandas/core/window.py\u001b[0m in \u001b[0;36m_apply\u001b[0;34m(self, func, **kwargs)\u001b[0m\n\u001b[1;32m 2167\u001b[0m int(self.ignore_na), int(self.min_periods))\n\u001b[1;32m 2168\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2169\u001b[0;31m \u001b[0mresults\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply_along_axis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2170\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2171\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_wrap_results\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresults\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mblocks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.6/site-packages/numpy/lib/shape_base.py\u001b[0m in \u001b[0;36mapply_along_axis\u001b[0;34m(func1d, axis, arr, *args, **kwargs)\u001b[0m\n\u001b[1;32m 378\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mStopIteration\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 379\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Cannot apply_along_axis when any iteration dimensions are 0'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 380\u001b[0;31m \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masanyarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc1d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minarr_view\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mind0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 381\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 382\u001b[0m \u001b[0;31m# build a buffer for storing evaluations of func1d.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;32m/opt/conda/lib/python3.6/site-packages/pandas/core/window.py\u001b[0m in \u001b[0;36mfunc\u001b[0;34m(arg)\u001b[0m\n\u001b[1;32m 2165\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2166\u001b[0m return cfunc(arg, self.com, int(self.adjust),\n\u001b[0;32m-> 2167\u001b[0;31m int(self.ignore_na), int(self.min_periods))\n\u001b[0m\u001b[1;32m 2168\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2169\u001b[0m \u001b[0mresults\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply_along_axis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalues\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mKeyboardInterrupt\u001b[0m: "]}]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X_save = train_X.copy\ntrain_y_save = train_y.copy\ntrain_y.head()\n#experiment with variations of time\ntrain_y = train_y**1.0\ntrain_X.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We will not train on the segments with a quake, because there are likely outliers\n# np.diff的功能是得到每一項減掉前一項的值，'> 0'使得發生quake的地方是True，其餘是False\n# np.nonzero取得每項非零的值的index，如下(array([ 36, 331, 694, 921, 1245, 1451, 1631, 2044, 2246, 2492, 2784, 3066, 3292, 3511, 3888, 4130]),)\ntrain_y_quake = np.nonzero(np.diff(train_y.time_to_failure) > 0)[0] + 1\nprint(len(train_y_quake))\nprint (len(train_y))\n\nfor idx in train_y_quake: \n train_y.drop([idx],inplace=True)\n train_X.drop([idx],inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print (len(train_y))\ntrain_y.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X.shape, train_y.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X.head(), train_y.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_X.to_csv('train_X.csv',index=False)\ntrain_y.to_csv('train_y.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"對訓練資料進行尺度縮放的處理"},{"metadata":{"trusted":true},"cell_type":"code","source":"scaler = StandardScaler()\nscaler.fit(train_X)\nscaled_train_X = pd.DataFrame(scaler.transform(train_X), columns=train_X.columns)\n#scaled_train_X = train_X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_train_X.head(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"temp = train_X.copy","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 處理測試資料\n對於測試資料做相同的處理"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('../input/sample_submission.csv', index_col='seg_id')\ntest_X = pd.DataFrame(columns=train_X.columns, dtype=np.float64, index=submission.index)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.shape, test_X.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"對測試資料取同樣的feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"for seg_id in tqdm_notebook(test_X.index):\n seg = pd.read_csv('../input/test/' + seg_id + '.csv')\n # convert to the log of the absolute values\n seg.acoustic_data =np.log1p(abs(seg.acoustic_data))\n create_features(seg_id, seg, test_X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_X.to_csv('test_X.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"對測試資料也進行使度縮放處理"},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_test_X = pd.DataFrame(scaler.transform(test_X), columns=test_X.columns)\n#scaled_test_X = test_X\nscaled_test_X.values[1117]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_test_X.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"scaled_test_X.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 模型部份\n先定義一個驗證的function，用以評估模型效能"},{"metadata":{},"cell_type":"markdown","source":"## 交叉驗證 Function"},{"metadata":{"trusted":true},"cell_type":"code","source":"n_fold = 5\ndef mae_cv (model):\n folds = KFold(n_splits=n_fold, shuffle=True, random_state=42).get_n_splits(scaled_train_X.values)\n mae = -cross_val_score (model, scaled_train_X.values, train_y, scoring=\"neg_mean_absolute_error\",\n verbose=0,\n cv=folds)\n return mae","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 開始定義各種不同模型，用以Ensemble"},{"metadata":{},"cell_type":"markdown","source":"## LGB Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''%%time\n\nlgb_params = {'num_leaves': 51,\n 'min_data_in_leaf': 10, \n 'objective':'regression',\n 'max_depth': -1,\n 'learning_rate': 0.001,\n \"boosting\": \"gbdt\",\n \"feature_fraction\": 0.91,\n \"bagging_freq\": 1,\n \"bagging_fraction\": 0.91,\n \"bagging_seed\": 42,\n \"metric\": 'mae',\n \"lambda_l1\": 0.1,\n \"verbosity\": -1,\n \"nthread\": -1,\n \"random_state\": 42}\n\n\nlgb_model = lgb.LGBMRegressor(objective='regression',num_leaves=5,\n learning_rate=0.01, n_estimators=720,\n bagging_freq = 5, feature_fraction = 0.2319,\n feature_fraction_seed=9, bagging_seed=9,\n min_data_in_leaf =6, min_sum_hessian_in_leaf = 11, n_jobs = -1)\n\n#score = mae_cv(lgb_model)\n#print(\"LGBM score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nlgb_model'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## XGB Model"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''%%time \nxgb_params = {'eta': 0.03,\n 'max_depth': 9,\n 'subsample': 0.85,\n 'objective': 'reg:linear',\n 'eval_metric': 'mae',\n 'silent': True,\n 'nthread': 4}\n \nxgb_model = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, \n learning_rate=0.05, max_depth=3, \n min_child_weight=1.7817, n_estimators=2200,\n reg_alpha=0.4640, reg_lambda=0.8571,\n subsample=0.5213, silent=1,\n random_state =7, nthread = -1, eval_metric = 'mae',)\n\n#score = mae_cv(xgb_model)\n#print(\"XGB score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nxgb_model\n\n# xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, \n# verbose_eval=500, params=xgb_params)'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Random Forest"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nrf_model = RandomForestRegressor(n_estimators=120, n_jobs=-1, min_samples_leaf=1, \n max_features = \"auto\",max_depth=15, )\n#score = mae_cv(rf_model)\n#print(\"Random Forest score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nrf_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Cat Boost"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\nparams = {'loss_function':'MAE',}\ncat_model = CatBoostRegressor(iterations=1000, eval_metric='MAE', verbose=False, **params)\n\n#score = mae_cv(cat_model)\n#print(\"Cat Boost score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\ncat_model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Kernel Ridge"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''%%time\nKRR_model = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)\n#score = mae_cv(KRR_model)\n#print(\"Kernel Ridge score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\n#print (score)\nKRR_model'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Elastic Net"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''%%time\n#ENet_model = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=0.9, random_state=3,max_iter=5000))\nENet_model = ElasticNet(alpha=0.0005, l1_ratio=0.9, random_state=3,max_iter=5000)\n#score = mae_cv(ENet_model)\n#print(\"Elastic Net score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nENet_model'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Lasso"},{"metadata":{"trusted":true},"cell_type":"code","source":"'''%%time\nlasso_model = Lasso(alpha =0.0005, random_state=1)\n#score = mae_cv(lasso_model)\n#print(\"Lasso score: {:.4f} ({:.4f})\\n\" .format(score.mean(), score.std()))\nlasso_model'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"trusted":true},"cell_type":"code","source":"class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):\n def __init__(self, models):\n self.models = models\n \n # we define clones of the original models to fit the data in\n def fit(self, X, y):\n self.models_ = [clone(x) for x in self.models]\n \n # Train cloned base models\n for model in self.models_:\n model.fit(X, y)\n\n return self\n \n #Now we do the predictions for cloned models and average them\n def predict(self, X):\n \n predictions = np.column_stack([\n model.predict(X) for model in self.models_\n ])\n return np.mean(predictions, axis=1) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):\n def __init__(self, base_models, meta_model, n_folds=5):\n self.base_models = base_models\n self.meta_model = meta_model\n self.n_folds = n_folds\n \n # We again fit the data on clones of the original models\n def fit(self, X, y):\n print (type(X))\n self.base_models_ = [list() for x in self.base_models]\n self.meta_model_ = clone(self.meta_model)\n kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)\n print (KFold)\n # Train cloned base models then create out-of-fold predictions\n # that are needed to train the cloned meta-model\n out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))\n for i, model in enumerate(self.base_models):\n for train_index, holdout_index in kfold.split(X, y):\n instance = clone(model)\n self.base_models_[i].append(instance)\n instance.fit(X.iloc[train_index], y.iloc[train_index])\n y_pred = instance.predict(X.iloc[holdout_index])\n out_of_fold_predictions[holdout_index, i] = y_pred\n \n # Now train the cloned meta-model using the out-of-fold predictions as new feature\n self.meta_model_.fit(out_of_fold_predictions, y)\n return self\n \n #Do the predictions of all base models on the test data and use the averaged predictions as \n #meta-features for the final prediction which is done by the meta-model\n def predict(self, X):\n meta_features = np.column_stack([\n np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)\n for base_models in self.base_models_ ])\n return self.meta_model_.predict(meta_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class StackingCVRegressorRetrained(BaseEstimator, RegressorMixin, TransformerMixin):\n def __init__(self, regressors, meta_regressor, n_folds=5, use_features_in_secondary=False):\n self.regressors = regressors\n self.meta_regressor = meta_regressor\n self.n_folds = n_folds\n self.use_features_in_secondary = use_features_in_secondary\n\n def fit(self, X, y):\n self.regr_ = [clone(x) for x in self.regressors]\n self.meta_regr_ = clone(self.meta_regressor)\n\n kfold = KFold(n_splits=self.n_folds, shuffle=True)\n\n out_of_fold_predictions = np.zeros((X.shape[0], len(self.regressors)))\n\n # Create out-of-fold predictions for training meta-model\n for i, regr in enumerate(self.regr_):\n for train_idx, holdout_idx in kfold.split(X, y):\n instance = clone(regr)\n instance.fit(X[train_idx], y[train_idx])\n out_of_fold_predictions[holdout_idx, i] = instance.predict(X[holdout_idx])\n\n # Train meta-model\n if self.use_features_in_secondary:\n self.meta_regr_.fit(np.hstack((X, out_of_fold_predictions)), y)\n else:\n self.meta_regr_.fit(out_of_fold_predictions, y)\n \n # Retrain base models on all data\n for regr in self.regr_:\n regr.fit(X, y)\n\n return self\n\n def predict(self, X):\n meta_features = np.column_stack([\n regr.predict(X) for regr in self.regr_\n ])\n\n if self.use_features_in_secondary:\n return self.meta_regr_.predict(np.hstack((X, meta_features)))\n else:\n return self.meta_regr_.predict(meta_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n#averaged_models = AveragingModels(models = (rf_model, xgb_model, KRR_model, lgb_model, ENet_model, cat_model, lasso_model))\n#averaged_models = AveragingModels(models = (rf_model, lgb_model, cat_model, lasso_model))\naveraged_models = AveragingModels(models = (rf_model,cat_model))\n\nscore = mae_cv(averaged_models)\nprint(\" Averaged base models score: {:.4f} ({:.4f})\\n\".format(score.mean(), score.std()))\n#averaged_models.fit (scaled_train_X.values, train_y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"averaged_models.fit (scaled_train_X.values, train_y)\naveraged_train_predict = averaged_models.predict(scaled_train_X.values)\nprint(mean_absolute_error(train_y, averaged_train_predict))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"averaged_prediction = np.zeros(len(scaled_test_X))\naveraged_prediction += averaged_models.predict(scaled_test_X.values)\naveraged_prediction","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n#stacked_predict = StackingAveragedModels(base_models =(rf_model, xgb_model, lgb_model, cat_model,ENet_model), \n# meta_model =lasso_model) \n#stacked_predict.fit(scaled_train_X, train_y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#stacked_train_pred = stacked_predict.predict(scaled_train_X)\n\n#print(mean_absolute_error(train_y, stacked_train_pred))\n\n#stacked_prediction = np.zeros(len(scaled_test_X))\n#stacked_prediction += stacked_predict.predict(scaled_test_X)**1.0\n#stacked_prediction[0:4]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"submission.time_to_failure = averaged_prediction\nsubmission.to_csv('submission_averaged_cat_rf_8.csv',index=True)\n#submission.time_to_failure = stacked_prediction\n#submission.to_csv('submissionV28_stacked.csv',index=True)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}