今天給你們介紹提高方法(Boosting), 提高算法是一種能夠用來減少監督式學習中誤差的機器學習算法。算法
面對的問題是邁可·肯斯（Michael Kearns）提出的：一組「弱學習者」的集合可否生成一個「強學習者」？json
弱學習者通常是指一個分類器，它的結果只比隨機分類好一點點。強學習者指分類器的結果很是接近真值。app
大多數提高算法包括由迭代使用弱學習分類器組成，並將其結果加入一個最終的成強學習分類器。加入的過程當中，一般根據它們的分類準確率給予不一樣的權重。加和弱學習者以後，數據一般會被從新加權，來強化對以前分類錯誤數據點的分類。dom
提高算法有種三個臭皮匠頂個諸葛亮的意思。在這裏將使用微軟的LightGBM這個提高算法，來預測分子的一個屬性，叫作耦合常數。機器學習
導入須要的庫
import os
import time
import datetime
import json
import gc
from numba import jit
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from tqdm import tqdm_notebook
from sklearn.preprocessing import StandardScaler
from sklearn.svm import NuSVR, SVR
from sklearn.metrics import mean_absolute_error
pd.options.display.precision = 15
import lightgbm as lgb
import xgboost as xgb
from catboost import CatBoostRegressor
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold
from sklearn import metrics
from sklearn import linear_model
import seaborn as sns
import warnings
warnings.filterwarnings("ignore")
from IPython.display import HTML
import altair as alt
from altair.vega import v5
from IPython.display import HTML
import networkx as nx
import matplotlib.pyplot as plt
%matplotlib inline
alt.renderers.enable('notebook')
讀取數據
讀取須要使用的csv文件，函數
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
structures = pd.read_csv('structures.csv')
csv文件中包含了如下這些數據，首先看下訓練集的數據，從pandas顯示的數據來看，訓練集的csv包含了，分子名，原子index，化學鍵類型和訓練集給出的耦合常數。學習
train.head()
在測試集數據中，除了耦合常數的數據其餘的數據跟訓練集的數據是一致的。測試
test.head()
再來看下structures的數據，也就是原子的數據，其中包含了，原子的索引，這個原子在哪一個分子中，和它的空間座標位置，這些數據後續將會用來作特徵挖掘。編碼
structures.head()
可數據化的數據
數據可視化一直是機器學習的重點，一個好的數據可視化，能夠幫助咱們理解數據的分佈、形態以及特性。atom
首先來看下耦合常數在已有的數據集中分佈，這裏展現了耦合常數的直方圖，以及每一個化學鍵類型下的耦合常數的數據大小的分佈狀況。
fig, ax = plt.subplots(figsize = (18, 6))
plt.subplot(1, 2, 1);
plt.hist(train['scalar_coupling_constant'], bins=20);
plt.title('Basic scalar_coupling_constant histogram');
plt.subplot(1, 2, 2);
sns.violinplot(x='type', y='scalar_coupling_constant', data=train);
plt.title('Violinplot of scalar_coupling_constant by type');
再來看下已有的數據，根據不一樣的化學鍵類型的原子鏈接形態。
從圖表中能夠發現，不一樣給的化學鍵類型，有着不一樣的形態，其中2JHH這個化學鍵的形態最爲特別。
fig, ax = plt.subplots(figsize = (20, 12))
for i, t in enumerate(train['type'].unique()):
train_type = train.loc[train['type'] == t]
G = nx.from_pandas_edgelist(train_type, 'atom_index_0', 'atom_index_1', ['scalar_coupling_constant'])
plt.subplot(2, 4, i + 1);
nx.draw(G, with_labels=True);
plt.title(f'Graph for type {t}')
特徵工程
提高算法，通常的數據形式都是基於表單的數據。好比說此次的案例中分子與原子的數據。
對於這類數據，最重要的一個工做就是特徵工程，如何從現有的特徵中挖掘出其餘有價值的特徵將是這類算法的最大的難點。下面咱們作一些簡單的特徵工程，來展現下特徵工程的過程，
首先，將結構表跟訓練數據表和測試數據表進行合併。
def map_atom_info(df, atom_idx):
df = pd.merge(df, structures, how = 'left',
left_on = ['molecule_name', f'atom_index_{atom_idx}'],
right_on = ['molecule_name', 'atom_index'])
df = df.drop('atom_index', axis=1)
df = df.rename(columns={'atom': f'atom_{atom_idx}',
'x': f'x_{atom_idx}',
'y': f'y_{atom_idx}',
'z': f'z_{atom_idx}'})
return df
train = map_atom_info(train, 0)
train = map_atom_info(train, 1)
test = map_atom_info(test, 0)
test = map_atom_info(test, 1)
觀察合併之後的數據，發現原子的類型以及空間位置都合併到了訓練集數據和測試集數據相應的行中。
train.head()
test.head()
計算距離特徵
分別計算訓練集和測試集數據集中，每一行中兩個原子間的距離特徵，其中包括，
-
原子距離；
-
原子x上的距離；
-
原子y上的距離；
-
原子z上的距離；
而且把這些距離合併到訓練集和測試集數據集上，
train_p_0 = train[['x_0', 'y_0', 'z_0']].values train_p_1 = train[['x_1', 'y_1', 'z_1']].values test_p_0 = test[['x_0', 'y_0', 'z_0']].values test_p_1 = test[['x_1', 'y_1', 'z_1']].values train['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1) test['dist'] = np.linalg.norm(test_p_0 - test_p_1, axis=1) train['dist_x'] = (train['x_0'] - train['x_1']) ** 2 test['dist_x'] = (test['x_0'] - test['x_1']) ** 2 train['dist_y'] = (train['y_0'] - train['y_1']) ** 2 test['dist_y'] = (test['y_0'] - test['y_1']) ** 2 train['dist_z'] = (train['z_0'] - train['z_1']) ** 2 test['dist_z'] = (test['z_0'] - test['z_1']) ** 2
計算每種化學鍵的距離的均值，
train['type_0'] = train['type'].apply(lambda x: x[0]) test['type_0'] = test['type'].apply(lambda x: x[0]) train['type_1'] = train['type'].apply(lambda x: x[1:]) test['type_1'] = test['type'].apply(lambda x: x[1:])
train['type_0'] = train['type'].apply(lambda x: x[0])
test['type_0'] = test['type'].apply(lambda x: x[0])
train['type_1'] = train['type'].apply(lambda x: x[1:])
test['type_1'] = test['type'].apply(lambda x: x[1:])
train['dist_to_type_mean'] = train['dist'] / train.groupby('type')['dist'].transform('mean')
test['dist_to_type_mean'] = test['dist'] / test.groupby('type')['dist'].transform('mean')
train['dist_to_type_0_mean'] = train['dist'] / train.groupby('type_0')['dist'].transform('mean')
test['dist_to_type_0_mean'] = test['dist'] / test.groupby('type_0')['dist'].transform('mean')
train['dist_to_type_1_mean'] = train['dist'] / train.groupby('type_1')['dist'].transform('mean')
test['dist_to_type_1_mean'] = test['dist'] / test.groupby('type_1')['dist'].transform('mean')
將化學鍵以分子爲單位進行分組並計算分子均值，
train[f'molecule_type_dist_mean'] = train.groupby(['molecule_name', 'type'])['dist'].transform('mean')
test[f'molecule_type_dist_mean'] = test.groupby(['molecule_name', 'type'])['dist'].transform('mean')
將原子與化學鍵類型編碼成數值，方便後續訓練。
for f in ['atom_0', 'atom_1', 'type_0', 'type_1', 'type']:
lbl = LabelEncoder()
lbl.fit(list(train[f].values) + list(test[f].values))
train[f] = lbl.transform(list(train[f].values))
test[f] = lbl.transform(list(test[f].values))
準備數據
首先排除一些咱們不須要的列數據，
X = train.drop(['id', 'molecule_name', 'scalar_coupling_constant'], axis=1) y = train['scalar_coupling_constant'] X_test = test.drop(['id', 'molecule_name'], axis=1)
而後將數據分拆成幾個fold，每一個fold的數據都會打亂。
n_fold = 5 folds = KFold(n_splits=n_fold, shuffle=True, random_state=11)
定義模型
特徵工程完成以後，就須要定義咱們的模型，這裏咱們使用LightGBM來訓練。
LigthGBM是boosting集合模型中的新進成員，由微軟提供，它和XGBoost同樣是對GBDT的高效實現，原理上它和GBDT及XGBoost相似，都採用損失函數的負梯度做爲當前決策樹的殘差近似值，去擬合新的決策樹。
LightGBM在不少方面會比XGBoost表現的更爲優秀，好比：更快的訓練效率；低內存使用；更高的準確率；支持並行化學習可處理大規模數據；支持直接使用category特徵等。
首先咱們須要設置LightGBM的參數
LightGBM的參數不少，在這裏對某幾個關鍵的參數作下解釋
max_depth：樹模型深度；
num_leaves：葉子節點數，數模型複雜度；
min_child_samples：表示一個葉子節點上包含的最少樣本數量；
objective：目標函數；
learning_rate：學習率；
max_depth：控制了樹的最大深度。該參數能夠顯式的限制樹的深度；
boosting_type：給出了基學習器模型算法；
params = {'num_leaves': 128,
'min_child_samples': 79,
'objective': 'regression',
'max_depth': 13,
'learning_rate': 0.2,
"boosting_type": "gbdt",
"subsample_freq": 1,
"subsample": 0.9,
"bagging_seed": 11,
"metric": 'mae',
"verbosity": -1,
'reg_alpha': 0.1,
'reg_lambda': 0.3,
'colsample_bytree': 1.0
}
定義訓練函數, 在訓練函數裏記錄了每一個特徵的重要性，並在訓練結束後經過圖表展現，有助於篩選特徵。
def group_mean_log_mae(y_true, y_pred, types, floor=1e-9):
maes = (y_true-y_pred).abs().groupby(types).mean()
return np.log(maes.map(lambda x: max(x, floor))).mean()
def train_model_regression(X, X_test, y, params, folds, model_type='lgb', eval_metric='mae', columns=None, plot_feature_importance=False, model=None,
verbose=10000, early_stopping_rounds=200, n_estimators=50000):
columns = X.columns if columns is None else columns
X_test = X_test[columns]
# to set up scoring parameters
metrics_dict = {'mae': {'lgb_metric_name': 'mae',
'catboost_metric_name': 'MAE',
'sklearn_scoring_function': metrics.mean_absolute_error},
'group_mae': {'lgb_metric_name': 'mae',
'catboost_metric_name': 'MAE',
'scoring_function': group_mean_log_mae},
'mse': {'lgb_metric_name': 'mse',
'catboost_metric_name': 'MSE',
'sklearn_scoring_function': metrics.mean_squared_error}
}
result_dict = {}
# out-of-fold predictions on train data
oof = np.zeros(len(X))
# averaged predictions on train data
prediction = np.zeros(len(X_test))
# list of scores on folds
scores = []
feature_importance = pd.DataFrame()
# split and train on folds
for fold_n, (train_index, valid_index) in enumerate(folds.split(X)):
print(f'Fold {fold_n + 1} started at {time.ctime()}')
if type(X) == np.ndarray:
X_train, X_valid = X[columns][train_index], X[columns][valid_index]
y_train, y_valid = y[train_index], y[valid_index]
else:
X_train, X_valid = X[columns].iloc[train_index], X[columns].iloc[valid_index]
y_train, y_valid = y.iloc[train_index], y.iloc[valid_index]
if model_type == 'lgb':
model = lgb.LGBMRegressor(**params, n_estimators = n_estimators, n_jobs = -1)
model.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_valid, y_valid)], eval_metric=metrics_dict[eval_metric]['lgb_metric_name'],
verbose=verbose, early_stopping_rounds=early_stopping_rounds)
y_pred_valid = model.predict(X_valid)
y_pred = model.predict(X_test, num_iteration=model.best_iteration_)
if model_type == 'xgb':
train_data = xgb.DMatrix(data=X_train, label=y_train, feature_names=X.columns)
valid_data = xgb.DMatrix(data=X_valid, label=y_valid, feature_names=X.columns)
watchlist = [(train_data, 'train'), (valid_data, 'valid_data')]
model = xgb.train(dtrain=train_data, num_boost_round=20000, evals=watchlist, early_stopping_rounds=200, verbose_eval=verbose, params=params)
y_pred_valid = model.predict(xgb.DMatrix(X_valid, feature_names=X.columns), ntree_limit=model.best_ntree_limit)
y_pred = model.predict(xgb.DMatrix(X_test, feature_names=X.columns), ntree_limit=model.best_ntree_limit)
if model_type == 'sklearn':
model = model
model.fit(X_train, y_train)
y_pred_valid = model.predict(X_valid).reshape(-1,)
score = metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid)
print(f'Fold {fold_n}. {eval_metric}: {score:.4f}.')
print('')
y_pred = model.predict(X_test).reshape(-1,)
if model_type == 'cat':
model = CatBoostRegressor(iterations=20000, eval_metric=metrics_dict[eval_metric]['catboost_metric_name'], **params,
loss_function=metrics_dict[eval_metric]['catboost_metric_name'])
model.fit(X_train, y_train, eval_set=(X_valid, y_valid), cat_features=[], use_best_model=True, verbose=False)
y_pred_valid = model.predict(X_valid)
y_pred = model.predict(X_test)
oof[valid_index] = y_pred_valid.reshape(-1,)
if eval_metric != 'group_mae':
scores.append(metrics_dict[eval_metric]['sklearn_scoring_function'](y_valid, y_pred_valid))
else:
scores.append(metrics_dict[eval_metric]['scoring_function'](y_valid, y_pred_valid, X_valid['type']))
prediction += y_pred
if model_type == 'lgb' and plot_feature_importance:
# feature importance
fold_importance = pd.DataFrame()
fold_importance["feature"] = columns
fold_importance["importance"] = model.feature_importances_
fold_importance["fold"] = fold_n + 1
feature_importance = pd.concat([feature_importance, fold_importance], axis=0)
prediction /= folds.n_splits
print('CV mean score: {0:.4f}, std: {1:.4f}.'.format(np.mean(scores), np.std(scores)))
result_dict['oof'] = oof
result_dict['prediction'] = prediction
result_dict['scores'] = scores
if model_type == 'lgb':
if plot_feature_importance:
feature_importance["importance"] /= folds.n_splits
cols = feature_importance[["feature", "importance"]].groupby("feature").mean().sort_values(
by="importance", ascending=False)[:50].index
best_features = feature_importance.loc[feature_importance.feature.isin(cols)]
plt.figure(figsize=(16, 12));
sns.barplot(x="importance", y="feature", data=best_features.sort_values(by="importance", ascending=False));
plt.title('LGB Features (avg over folds)');
result_dict['feature_importance'] = feature_importance
return result_dict
開始訓練
在訓練結束的時候，會輸出每一個特徵值的重要程度，柱狀圖越長，代表這個特徵在預測時候所起到的做用越大。
經過重要性的觀察，能夠幫助咱們篩選出較好的特徵，剔除一些不過重要的特徵。一來特徵越少訓練速度越快，二來保留重要的特徵不會對結果有較大影響。
result_dict_lgb = train_model_regression(X=X, X_test=X_test, y=y, params=params, folds=folds, model_type='lgb', eval_metric='group_mae', plot_feature_importance=True,
verbose=1000, early_stopping_rounds=200, n_estimators=10000)
預測結果分佈
圖表中顯示了每一個化學鍵類型的預測結果分佈狀況，數據越集中，說明咱們的預測結果越一致，相比來講，效果也越好。
plot_data = pd.DataFrame(y)
plot_data.index.name = 'id'
plot_data['yhat'] = result_dict_lgb['oof']
plot_data['type'] = lbl.inverse_transform(X['type'])
def plot_oof_preds(ctype, llim, ulim):
plt.figure(figsize=(6,6))
sns.scatterplot(x='scalar_coupling_constant',y='yhat',
data=plot_data.loc[plot_data['type']==ctype,
['scalar_coupling_constant', 'yhat']]);
plt.xlim((llim, ulim))
plt.ylim((llim, ulim))
plt.plot([llim, ulim], [llim, ulim])
plt.xlabel('scalar_coupling_constant')
plt.ylabel('predicted')
plt.title(f'{ctype}', fontsize=18)
plt.show()
plot_oof_preds('1JHC', 0, 250)
plot_oof_preds('1JHN', 0, 100)
plot_oof_preds('2JHC', -50, 50)
plot_oof_preds('2JHH', -50, 50)
plot_oof_preds('2JHN', -25, 25)
plot_oof_preds('3JHC', -25, 100)
plot_oof_preds('3JHH', -20, 20)
plot_oof_preds('3JHN', -15, 15)
「預測分子屬性」案例將在矩池雲上線，感興趣的用戶能夠在矩池雲Demo鏡像中體驗。