15. Machine Learning with Scikit¶

# 生成 10x2 大小的 X，以及對應的 10 個 y 數列
X = np.arange(20).reshape((10, 2)) 
y = np.zeros(10)
y[6:] = 1.0

print('X =\n', X)
print('y =\n', y)

from sklearn.model_selection import train_test_split

# 刻意指定一個很不剛好的比例
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)

# 指定的比例如果沒有剛好整數，切割時還是會調整成整數
print('X_train =\n', X_train)
print('y_train =\n\n', y_train)
print('X_test =\n', X_test)
print('y_test =\n', y_test)

# 沒有 shuffle 的話，切割是照原本順序
train_test_split(y, test_size=0.3, shuffle=False)

# 指定訓練集測試集的類別分布要一致
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y)

print('X_train =\n', X_train)
print('y_train =\n\n', y_train)
print('X_test =\n', X_test)
print('y_test =\n', y_test)

# 典型的模型訓練流程圖
Image(url='https://scikit-learn.org/stable/_images/grid_search_workflow.png', embed=True, width=480)

# k-fold 交叉驗證
Image(url='https://scikit-learn.org/stable/_images/grid_search_cross_validation.png', embed=True, width=480)

from sklearn.model_selection import KFold

kfold = KFold(n_splits=3)
# 使用之前生成的 10x2 大小的 X，以及 10 個 y 的數列當訓練集
for k, (train_index, test_index) in enumerate(kfold.split(X, y)):
    print('>> Fold #', k)
    X_train, X_test = X[train_index], X[test_index]
    y_train, y_test = y[train_index], y[test_index]
    print(' - Train index:', train_index, '\n - X:', X_train, '\n - y:', y_train, '\n')
    print(' - Test index:', test_index, '\n - X:', X_test, '\n - y:', y_test, '\n')

# 生成 100 個新的範例資料點
n_points = 100
X = np.random.randn(100, 10)

# y 包含三個分布非常不平衡的類別
percentiles_classes = [0.1, 0.3, 0.6]
y = np.hstack([[c] * int(100 * perc) for c, perc in enumerate(percentiles_classes)])

# 要顯示的不同類別
from sklearn.model_selection import (
    KFold, ShuffleSplit, StratifiedKFold, StratifiedShuffleSplit
)
cvs = [KFold, ShuffleSplit, StratifiedKFold, StratifiedShuffleSplit]
n_splits = 4

# 用圖顯示不同 k-fold 類別的分割分配
fig, axs = plt.subplots(2, 2, figsize=(12, 8))

cmap_data = plt.cm.Paired
cmap_cv = plt.cm.coolwarm
x_lim = [0, n_points]
y_lim = [n_splits+1.2, -0.2]
y_ticks = np.arange(n_splits+1) + 0.5
y_ticklabels = list(range(n_splits)) + ['class']

from matplotlib.patches import Patch

for n, (cv, ax) in enumerate(zip(cvs, axs.reshape(-1))):
    cv_n = cv(n_splits=n_splits)
    for k, (idxtrain, idxtest) in enumerate(cv_n.split(X=X, y=y)):
        indices = np.zeros(X.shape[0])
        indices[idxtest] = 1.0
        ax.scatter(range(X.shape[0]), [k+0.5] * X.shape[0],
                   c=indices, marker='_', lw=20, cmap=cmap_cv,
                   vmin=-0.2, vmax=1.2)

    # Plot the data classes at the end
    ax.scatter(range(X.shape[0]), [k + 1.5] * X.shape[0],
               c=y, marker='_', lw=20, cmap=cmap_data)

    ax.set(xlim=x_lim, ylim=y_lim, yticks=y_ticks, yticklabels=y_ticklabels)
    if n in [0, 2]:
        ax.set_ylabel('CV iteration', fontsize=14)
    if n in [2, 3]:
        ax.set_xlabel('Sample index', fontsize=14)

    ax.set_title(type(cv_n).__name__, fontsize=16, fontweight='bold')

    ax.legend([Patch(color=cmap_cv(.8)), Patch(color=cmap_cv(.02))],
              ['Testing set', 'Training set'], loc='lower right')

plt.tight_layout()
plt.show()

# 取得訓練集與測試集
dsWdbc = WdbcDataset()
X_train, X_test, Y_train, Y_test = dsWdbc.get_xy()

# 顯示切割後的大小
print('X: train shape = {}, test shape = {}'.format(X_train.shape, X_test.shape))
print('Y: train shape = {}, test shape = {}'.format(Y_train.shape, Y_test.shape))

# 顯示切割後的類別分布
train_labels, train_counts = np.unique(Y_train, return_counts=True)
print('Train Classes {}:{} = {}:{}'.format(train_labels[0], train_labels[1], train_counts[0], train_counts[1]))
test_labels, test_counts = np.unique(Y_test, return_counts=True)
print('Test Classes {}:{} = {}:{}'.format(test_labels[0], test_labels[1], test_counts[0], test_counts[1]))

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_validate

sklr1 = LogisticRegression(max_iter=100)
# 10-fold 交叉驗證，使用所有處理核心平行運算
scores1 = cross_validate(sklr1, X_train, Y_train, cv=10, n_jobs=-1)

# 顯示每個 fold 的 accuracy 分數
# 轉成 DataFrame 方便觀察，也方便輸出成檔案做記錄
dfscores1 = pd.DataFrame(scores1)
dfscores1

# 不同超參數訓練的模型比較，通常觀察平均分數
print('[Average scores]:\n', dfscores1.mean(axis=0))

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedShuffleSplit
from sklearn.preprocessing import StandardScaler

sklr2 = LogisticRegression(max_iter=100)
ss10fold = StratifiedShuffleSplit(n_splits=10)

scores2 = []
for k, (idx_train, idx_validate) in enumerate(ss10fold.split(X_train, Y_train)):
    # fold k 分割
    Xk_train, Xk_validate = X_train[idx_train], X_train[idx_validate]
    Yk_train, Yk_validate = Y_train[idx_train], Y_train[idx_validate]
    # 正規化
    scaler_std = StandardScaler()
    # 注意 fit() 只套用在這個 fold 的訓練分割
    scaler_std.fit(Xk_train)
    Xk_train = scaler_std.transform(Xk_train)
    Xk_validate = scaler_std.transform(Xk_validate)
    # 訓練模型
    sklr2.fit(Xk_train, Yk_train)
    # 驗證這個模型的正確率
    accuracy = sklr2.score(Xk_validate, Yk_validate)
    # 記錄這個模型的驗證正確率
    scores2.append(accuracy)
    # 結束這個 fold，顯示訊息
    print('Fold #{} accuracy = {:.3f}'.format(k, accuracy))
# 交叉驗證完成，顯示平均分數
print('\n10-fold cross-validation result: average accuracy = {:.3f}'.format(np.mean(scores2)))

# 再把剛剛不好處理的 cross-validate 用 Pipeline 做一次
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedShuffleSplit, cross_validate

# 定義包含正規化前處理的 pipeline
pipelr3 = Pipeline(
    steps=[
        ('zscaler', StandardScaler()),
        ('lr', LogisticRegression(max_iter=100))
    ]
)

# 10-fold 交叉驗證，使用所有處理核心平行運算
scores3 = cross_validate(pipelr3, X_train, Y_train, cv=StratifiedShuffleSplit(n_splits=10), n_jobs=-1)

# 顯示每個 fold 的 accuracy 分數
dfscores3 = pd.DataFrame(scores3)
dfscores3

# 觀察平均分數
print('[Average scores]:\n', dfscores3.mean(axis=0))

Image(filename='images/ConfusionMatrix.png', width=480)

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedShuffleSplit, cross_validate
# 示範各 metrics 的使用
from sklearn.metrics import (
    make_scorer,
    accuracy_score,
    recall_score,
    roc_auc_score,
    roc_curve,
    plot_roc_curve,
    confusion_matrix
)

# 與前面範例相同的 pipeline 定義
pipelr4 = Pipeline(
    steps=[
        ('zscaler', StandardScaler()),
        ('lr', LogisticRegression(max_iter=100))
    ]
)

# 交叉驗證過程可以指定多項評估指標
combo_scores = {
    'roc_auc': make_scorer(roc_auc_score),
    'accuracy': make_scorer(accuracy_score),
    'sensitivity': make_scorer(recall_score, pos_label=1),
    'specificity': make_scorer(recall_score, pos_label=0)
}

# 10-fold 交叉驗證，使用所有處理核心平行運算
scores4 = cross_validate(
    pipelr4,
    X_train,
    Y_train,
    scoring=combo_scores,
    cv=StratifiedShuffleSplit(n_splits=10),
    n_jobs=-1
)

# 顯示每個 fold 的指標分數，明確指定需要的欄位
col_names = ['test_roc_auc', 'test_accuracy', 'test_sensitivity', 'test_specificity']
dfscores4 = pd.DataFrame({k:v for k,v in scores4.items() if k in col_names})
dfscores4

# 觀察平均分數
print('[Average scores]:\n', dfscores4.mean(axis=0))

# 與前面範例相同的 pipeline 定義，假設超參數就是最佳
pipelr5 = Pipeline(
    steps=[
        ('zscaler', StandardScaler()),
        ('lr', LogisticRegression(max_iter=100))
    ]
)

# 直接訓練整個訓練集
pipelr5.fit(X_train, Y_train)
# 已訓練模型用來預測測試集
Y_predict = pipelr5.predict(X_test)

# confusion matrix
TN, FP, FN, TP = confusion_matrix(Y_test, Y_predict).ravel()
print('[Confusion Matrix]:\n  | TP: {} | FP: {} |\n  | FN: {} | TN: {} |\n'.format(TP, FP, FN, TN))

# sensitivity and specificity
sensitivity = recall_score(Y_test, Y_predict, pos_label=1)
specificity = recall_score(Y_test, Y_predict, pos_label=0)
print('[Testing Scores]:\n  * Sensitivity: {:.3f}\n  * Specificity: {:.3f}'.format(sensitivity, specificity))

# accuracy and ROC AUC
accuracy = accuracy_score(Y_test, Y_predict)
rocauc = roc_auc_score(Y_test, pipelr5.predict_proba(X_test)[:,1])
print('  * Accuracy: {:.3f}\n  * ROC AUC: {:.3f}'.format(accuracy, rocauc))

# 畫出 ROC 曲線
_, ax5 = plt.subplots(figsize=(8, 7))
plot_roc_curve(pipelr5, X_test, Y_test, name='5th Logistic Regression', ax=ax5)
plt.show()

# 維基百科範例： Logistic Regression 模型由念書的時間預測考試會pass的機率
Image(url='https://upload.wikimedia.org/wikipedia/commons/6/6d/Exam_pass_logistic_curve.jpeg', embed=True, width=600)

# 取得模型預測測試集的機率 
predict_proba = pipelr5.predict_proba(X_test)
# 機率 p(y=0), p(y=1)
py0, py1 = predict_proba[:,0], predict_proba[:,1]
pd.DataFrame({'y label': Y_test, 'p(y=0)': py0, 'p(y=1)': py1})

# 二元分類模型中，擬合的機率固定就是 p(y=0) + p(y=1) = 1
py0 + py1

# 取 p(y=1) 由大到小遞減排序，回傳排序後的原序號
py1_dec_index = np.argsort(py1, kind='stable')[::-1]

# 取得遞減的 p(y=1) 
py1_dec = py1[py1_dec_index]
print('p(y=1) decreasing sorted:\n', py1_dec)

# y=1 的標籤轉成 True/False 的陣列
y_true = (Y_test == 1)

# 標籤對應 p(y=1) 的序號
y_true = y_true[py1_dec_index]
print('corresponding label y=1:\n', y_true)

# 每個 p(y=1) 都是 threshold，取累計 True Positive
tps = np.cumsum(y_true, dtype=np.float64)
print('Cumulative True Positives:\n', tps)

# 取累計 False Positive
fps = np.cumsum((1 - y_true), dtype=np.float64)
print('Cumulative False Positives:\n', fps)

tpr = tps / tps[-1]
print('True Positive Rates:\n', tpr)

fpr = fps / fps[-1]
print('False Positive Rates:\n', fpr)

# 畫 ROC 曲線前補上 (0,0) 的點
tpr = np.r_[0.0, tpr]
print('True Positive Rates:\n', tpr)

fpr = np.r_[0.0, fpr]
print('False Positive Rates:\n', fpr)

_, ax6 = plt.subplots(figsize=(8, 7))
# plot ROC curve
ax6.plot(fpr, tpr, label='5th Logistic Regression (AUC = {:.3f})'.format(rocauc), lw=2)

# the diagonal line
ax6.plot([0, 1], [0, 1], linestyle='--', lw=2, color='#808080', alpha=0.8)

ax6.set(xlim=[-0.05, 1.05], ylim=[-0.05, 1.05])
ax6.set_xlabel('False Positive Rate', fontsize=16)
ax6.set_ylabel('True Positive Rate', fontsize=16)
ax6.legend(loc="lower right", fontsize=14)
plt.show()