概率校准

用户3577892

发布于 2020-06-10 17:48:54

2.6K0

发布于 2020-06-10 17:48:54

文章被收录于专栏：数据科学CLUB

使用sklearn自动生成二分类数据集，划分训练集、验证集和测试集对不同的分类器，画出可靠性曲线在训练集上：在验证集上如何进行概率校准（probability calibration）方法1：Platt Scaling方法2：Isotonic Regression 保序回归scikit-learn 提供了执行概率预测校准的两种方法的API评价：Brier score

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.linear_model import LogisticRegression
from sklearn.isotonic import IsotonicRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import LinearSVC
from sklearn.calibration import calibration_curve, CalibratedClassifierCV

对于一个二分类问题，我们希望：分类器预测某个样本属于正类的概率是0.8，那么就应当说明有80%的把握认为该样本属于正类，或者100个概率为0.8的里面有80个确实属于正类。。这个目的也是出于实际业务的考虑。（例如，在信贷风控中，将预测的客户违约概率与真实违约概率对标，即模型风险概率能够代表真实的风险等级。）

由于我们无法获知真实的条件概率，通常用观测样本的标签来统计代替，并用可靠性曲线（Reliability Curve）来直观展示当前模型的输出结果与真实结果的偏差。画出可靠性曲线的方法步骤是：

选定一个阈值，将[0,1]划分成若干小区间
计算在这个区间上正样本率
描点连线。曲线与越接近，那么说明模型被校准得很好。（分类器输出的概率能够代表真实的概率）

下面使用使用sklearn自动生成的二分类数据集画出几种基本的二分类模型的可靠性曲线。

使用sklearn自动生成二分类数据集，划分训练集、验证集和测试集

X, y = datasets.make_classification(n_samples=9600, n_features=20,
                                    n_informative=2, n_redundant=2, random_state=2020)
train_samples = 5600  # Samples used for training the models

X_train = X[:train_samples]
X_test = X[train_samples:]
y_train = y[:train_samples]
y_test = y[train_samples:]

X_val, X_test, y_val, y_test = \
    train_test_split(X_test, y_test, test_size=0.4, random_state=2020)

X_train.shape, y_train.shape, X_test.shape, y_test.shape, \
    X_val.shape, y_val.shape

((5600, 20), (5600,), (1600, 20), (1600,), (2400, 20), (2400,))

对不同的分类器，画出可靠性曲线

# Create classifiers
lr = LogisticRegression(solver='liblinear')
gnb = GaussianNB()
svc = LinearSVC(C=1.0)
rfc = RandomForestClassifier()

# Plot calibration plots
plt.figure(figsize=(8, 8))
ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
ax2 = plt.subplot2grid((3, 1), (2, 0))

ax1.plot([0, 1], [0, 1], "k:", label="Perfectly calibrated")
for clf, name in [(lr, 'Logistic'),
                  (gnb, 'Naive Bayes'),
                  (svc, 'Linear SVC'),
                  (rfc, 'Random Forest')]:
    clf.fit(X_train, y_train)
    if hasattr(clf, "predict_proba"):
        prob_pos = clf.predict_proba(X_test)[:, 1]
    else:  # use decision function   svc 没有 predict_proba 方法
        prob_pos = clf.decision_function(X_test)
        prob_pos = \
            (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())
    fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, prob_pos, n_bins=10)

    ax1.plot(mean_predicted_value, fraction_of_positives, "s-",
             label="%s" % (name, ))

    ax2.hist(prob_pos, range=(0, 1), bins=10, label=name,
             histtype="step", lw=2)

ax1.set_ylabel("Fraction of positives")
ax1.set_ylim([-0.05, 1.05])
ax1.legend(loc="lower right")
ax1.set_title('Calibration plots  (reliability curve)')

ax2.set_xlabel("Mean predicted value")
ax2.set_ylabel("Count")
ax2.legend(loc="upper center", ncol=2)

plt.tight_layout()
plt.show()

曲线分析：

Logistic回归预测概率比较准（模型本身的特点，经验如此）。
朴素贝叶斯过于自信（可能由于冗余特征所致，违背了特征独立性前提）呈反sigmoid曲线。
SVM很不自信呈sigmoid曲线，随机森林也是。

一般来说，boosting 的树模型和朴素贝叶斯呈现比较差的可靠性曲线，而神经网络（逻辑回归）和 bagging 的树模型相对较好。

分类准确率：

在训练集上：

for clf, name in [(lr, 'Logistic'),
                  (gnb, 'Naive Bayes'),
                  (svc, 'Linear SVC'),
                  (rfc, 'Random Forest')]:
    print(name, 'Train Set: ', clf.score(X_train, y_train))

Logistic Train Set:  0.8455357142857143
Naive Bayes Train Set:  0.8448214285714286
Linear SVC Train Set:  0.8457142857142858
Random Forest Train Set:  0.9932142857142857

在验证集上

for clf, name in [(lr, 'Logistic'),
                  (gnb, 'Naive Bayes'),
                  (svc, 'Linear SVC'),
                  (rfc, 'Random Forest')]:
    print(name, 'Test Set: ', clf.score(X_test, y_test))

Logistic Test Set:  0.821875
Naive Bayes Test Set:  0.821875
Linear SVC Test Set:  0.820625
Random Forest Test Set:  0.885

可以看到，即便这些分类器的准确率相似，但是可靠性曲线却大相径庭

那么，如何对分类器进行概率校准呢，使得模型输出的概率能够近似代表实际样本为正的概率？

如何进行概率校准（probability calibration）

方法1：Platt Scaling

适用于呈现sigmoid型可靠性曲线的分类器。

将模型输出放入逻辑回归中训练，最后将逻辑回归的结果作为模型的校准结果。如果原来的二分类器得到的结果是，那么就把当做新的训练数据用Logistic回归训练（因为经验告诉我们LR的表现总是好的）。

为了避免过拟合，这两次训练要使用不同的数据集！

在提出Platt Scaling的论文中，作者还建议将Logistic回归的标签从0，1变成和，其中和分别是正样本和负样本的数量。（分子分母的常数以拉普拉斯平滑得出）

下面对Linear SVC模型应用改方法进行概率校准

# 这里用验证集上的数据训练 Logistic回归
prob_pos = svc.decision_function(X_val)
prob_pos = (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())

# 训练 lr
svc_lr = LogisticRegression(solver='liblinear')
svc_lr.fit(prob_pos.reshape(-1, 1), y_val)

LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
                   intercept_scaling=1, l1_ratio=None, max_iter=100,
                   multi_class='warn', n_jobs=None, penalty='l2',
                   random_state=None, solver='liblinear', tol=0.0001, verbose=0,
                   warm_start=False)

prob_pos = svc.decision_function(X_test)
prob_pos = (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())
prob_pos = svc_lr.predict_proba(prob_pos.reshape(-1, 1))[:, 1]

plt.figure(figsize=(12, 4.5))
plt.subplot(121)
# 先算出在测试集上的预测概率
prob_pos = svc.decision_function(X_test)
prob_pos = (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())
fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, prob_pos, n_bins=10)
plt.plot([0, 1], [0, 1], 'k:', mean_predicted_value, fraction_of_positives, 's-')
# 用测试集上的预测概率得到校准模型的输出
prob_pos = svc_lr.predict_proba(prob_pos.reshape(-1, 1))[:, 1]
plt.title('before calibration')
plt.subplot(122)
fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, gnb.predict_proba(X_test)[:, 1], n_bins=10)
plt.plot([0, 1], [0, 1], 'k:', mean_predicted_value, fraction_of_positives, 's-')
plt.title('after calibration')
plt.show()

左边是原可靠性曲线，右边是经过校准后的可靠性曲线

方法2：Isotonic Regression 保序回归

保序回归本身是回归分析的一种，不过不太常见，下面对保序回归做一个简单介绍给定学习样本和标签，其中

希望：取最小值

保序回归的结果是分段函数。

看一个例子：假设药物使用量为数组，病人对药物的反应量为，而由于个体的原因，不是一个单调函数(即：存在波动)，如果我们按照药物反应排序，对应的X就会成为乱序，失去了研究的意义。而我们的研究的目的是观察随着药物使用量的递增，病人的平均反应状况。在这种情况下，使用保序回归，既不改变X的排列顺序，又求的Y的平均值状况。

x = np.arange(0, 100, 2)
y = np.random.randint(-10, 10, size=x.size) + 8 * np.log(x+2)

ir = IsotonicRegression()
y_ = ir.fit_transform(x, y)

plt.plot(x, y, '.', x, y_, 'r.')
plt.show()

生成一列总体来说递增的数据，如上图红线的结果就是对原样本进行保序回归的拟合结果，是一个分段函数

那么，如何应用保序回归来进行概率校准呢？假设有数据集和预测结果。 对预测结果按 排序再进行保序回归，就是校准好的模型

保序回归其实是做出了[0,1]上的一条分段单调分减的曲线，然后输出 (x=原始模型的预测值) 的y值。和 Platt Scaling 类似，为了保证不引入偏差，用作校准的数据集应该和训练模型的数据集不同。

下面对之前的朴素贝叶斯模型应用保序回归方法进行概率校准：训练集用来训练朴素贝叶斯模型，验证集用来训练保序回归，画出测试集上的可靠性曲线

data = sorted(zip(gnb.predict_proba(X_val)[:, 1], y_val), key=lambda x: x[0])

ix = np.array([i[0] for i in data])
iy = np.array([i[1] for i in data])

ir = IsotonicRegression()
iy_ = ir.fit_transform(ix, iy)

校准后的概率：

prob = ir.transform(gnb.predict_proba(X_test)[:, 1])

plt.figure(figsize=(12, 4.5))
plt.subplot(121);plt.plot([0, 1], [0, 1], 'k:')
fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, prob, n_bins=10)
plt.plot(mean_predicted_value, fraction_of_positives, 's-')
plt.title('after calibration')
plt.subplot(122);plt.plot([0, 1], [0, 1], 'k:')
fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, gnb.predict_proba(X_test)[:, 1], n_bins=10)
plt.title('before calibration')
plt.plot(mean_predicted_value, fraction_of_positives, 's-');plt.show()

左图是用保序回归校准后的可靠性曲线，右边是原始的可靠性曲线

经验来说，Isotonic Regression方法有更广阔的适应性，而 Platt Scaling 方法比较适用于呈 sigmoid 型可靠性曲线的分类器。不幸的是，当数据量小时，Isotonic Regression方法更容易过拟合。

以上介绍了概率校准的两种方式并且用代码实践了。

scikit-learn 提供了执行概率预测校准的两种方法的API

sklearn.calibration 是用来进行概率校准的模块

内含两个函数：

sklearn.calibration.CalibratedClassifierCV: Probability calibration with isotonic regression or sigmoid.
sklearn.calibration.calibration_curve: Compute true and predicted probabilities for a calibration curve.

calibration_curve 函数用来画出可靠性曲线（上面已经用过）

CalibratedClassifierCV是一个类，接收参数：

base_estimator: 基础模型
method: 'sigmoid' 或者 'isotonic'，默认'sigmoid'
cv: 交叉验证数

可以使用CalibratedClassifierCV来进行概率校准，下面以method='isotonic'为例

lr = CalibratedClassifierCV(lr, method='isotonic')
gnb = CalibratedClassifierCV(gnb, method='isotonic')
svc = CalibratedClassifierCV(svc, method='isotonic')
rfc = CalibratedClassifierCV(rfc, method='isotonic')

# Plot calibration plots
plt.figure(figsize=(8, 8))
ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
ax2 = plt.subplot2grid((3, 1), (2, 0))

ax1.plot([0, 1], [0, 1], "k:", label="Perfectly calibrated")
for clf, name in [(lr, 'Logistic'),
                  (gnb, 'Naive Bayes'),
                  (svc, 'Linear SVC'),
                  (rfc, 'Random Forest')]:
    clf.fit(X_train, y_train)
    if hasattr(clf, "predict_proba"):
        prob_pos = clf.predict_proba(X_test)[:, 1]
    else:  # use decision function
        prob_pos = clf.decision_function(X_test)
        prob_pos = \
            (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())
    fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, prob_pos, n_bins=10)

    ax1.plot(mean_predicted_value, fraction_of_positives, "s-",
             label="%s" % (name, ))

    ax2.hist(prob_pos, range=(0, 1), bins=10, label=name,
             histtype="step", lw=2)

ax1.set_ylabel("Fraction of positives")
ax1.set_ylim([-0.05, 1.05])
ax1.legend(loc="lower right")
ax1.set_title('reliability curve (using isotonic calibration)')

ax2.set_xlabel("Mean predicted value")
ax2.set_ylabel("Count")
ax2.legend(loc="upper center", ncol=2)

plt.tight_layout()
plt.show()