Python机器学习(原书第3版)- 第3章scikit-learn 机器学习分类器
3.1 选择分类算法 监督机器学习算法训练的主要步骤:1)选择特征并收集训练样本;2)选择度量性能的指标;3)选择分类器并优化算法;4)评估模型的性能;5)调整算法。3.2 了解scikit-learn第一步——训练感知器 我们仍然使用鸢尾花数据集来训练感知器,首先使用数据集的两个特征来实现可视化:将150个鸢尾花样本的花瓣长度和宽度存入特征矩阵X,把相应的品种分类标签存入向量y,实现
3.1 选择分类算法
监督机器学习算法训练的主要步骤:
1)选择特征并收集训练样本;
2)选择度量性能的指标;
3)选择分类器并优化算法;
4)评估模型的性能;
5)调整算法。
3.2 了解scikit-learn第一步——训练感知器
我们仍然使用鸢尾花数据集来训练感知器,首先使用数据集的两个特征来实现可视化:将150个鸢尾花样本的花瓣长度和宽度存入特征矩阵X,把相应的品种分类标签存入向量y,实现代码如下:
from sklearn import datasets
import numpy as np
iris = datasets.load_iris()
X = iris.data[:, [2, 3]]
y = iris.target
print('Class labels:',np.unique(y))
可以看到输出结果为:
Class labels: [0 1 2]
代码中的np.unique(y)函数把返回的三个独立分类标签存储在 iris.target中,以整数(0、1、2)分别存储鸢尾花Iris-setosa、Iris-versicolor和Iris-virginica的标签。
虽然许多scikit-learn函数和分类方法也能处理字符串形式的分类标签,但推荐采用整数标签。其优点在于可以避免技术故障、所占内存更小同时提高模型的计算性能。
这之后我们进一步将数据集分割成单独的训练数据集和测试数据集。实现代码为:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1, stratify=y)
在这段代码中,我们把X和y数组随机拆分为30%的测试数据集(45个样本)和70%的训练数据集(105个样本)。
最后,再定义stratify=y 获得内置的分层支持。这种分层可以返回与输入数据集的分类标签相同比例的训练数据集和测试训练集。可以调用Numpy和bincount函数统计数组中每个值出现数,以验证数据,实现代码为:
print('Labels counts in y:', np.bincount(y))
print('Lbels counts in y_train:', np.bincount(y_train))
print('Lbels counts in y_test:', np.bincount(y_test))
执行后输出结果为:
Labels counts in y: [50 50 50]
Lbels counts in y_train: [35 35 35]
Lbels counts in y_test: [15 15 15]
之前我们学到许多机器学习和优化算法需要进行特征缩放来获得最佳性能,在这里我们使用scikit-learn的preprocessing模块中的StandardScaler类来标准化特征,实现代码如下:
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)
在上述代码中,我们加载了preprocessing模块中的StandardScaler类初始化一个新的StandardScaler对象,然后将其分配给变量sc,再调用StandardScaler的fit方法对训练数据的每个特征维度估计参数μ\muμ(样本均值)和σ\sigmaσ(标准差)进行估算。最后调用transform方法,利用估计的参数μ\muμ和σ\sigmaσ标准化训练数据。
在完成以上步骤后,我们就可以开始训练感知器模型。通过调用一对其余(OvR) 方法,把三类花的数据同时提交给感知器。具体实现代码如下:
from sklearn.linear_model import Perceptron
ppn = Perceptron(eta0=0.1, random_state=1)
ppn.fit(X_train_std, y_train)
这便是从linear_model模块中加载Perceptron类之后,初始化新的Perceptron对象,然后调用fit方法对模型进行训练。依然是再每次迭代后对训练数据集用random_state参数进行洗牌,以确保初始结果可以重现。
训练完模型之后,就可以调用predict方法做预测,具体代码如下:
y_pred = ppn.predict(X_test_std)
print('Misclassified examples: %d' % (y_test != y_pred).sum())
输出结果得到:
Misclassified examples: 1
意为在感知器处理数据时出现过1次错误分类,因此测试数据集上的分类错误率大概为2.2%(1/45≈0.022)。分类准确率即为$1-2.2%=97.8% $。在scikit-learn中也可以直接实现性能指标,例如我们可以通过调用metrics模块计算测试数据集上感知器的分类准确率,代码为:
from sklearn.metrics import accuracy_score
print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))
可以直接得到结果:
Accuracy: 0.978
或者,也可以直接使用scikit-learn分类其中的评分方法,综合调用predict和accuracy_score计算出分类器的预测准确率:
print('Accuracy: %.3f' % ppn.score(X_test_std, y_test))
同样可以得到上述结果。
最后,我们再利用第2章中的plot_decision_regions函数绘制新训练的感知器模型的决策区,并以可视化的方式展示出不同花朵样本的效果。但要略加修改通过圆圈来突出显示来自测试数据集的样本,代码如下:
#绘制决策区
from matplotlib.colors import ListedColormap
import matplotlib.pyplot as plt
def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):
#setup marker generator and color map
markers = ('s', 'x', 'o', '^', 'v')
colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
cmap = ListedColormap(colors[:len(np.unique(y))])
#plot the decision surface
x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
np.arange(x2_min, x2_max, resolution))
Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
Z = Z.reshape(xx1.shape)
plt.contourf(xx1, xx2, Z, alpha=0.3, cmap=cmap)
plt.xlim(xx1.min(), xx1.max())
plt.ylim(xx2.min(), xx2.max())
for idx, cl in enumerate(np.unique(y)):
plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=colors[idx], marker=markers[idx], label=cl, edgecolor='black')
#highlight test examples
if test_idx:
#plot all examples
X_test, y_test = X[test_idx, :], y[test_idx]
plt.scatter(X_test[:, 0], X_test[:, 1], facecolors='none', edgecolor='black', alpha=1, linewidth=1, marker='o', s=100, label='test set')
#修改plot_decision_regions函数
X_combined_std = np.vstack((X_train_std, X_test_std))
y_combined = np.hstack((y_train, y_test))
plot_decision_regions(X=X_combined_std, y=y_combined, classifier=ppn, test_idx=range(105, 150))
plt.xlabel('petal length[standardized]')
plt.ylabel('petal width[standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
得到结果如下图所示:
从上图可以看到,此时三类花并不能被线性决策边界完全分离。第二章曾讨论过感知器算法从不收敛于不完全线性可分离的数据集,这就是为什么在实践中不推荐使用感知器算法。接下来我们研究更强大的可以收敛到代价最小值的线性分类器,即使这些类不是完全线性可分的。
3.3 基于逻辑回归的分类概率建模
虽然感知器规则提供了良好且易用的入门级机器学习分类算法,但其最大的缺点是:如果类不是完全线性可分的,那么它将永远不收敛。直观的说,原因是权重在不断更新,因为每次迭代至少会有一个错误分类样本存在,虽然可以通过改变学习速率、增加迭代次数来缓解这个问题,但仍然可能会出现永远不会收敛的问题。
为了提高效率,我们学习逻辑回归。这是一种解决线性二元分类问题的算法。(注意:逻辑回归是一种分类模型,而不是回归模型。)
3.3.1 逻辑回归和条件概率
逻辑回归是一种很容易实现的分类模型,但仅在线性可分类上表现不错。是一种应用广泛的 二分类模型,而且可以利用 OVR 技术扩展到多元分类。
在了解逻辑回归远离之前,首先要了解让步比(oods ratio)。
让步比的定义为:某一特定时间发生的概率。数学公式定义为 p1−p\displaystyle\frac{p}{1-p}1−pp ,其中 ppp 代表正事件的概率(正事件并不一定意味着好的事件,指的是要预测的事件。例如:病人有某种疾病的可能性)。可以认为正事件的分类标签为 y=1y=1y=1。 可以进一步定义logit函数,这仅仅是让步比的对数形式(log-odds):
logit(p)=logp(1−p)logit(p)=log\frac{p}{(1-p)}logit(p)=log(1−p)p
这里的log是自然对数,与计算机科学中的通用惯例一致。logit函数的输入值取值在0到1之间,并将其转换为整个实数范围的值,可以用它来表示特征值和对数概率(log-odds)之间的线性关系,写为:
logit(p(y=1∣x))=ω0x0+ω1x1+…+ωmxm=∑i=0mωixi=ωTxlogit(p(y=1|x))=\omega_0x_0+\omega_1x_1+…+\omega_mx_m=\displaystyle\sum_{i=0}^m\omega_ix_i=\omega^Txlogit(p(y=1∣x))=ω0x0+ω1x1+…+ωmxm=i=0∑mωixi=ωTx
其中p(y=1∣x)p(y=1|x)p(y=1∣x)是给定特征xxx ,某个特定样本属于类1的条件概率。
实际上,我们感兴趣的是预测某个样本属于某个特定类的概率,它是logit函数的你函数,这也被称为逻辑Sigmoid函数,属于Sigmoid函数,被称为Logistic函数(大多时间大家直接称为Sigmoid函数,为和同属于Sigmoid函数的Tanh函数区分,这里只称为Logistic函数。)其公式表示为:
ϕ(z)=11+e−z\phi(z)=\displaystyle\frac{1}{1+e^{-z}}ϕ(z)=1+e−z1
其中z为净输入,是权重和样本特征的线性组合。具体表示为:
z=ωTx=ω0x0+ω1x1+…+ωmxmz=\omega^Tx=\omega_0x_0+\omega_1x_1+…+\omega_mx_mz=ωTx=ω0x0+ω1x1+…+ωmxm
使用-7到7之间的一些值简单地绘制出Logistic函数来观察具体的情况,代码如下:
#Logistic
import numpy as np
import matplotlib.pyplot as plt
# 定义Logistic激活函数
def logistic(z):
return 1.0 / (1.0 + np.exp(-z))
# 生成输入值范围
z = np.arange(-7, 7, 0.1)
phi_z = logistic(z)
plt.plot(z, phi_z)
plt.axvline(0, color='k')
plt.ylim(-0.1, 1.1)
plt.xlabel('z')
plt.ylabel('$\phi(z)$')
plt.yticks([0, 0.5, 1])
ax = plt.gca()
ax.yaxis.grid(True)
plt.tight_layout()
plt.show()
看到S型曲线如下图所示:
从图中可以看到,当z趋向无限大时,ϕ(z)\phi(z)ϕ(z)的值接近于1,因为当z值很大时,e−ze^{-z}e−z的值会变得很小。因此可以得到这样的结论:Logistic函数将以实数值作为输入,并在截距为 ϕ(z)=0.5\phi(z)=0.5ϕ(z)=0.5时转换为[0,1]范围内的值。
为了直观地理解逻辑回归模型,我们将其与第二章中介绍的Adaline联系起来,在Adaline中用恒等函数ϕ(z)=z\phi(z)=zϕ(z)=z作为激活函数,在逻辑回归中,只是简单地将Logistic函数作为激活函数。二者的区别如下图所示:
Logistic函数的输出被解释为特定样本属于类1的概率,ϕ(z)=P(y=1∣x;ω)\phi(z)=P(y=1|x; \omega)ϕ(z)=P(y=1∣x;ω),这里假设特征x被权重ω\omegaω参数化。例如,对于某种花的样本,计算出ϕ(z)=0.8\phi(z)=0.8ϕ(z)=0.8,说明该样本属于Iris-versicolor的概率是80%。因此该样本属于Iris-setosa的概率可以计算为P(y=0∣x;ω)=1−P(y=1+x;ω)=0.2P(y=0|x; \omega)=1-P(y=1+x; \omega)=0.2P(y=0∣x;ω)=1−P(y=1+x;ω)=0.2。此时,预测概率从可以通过阈值函数简单地转换为二元输出:
y^={1,ϕ(z)≥0.50,else\hat{y}=\begin{cases} 1, \phi(z)\geq 0.5 \\ 0, else \\ \end{cases}y^={1,ϕ(z)≥0.50,else
那么Logistic函数图就可以等同于:
y^={1,z≥00,else\hat{y}=\begin{cases} 1, z\geq 0 \\ 0, else \\ \end{cases}y^={1,z≥00,else
3.3.2 学习逻辑代价函数的权重
在第二章中,我们将代价函数定义为误差平方和(SSE),计算公式为:
J(ω)=∑i12(ϕ(z(i)−y(i))2J(\omega)=\sum_i\frac{1}{2}(\phi(z^{(i)}-y^{(i)})^2J(ω)=i∑21(ϕ(z(i)−y(i))2
为了在Adaline分类模型中学习权重ω\omegaω,我们简化了函数。那么要如何才能得到逻辑回归的代价函数,首先我们需要定义最大似然函数L。假设数据集中的每个样本都是相互独立的,有最大似然函数L公式为:
L(ω)=P(y∣x;ω)=∏i=1nP(y(i)∣x(i);ω)=∏i=1n(ϕ(z(i)))y(i)(1−ϕ(z(i)))1−y(i)L(\omega)=P(y|x; \omega)=\prod_{i=1}^nP(y^{(i)}|x^{(i)}; \omega)=\prod_{i=1}^n(\phi(z^{(i)}))^{y^{(i)}}(1-\phi(z^{(i)}))^{1-y^{(i)}}L(ω)=P(y∣x;ω)=i=1∏nP(y(i)∣x(i);ω)=i=1∏n(ϕ(z(i)))y(i)(1−ϕ(z(i)))1−y(i)
因为在实践中很容易最大化该方程的自然对数(求其最大值),所以我们再次定义对数似然函数:
l(ω)=logL(ω)=∑i=1n[y(i)log(ϕ(z))+(1−y(i))log(1−ϕ(z(i)))]l(\omega)=logL(\omega)=\sum_{i=1}^n[y^{(i)}log(\phi(z))+(1-y^{(i)})log(1-\phi(z^{(i)}))]l(ω)=logL(ω)=i=1∑n[y(i)log(ϕ(z))+(1−y(i))log(1−ϕ(z(i)))]
改写对数似然函数作为代价函数JJJ,用梯度下降方法最小化代价函数有:
J(ω)=∑i=1n[−y(i)log(ϕ(z(i)))−(1−y(i))log(1−ϕ(z(i)))]J(\omega)=\sum_{i=1}^n[-y^{(i)}log(\phi(z^{(i)}))-(1-y^{(i)})log(1-\phi(z^{(i)}))]J(ω)=i=1∑n[−y(i)log(ϕ(z(i)))−(1−y(i))log(1−ϕ(z(i)))]
为了更好地理解这个代价函数,我们计算一个训练样本的代价来尝试理解:
J(ϕ(z),y;ω)=−y(i)log(ϕ(z))−(1−y)log(1−ϕ(z))J(\phi(z), y; \omega)=-y^{(i)}log(\phi(z))-(1-y)log(1-\phi(z))J(ϕ(z),y;ω)=−y(i)log(ϕ(z))−(1−y)log(1−ϕ(z))
从上面的方程可以看到,如果y=0y=0y=0,第一项为0,如果y=1y=1y=1,第二项为0:
J(ϕ(z),y;ω)={−log(ϕ(z)),y=1−log(1−ϕ(z)),y=0J(\phi(z), y; \omega)=\begin{cases} -log(\phi(z)), y=1 \\ -log(1-\phi(z)), y=0 \\ \end{cases}J(ϕ(z),y;ω)={−log(ϕ(z)),y=1−log(1−ϕ(z)),y=0
我们也可以绘制一张图,来说明对于不同的ϕ(z)\phi(z)ϕ(z)值,分类一个训练样本的代价,实现代码如下:
#分类训练样本的代价
def cost_1(z):
return - np.log(logistic(z))
def cost_0(z):
return - np.log(1 - logistic(z))
z = np.arange(-10, 10, 0.1)
phi_z = logistic(z)
c1 = [cost_1(x) for x in z]
plt.plot(phi_z, c1, label='J(w) if y=1')
c0 = [cost_0(x) for x in z]
plt.plot(phi_z, c0, linestyle='--', label='J(w) if y=0')
plt.ylim(0, 5.1)
plt.xlim([0,1])
plt.xlabel('$\phi(z)$')
plt.ylabel('J(w)')
plt.legend(loc='best')
plt.tight_layout()
plt.show()
运行代码得到下图:
上图绘制了0到1范围内,y轴是x轴的logistic激活函数相关联的逻辑代价。可以看到,如果正确地预测样本属于类1,代价就回接近0;可以在y轴上看到,如果正确地预测y=0,那么代价也接近于0.然而, 如果预测错误,代价就会趋于无穷大。关键就在于用越来越大的代价惩罚错误的预测。
3.3.3 将Adaline实现转换为一个逻辑回归算法
如果想自己动手实现逻辑回归,可以直接用新的代价函数取代第二章Adaline实现中的代价函数J,例如:
J(ω)=−∑iy(i)log(ϕ(z(i)))+(1−y(i))log(1−ϕ(z(i)))J(\omega)=-\sum_iy^{(i)}log(\phi(z^{(i)}))+(1-y^{(i)})log(1-\phi(z^{(i)}))J(ω)=−i∑y(i)log(ϕ(z(i)))+(1−y(i))log(1−ϕ(z(i)))
在对训练样本进行分类的过程中,用该公式来计算每次迭代的额代价。另外还需要用Logistic激活函数来替换线性激活函数,同时改变阈值函数,返回类标签0和1,而不是返回-1和1.由此实现逻辑回归的代码,如下:
class LogisticRegressionGD(object):
"""Logistic Regression Classifier using gradient descent.
Parameters
------------
eta : float
Learning rate (between 0.0 and 1.0)
n_iter : int
Passes over the training dataset.
random_state : int
Random number generator seed for random weight
initialization.
Attributes
-----------
w_ : 1d-array
Weights after fitting.
cost_ : list
Logistic cost function value in each epoch.
"""
def __init__(self, eta=0.05, n_iter=100, random_state=1):
self.eta = eta
self.n_iter = n_iter
self.random_state = random_state
def fit(self, X, y):
""" Fit training data.
Parameters
----------
X : {array-like}, shape = [n_examples, n_features]
Training vectors, where n_examples is the number of examples and
n_features is the number of features.
y : array-like, shape = [n_examples]
Target values.
Returns
-------
self : object
"""
rgen = np.random.RandomState(self.random_state)
self.w_ = rgen.normal(loc=0.0, scale=0.01, size=1 + X.shape[1])
self.cost_ = []
for i in range(self.n_iter):
net_input = self.net_input(X)
output = self.activation(net_input)
errors = (y - output)
self.w_[1:] += self.eta * X.T.dot(errors)
self.w_[0] += self.eta * errors.sum()
# note that we compute the logistic `cost` now
# instead of the sum of squared errors cost
cost = -y.dot(np.log(output)) - ((1 - y).dot(np.log(1 - output)))
self.cost_.append(cost)
return self
def net_input(self, X):
"""Calculate net input"""
return np.dot(X, self.w_[1:]) + self.w_[0]
def activation(self, z):
"""Compute logistic sigmoid activation"""
return 1. / (1. + np.exp(-np.clip(z, -250, 250)))
def predict(self, X):
"""Return class label after unit step"""
return np.where(self.net_input(X) >= 0.0, 1, 0)
当拟合逻辑回归模型时,时刻谨记该模型只适用于二元分类,所以这里我们只考虑Iris-setosa和Iris-versicolor量中华,并验证逻辑回归的有效性,代码如下:
X_train_01_subset = X_train[(y_train == 0) | (y_train == 1)]
y_train_01_subset = y_train[(y_train == 1) | (y_train == 1)]
lrgd = LogisticRegressionGD(eta=0.05, n_iter=1000, random_state=1)
lrgd.fit(X_train_01_subset, y_train_01_subset)
plot_decision_regions(X=X_train_01_subset, y=y_train_01_subset, classifier=lrgd)
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
绘制出图像为:
3.3.4 用scikit-learn训练逻辑回归模型
在以上逻辑回归模型中,只支持二元分类。接下来我们尝试如何使用scikit-learn或者更优化的逻辑回归来实现多元分类的场景。在scikit-learn的新版中,将会自动选择多元分类、多项式或者OVR技术。我们尝试使用scikit-learn中的方法来在三种花的标准化训练数据集上训练模型。
注意,在这里,为了方便解释,我们将设置multi_class=‘ovr’。示例代码如下:
#sklearn实现多元分类
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(C=100, random_state=1, solver='lbfgs', multi_class='ovr')
lr.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=lr, test_idx=range(105, 150))
plt.xlabel('petal legngth [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
运行代码后,得到下图:
对于分类结果,我们可以用predict_proba来计算训练数据集中某个样本属于某个特定类的概率。例如可以使用以下代码预测测试数据集中前三类的概率:
lr.predict_proba(X_test_std[:3, :1])
执行后返回结果:
array([[3.81527885e-09, 1.44792866e-01, 8.55207131e-01],
[8.34020679e-01, 1.65979321e-01, 3.25737138e-13],
[8.48831425e-01, 1.51168575e-01, 2.62277619e-14]])
以上数组中,每一行代表第几个样本的概率,所有列数据之和为1.而每列代表第几种花。例如,第一行的最大值为第三列的0.85,这意味着第一个样本属于第三种话(Iris-virginica)的预测概率为85%。为了快速找到每一行的最大值,可以使用numpy中的argmax函数实现,代码为:
lr.predict_proba(X_test_std[:3, :1].argmax(axis=1))
输出结果为:
array([2, 0, 0], dtype=int64)
这里其实是因为我们在前面计算了条件概率,直接用numpy中的argmax函数将其手动转换为分类标签。在scikit-learn中我们可以使用更为方便的predict方法,代码为:
lr.predict(X_test_std[:3, :1])
输出结果同上,这表示所分出的三个类别中,第一个样本类别被预测为类别2,即第三种花Iris-setosa;第二及第三个样本类别被预测为类别0,即第一种花Iris-virginica.
如果想要单独预测话样本的分类标签,可以在以上代码中输入一个二维数组,调用reshape方法增加一个新维度可以实现将一行数据转换为二维数组,代码如下:
lr.predict(X_test_std[0, :].reshape(1, -1))
此时输出结果为:
array([2])
这表示第一个样本被预测为类别2,即第三种花Iris-setosa。
3.3.5 通过正则化解决过拟合问题
- 过拟合:是指学习时选择的模型所包含的参数过多,以至于出现这一模型对已知数据预测的很好,但对未知数据预测得很差的现象。这种情况下模型可能只是记住了训练集数据,而不是学习到了数据特征。
- 欠拟合:模型描述能力太弱,以至于不能很好地学习到数据中的规律。产生欠拟合的原因通常是模型过于简单。
过拟合和欠拟合状态如下图所示:
找到一个好的偏置-方差权衡的方法是通过正则化来调整模型的复杂性。正则化是处理共线性、滤除数据中的噪声并最终防止过拟合的一种非常有效的方法。
正则化背后的逻辑是引入额外的信息(偏置)来惩罚极端的参数值(权重)。最常见的正则化则是所谓的L2正则化(有时也称作L2收缩或者权重衰减),可写作:
λ2∣∣ω∣∣2=λ2∑j=1mωj2\frac{\lambda}{2}||\omega||^2=\frac{\lambda}{2}\sum_{j=1}^m\omega_j^22λ∣∣ω∣∣2=2λj=1∑mωj2
其中λ\lambdaλ为正则化参数。
在其中增加一个简单的正则项,就可以正则化逻辑回归的代价函数:
J(ω)=∑i=1n[−y(i)log(ϕ(z(i)))−(1−y(i))log(1−ϕ(z(i)))]+λ2∣∣ω∣∣2J(\omega)=\sum_{i=1}^n[-y^{(i)}log(\phi(z^{(i)}))-(1-y^{(i)})log(1-\phi(z^{(i)}))]+\frac{\lambda}{2}||\omega||^2J(ω)=i=1∑n[−y(i)log(ϕ(z(i)))−(1−y(i))log(1−ϕ(z(i)))]+2λ∣∣ω∣∣2
这可以在模型训练的过程中缩小权重。通过正则化参数λ\lambdaλ,保持权重较小时,我们就可以控制模型与训练数据的拟合程度,加大λ\lambdaλ值可以增强正则化的强度。
上一节中用scikit-learn训练逻辑回归模型中的参数C就是正则化参数λ\lambdaλ的倒数,与λ\lambdaλ直接相关。因此,减小逆正则化参数C意味着增加正则化的强度,在这里,我们通过绘制对两个权重系数进行L2正则化后的图像来展示,代码为:
#正则化两个权重系数
weights, params = [], []
for c in np.arange(-5,5):
lr = LogisticRegression(C=10.**c, random_state=1, solver='lbfgs', multi_class='ovr')
lr.fit(X_train_std, y_train)
weights.append(lr.coef_[1])
params.append(10.**c)
weights = np.array(weights)
plt.plot(params, weights[:, 0], label='petal length')
plt.plot(params, weights[:, 1], linestyle='--', label='petal width')
plt.ylabel('weight coefficient')
plt.xlabel('C')
plt.legend(loc='upper left')
plt.xscale('log')
plt.show()
在以上代码中,我们用不同的逆正则化参数C值拟合10个逻辑回归模型,得到下图:
3.4 使用支持向量机最大化分类间隔
支持向量机(SVM) 是一种强大且广泛使用的机器学习算法,可以看作是感知器的扩展。感知器算法的目标是最小化分类误差,而支持向量机算法的优化目标是最大化分类间隔。所谓的间隔就是可分离的超平面(决策边界)和其最近的训练样本之间的距离。最靠近决策边界的训练样本就被称为支持向量。
3.4.1 对分类间隔最大化的直观认识

支持向量机的原理是:决策边界间隔较大往往会产生较低的泛化误差,而间隔较小的模型则更容易产生过拟合。
首先让我们来看看与决策边界平行的正超平面和负超平面,其数学表示为:
ω0+ωTxpos=1(1)ω0+ωTxneg=−1(2)\omega_0+\omega^Tx_{pos}=1\quad\quad\quad(1)\\ \omega_0+\omega^Tx_{neg}=-1\quad\quad(2)ω0+ωTxpos=1(1)ω0+ωTxneg=−1(2)
xposx_{pos}xpos(positive support vector):正类支持向量
xnegx_{neg}xneg(negative support vector):负类支持向量
将以上两个等式相减得到:
ωT(xpos−xneg)=2\omega^T(x_{pos}-x_{neg})=2ωT(xpos−xneg)=2
通过求解向量ω\omegaω的长度来规范化该方程,具体计算如下:
∣∣ω∣∣=∑j=1mωj2ωT(xpos−xneg)∣∣ω∣∣=2∣∣ω∣∣||\omega||=\displaystyle\sqrt{\sum_{j=1}^m\omega_j^2}\\ \displaystyle\frac{\omega^T(x_{pos}-x_{neg})}{||\omega||}=\frac{2}{||\omega||}∣∣ω∣∣=j=1∑mωj2∣∣ω∣∣ωT(xpos−xneg)=∣∣ω∣∣2
在最后的简化结果中,方程的左边即为正、负超平面之间的距离,也就是所谓的最大化间隔。
在样本分类正确的条件约束下,最大化分类间隔即2∣∣ω∣∣\displaystyle\frac{2}{||\omega||}∣∣ω∣∣2最大化,这就是支持向量机的目标函数,可以将其表示为:
ω0+ωTx(i)≥1,y(i)=1ω0+ωTx(i)≤−1,y(i)=−1\omega_0+\omega^Tx^{(i)}\geq1,\quad\quad\quad y^{(i)}=1\\ \omega_0+\omega^Tx^{(i)}\leq-1,\quad\quad y^{(i)}=-1ω0+ωTx(i)≥1,y(i)=1ω0+ωTx(i)≤−1,y(i)=−1
其中i=1……Ni=1……Ni=1……N.
这两个方程可以解释为:所有的负类样本基本上都落在负超平面一侧,所有的正类样本都落在正超平面一侧,我们可以使用更紧凑的方式表达为:
y(i)(ω0+ωTx(i))≥1∀iy^{(i)}(\omega_0+\omega^Tx^{(i)})\geq1\quad\quad\forall iy(i)(ω0+ωTx(i))≥1∀i
3.4.2 用松弛变量解决非线性可分问题
Vladimir Vapnik与1995年提出松弛变量ξ\xiξ,它引出了所谓的软间隔分类。对于非线性可分数据来说,需要放松线性约束以允许在分类错误存在的情况下通过适当代价的惩罚来确保优化可以收敛。
我们可以直接把取值为正的松弛变量加入线性约束,得到:
ω0+ωTx(i)≥1−ξ(i),y(i)=1ω0+ωTx(i)≤−1+ξ(i),y(i)=−1\omega_0+\omega^Tx^{(i)}\geq1-\xi^{(i)},\quad\quad\quad y^{(i)}=1\\ \omega_0+\omega^Tx^{(i)}\leq-1+\xi^{(i)},\quad\quad y^{(i)}=-1ω0+ωTx(i)≥1−ξ(i),y(i)=1ω0+ωTx(i)≤−1+ξ(i),y(i)=−1
这样我们就得到新的最小化(有约束)目标为:
12∣∣ω∣∣2+C(∑iξ(i))\frac{1}{2}||\omega||^2+C(\sum_i\xi^{(i)})21∣∣ω∣∣2+C(i∑ξ(i))
在这里,我们可以通过控制变量C来控制对分类错误的惩罚,C的值越大相应的错误惩罚就越大。如果选择的目标较小,对分类错误的要求就不那么严格。因此可以用参数C来控制间隔的宽度来权衡偏置-方差,如下图所示:
现在,我们就可以训练一个支持向量机模型来对鸢尾花数据集中的不同种类的花进行分类,代码如下:
#支持向量机实现鸢尾花分类
from sklearn.svm import SVC
svm = SVC(kernel='linear', C=1, random_state=1)
svm.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=svm, test_idx=range(105, 150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
执行代码得到以下图:
3.5 用核支持向量机求解非线性问题
支持向量机的优势在于,它可以很容易使用核技巧来解决非线性分类问题。
3.5.1 处理线性不可分数据的核方法
核技巧Kernel Trick是一种在支持向量机和其他基于核方法的机器学习算法中常用的技术。核技巧的目的是通过在高维特征空间中进行计算,而不需要显式地计算高维特征空间中的数据点的转换。
在传统的SVM中,通过将低维数据映射到高维特征空间中来解决数据线性不可分的问题,但是显式地进行高维转换需要计算大量的特征,这可能导致计算上的复杂性和存储开销的增加。核技巧通过利用核函数(kernel function)的特性,将原始低维空间中的数据点在高维特征空间中的内积表示为低维空间中数据点之间的核函数计算结果。这种技巧允许我们在计算时直接使用核函数,而不需要进行显式的高维特征转换,从而大大减少了计算上的复杂性。
在讨论核支持向量机的原理之前,我们先通过创建一个样本数据集来认识一下何为非线性分类问题。执行下述代码,调用numpy的logical_or函数创建一个经过“异或”操作的数据集,其中有100个样本的分类标签为1,100个样本的分类标签为-1:
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(1)
X_xor = np.random.randn(200, 2)
y_xor = np.logical_xor(X_xor[:, 0] > 0, X_xor[:, 1] > 0)
y_xor = np.where(y_xor, 1, -1)
plt.scatter(X_xor[y_xor == 1, 0], X_xor[y_xor == 1, 1], c='b', marker='x', label='1')
plt.scatter(X_xor[y_xor == -1, 0], X_xor[y_xor == -1, 1], c='r', marker='s', label='-1')
plt.xlim([-3, 3])
plt.ylim([-3, 3])
plt.legend(loc='best')
plt.tight_layout()
plt.show()
执行代码得到具有随机噪声的“异或”数据集,如下图所示:
显而易见,使用之前讨论的线性逻辑回归或线性支持向量机模型,并将线性超平面作为决策边界,是无法将上图中的样本正确地划分为正类和负类的。
核方法的基本思想是针对线性不可分数据,建立非线性组合,通过映射函数ϕ\phiϕ把原始特征投影到一个高维空间,使特征在该空间变得线性可分。
如上图所示,可以把一个二维数据集转换为一个新的三维特征空间,再通过投影
ϕ(x1,x2)=(z1,z2,z3)=(x1,x2,x12+x22)\phi(x_1,x_2)=(z_1,z_2,z_3)=(x_1,x_2,x_1^2+x_2^2)ϕ(x1,x2)=(z1,z2,z3)=(x1,x2,x12+x22)
使得样本可分。如上图所示,通过线性超平面,我们可以把图中两个类分开,如果再把它投影回原始特征空间上,就可以形成非线性的决策边界。
3.5.2 利用核技巧发现高维空间的分离超平面
为了使用SVM解决非线性问题,需要调用映射函数ϕ\phiϕ将训练数据变换到高维特征空间,然后训练线性SVM模型对新特征空间里的数据进行分类。可以用相同的映射函数ϕ\phiϕ对未知新数据进行变换,用线性支持向量机模型进行分类。
但这种映射方法构建新特征的计算成本太高,特别是在处理高维数据的时候。此时我们就可以使用核技巧来发挥作用。只要用ϕ(x(i))Tϕ(x(j))\phi(x^{(i)})^T\phi(x^{(j)})ϕ(x(i))Tϕ(x(j))替换点乘x(i)Tx(j)x^{(i)T}x^{(j)}x(i)Tx(j)。为显著降低计算两点间点成的昂贵计算成本,定义所谓的核函数如下:
κ(x(i)x(j))=ϕ(x(i))Tϕ(x(j))\kappa(x^{(i)}x^{(j)})=\phi(x^{(i)})^T\phi(x^{(j)})κ(x(i)x(j))=ϕ(x(i))Tϕ(x(j))
其中使用最为广泛的核函数就是径向基函数(RBF)核或者简称为高斯核:
κ(x(i),x(j))=e−∣∣x(i)−x(j)∣∣22σ2\kappa(x^{(i)},x^{(j)})=e^{-\displaystyle\frac{||x^{(i)}-x^{(j)}||^2}{2\sigma^2}}κ(x(i),x(j))=e−2σ2∣∣x(i)−x(j)∣∣2
简化得到:
κ(x(i),x(j))=e−γ∣∣x(i)−x(j)∣∣2\kappa(x^{(i)},x^{(j)})=e^{-\gamma\displaystyle||x^{(i)}-x^{(j)}||^2}κ(x(i),x(j))=e−γ∣∣x(i)−x(j)∣∣2
其中,γ=12σ2\gamma=\displaystyle\frac{1}{2\sigma^2}γ=2σ21是要优化的自由参数。
总之,所谓的“核”可以理解为一对样本之间的相似函数。公式中的负号把距离转换为相似性得分,而指数运算把由此产生的相似性得分控制在1(完全相似)核0(完全不同)之间。
现在我们来训练一个核支持向量机来对异或数据进行分类。这里只需要将SVC类中的参数kernel='rbf’替换掉kernel=‘linear’.代码如下:
svm = SVC(kernel='rbf', random_state=1, gamma=0.1, C=10)
svm.fit(X_xor, y_xor)
plot_decision_regions(X_xor, y_xor, classifier=svm)
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
运行代码,我们得到下图:
代码中的γ\gammaγ是高斯径向基函数的参数,用于控制模型的复杂度和决策边界的灵活性,控制了模型的泛化能力。
高斯径向基函数用于将输入特征映射到高维空间,从而使得数据在新的特征空间中更容易分离。γ\gammaγ参数控制了径向基函数的宽度.具体来说,较小的γ\gammaγ值会导致径向基函数更宽,决策边界较为平滑;而较大的γ\gammaγ值会导致径向基函数更窄,决策边界更为准确地围绕每个样本点。需要注意的是,γ\gammaγ参数的选择会对SVM的性能和预测结果产生重要影响。过小或过大的γ\gammaγ值都可能导致过拟合或欠拟合。通常情况下,通过交叉验证或网格搜索等方法来选择最佳的γ\gammaγ值,以获得更好的模型性能。
既然增大γ\gammaγ值可以增大训练样本的影响范围,从而导致决策边界紧缩和波动。为了能够更好地理解γ\gammaγ,我们把RBF核支持向量机应用于鸢尾花数据集:
svm = SVC(kernel='rbf', random_state=1, gamma=0.2, C=1)
svm.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=svm, test_idx=range(105, 150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
此时我们还是保守地选择一个较小地γ\gammaγ值,得到的RBF核SVM模型的决策边界相对偏松,可以看到如下图所示:
现在加大γ\gammaγ的值,观察其对决策边界的影响:
svm = SVC(kernel='rbf', random_state=1, gamma=100, C=1)
svm.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=svm, test_idx=range(105, 150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
此时我们得到图像:
可以看到,此时类别周围的决策边界更为紧密。虽然模型对训练数据集拟合效果非常好,但是这种分类器对未知数据可能会有一个很高的泛化误差,这说明当算法对训练数据集中的波动过于敏感时,参数γ\gammaγ在控制过拟合或方差方面也起着重要的作用。
3.6 决策树学习
决策树分类器在可解释性方面是非常有效的模型。我们也可以把决策树算法使用在鸢尾花分类之中。使用决策树算法,我们从树根开始,在信息增益最大的特征上分裂数据,各子节点在迭代过程中重复该分裂过程,直至只剩下叶子节点为止,这就意味着所有节点上的样本都属于同一类。在实践中,这可能会出现根深叶茂的树,这样会容易导致过拟合。因此我们通常希望限制树的最大深度来对树进行修剪。
划分数据集的之前之后信息发生的变化称为信息增益,它定义为一个特征能够为分类系统带来多少信息,带来的信息越多,说明该特征越重要,相应的信息增益也就越大。
决策树修剪是一种提高决策树模型泛化能力的方法,可以有效防止决策树过拟合的问题。
3.6.1 最大化信息增益————获得最大收益
为保证特征信息信息增益最大的情况下分裂节点,我们需要先定义目标函数,然后进行决策树学习和算法优化。定义以下目标函数以实现在每次分裂的时候最大化信息增益:
IG(Dp,f)=I(Dp)−∑j=1mNjNpI(DJ)IG(D_p, f)=I(D_p)-\sum_{j=1}^m\frac{N_j}{N_p}I(D_J)IG(Dp,f)=I(Dp)−j=1∑mNpNjI(DJ)
这里fff是分裂数据的特征依据,DpD_pDp和DjD_jDj为父节点和第jjj个子节点,III为杂质含量,NjN_jNj是第jjj个子节点的样本数。可以看到,父节点和子节点之间的信息增益仅仅在杂质含量方面存在着差异,即子节点杂质含量越低信息增益越大。
决策树模型中的杂质是指数据集中的不纯度或混杂程度。决策树算法通常使用一些评估指标来度量数据集的杂质,常见的指标包括:
- 熵(Entropy):衡量数据集的混乱程度。熵越高,表示数据集的混杂程度越高。
- 基尼不纯度(Gini impurity):衡量从数据集中随机选择两个样本,其类别标签不一致的概率。基尼不纯度越高,表示数据集的混杂程度越高。
- 分类误差(Classification Error):衡量将数据集中最常见的类别标签错误地分类为其他类别的概率。
为方便起见,同时考虑到减少组合搜索空间,大多数软件库只实现二元决策树,这意味着每个父节点只有DleftD_{left}Dleft和DrightD_{right}Dright两个子节点,即目标函数变为:
IG(Dp,f)=I(Dp)−NleftNpI(Dleft)−NrightNpI(Dright)IG(D_p, f)=I(D_p)-\frac{N_{left}}{N_p}I(D_{left})-\frac{N_{right}}{N_p}I(D_{right})IG(Dp,f)=I(Dp)−NpNleftI(Dleft)−NpNrightI(Dright)
在二元决策树中,质量杂质含量或者分裂标准的三个常用指标分别为基尼杂质度(IGI_GIG)、熵(IHI_HIH)和分类误差(IEI_EIE)。
其中,非空类(p(i∣t)log2p(i∣t)p(i|t)log_2p(i|t)p(i∣t)log2p(i∣t))熵定义为:
IH(t)=−∑i=1cp(i∣t)log2p(i∣t)I_H(t)=-\displaystyle\sum_{i=1}^cp(i|t)log_2p(i|t)IH(t)=−i=1∑cp(i∣t)log2p(i∣t)
其中p(i∣t)≠0p(i|t)\neq0p(i∣t)=0为某节点t属于i类样本的概率。如果节点上的所有样本都属于同一个类,则熵为0;如果类的分布均匀,则熵值最大。例如,对于二元分类,如果p(i=1∣t)=1p(i=1|t)=1p(i=1∣t)=1或p(i=0∣t)=0p(i=0|t)=0p(i=0∣t)=0,则熵为0.如果类分布均匀,p(i=1∣t)=0.5p(i=1|t)=0.5p(i=1∣t)=0.5,p(i=0∣t)=0.5p(i=0|t)=0.5p(i=0∣t)=0.5,则熵为1.因此可以说熵准则试图减少错误分类概率的判断标准。
而我们可以把基尼杂质理解为尽量减少错误分类概率的判断标准,定义为:
IG(t)=∑i=1cp(i∣t)(1−p(i∣t))=1−∑i=1cp(i∣t)2I_G(t)=\displaystyle\sum_{i=1}^cp(i|t)(1-p(i|t))=1-\sum_{i=1}^cp(i|t)^2IG(t)=i=1∑cp(i∣t)(1−p(i∣t))=1−i=1∑cp(i∣t)2
与熵类似,如果类是完全混合的,那么基尼杂质最大,例如在二元分类中(c=2),此时有:
IG(t)=1−∑i=1c0.52=0.5I_G(t)=1-\sum_{i=1}^c0.5^2=0.5IG(t)=1−i=1∑c0.52=0.5
然而,基尼杂质和熵在实践中经常会产生非常相似的结果,而且通常并不值得花很多时间用不同的杂质标准来评估树,更好的选择是尝试不同的修剪方法。
分类误差定义为:
IE(t)=1−max{p(i∣t)}I_E(t)=1-max\{p(i|t)\}IE(t)=1−max{p(i∣t)}
这对修剪树枝来说是个非常有用的判断标准,但在决策树中并不推荐,因为它对节点的分类概率变化不太敏感。可以从下图所示的两种可能的分裂场景来理解。
从数据集的父节点DpD_pDp开始,它包含40个分类标签为1的样本和40个分类标签为2的样本,要分裂成DleftD_{left}Dleft和DrightD_{right}Dright两个数据集。在A和B两种场景下,用分类误差作为分裂标准得到信息增益(IGE=0.25IG_E=0.25IGE=0.25)相同,计算过程如下:
IE(Dp)=1−0.5=0.5A:IE(Dleft)=1−34=0.25A:IE(Dright)=1−34=0.25A:IGE=0.5−480.25−480.25=0.25B:IE(Dleft)=1−46=13B:IE(Dright)=1−1=0B:IGE=0.5−68×13−0=0.25I_E(D_p)=1-0.5=0.5\\ A: I_E(D_{left})=1-\frac{3}{4}=0.25\\ A: I_E(D_{right})=1-\frac{3}{4}=0.25\\ A: IG_E=0.5-\frac{4}{8}0.25-\frac{4}{8}0.25=0.25\\ B: I_E(D_{left})=1-\frac{4}{6}=\frac{1}{3}\\ B: I_E(D_{right})=1-1=0\\ B:IG_E=0.5-\frac{6}{8}\times\frac{1}{3}-0=0.25IE(Dp)=1−0.5=0.5A:IE(Dleft)=1−43=0.25A:IE(Dright)=1−43=0.25A:IGE=0.5−840.25−840.25=0.25B:IE(Dleft)=1−64=31B:IE(Dright)=1−1=0B:IGE=0.5−86×31−0=0.25
然而,与场景A(IGG=0.125IG_G=0.125IGG=0.125)相比,基尼杂质有利于分裂场景B(IGG=0.16IG_G=0.16IGG=0.16),该场景确实是纯度更高,计算过程如下:
IG(Dp)=1−(0.52+0.52)=0.5A:IG(Dleft)=1−((34)2+(14)2)=38=0.375A:IG(Dright)=1−((14)2+(34)2)=38=0.375A:IGG=0.5−480.375−480.375=0.125B:IG(Dleft)=1−((26)2+(46)2)=49=0.4B:IG(Dright)=1−(12+02)=0B:IGG=0.5−680.4−0=0.16I_G(D_p)=1-(0.5^2+0.5^2)=0.5\\ A: I_G(D_{left})=1-((\frac{3}{4})^2+(\frac{1}{4})^2)=\frac{3}{8}=0.375\\ A: I_G(D_{right})=1-((\frac{1}{4})^2+(\frac{3}{4})^2)=\frac{3}{8}=0.375\\ A: IG_G=0.5-\frac{4}{8}0.375-\frac{4}{8}0.375=0.125\\ B: I_G(D_{left})=1-((\frac{2}{6})^2+(\frac{4}{6})^2)=\frac{4}{9}=0.4\\ B: I_G(D_{right})=1-(1^2+0^2)=0\\ B: IG_G=0.5-\frac{6}{8}0.4-0=0.16\\ IG(Dp)=1−(0.52+0.52)=0.5A:IG(Dleft)=1−((43)2+(41)2)=83=0.375A:IG(Dright)=1−((41)2+(43)2)=83=0.375A:IGG=0.5−840.375−840.375=0.125B:IG(Dleft)=1−((62)2+(64)2)=94=0.4B:IG(Dright)=1−(12+02)=0B:IGG=0.5−860.4−0=0.16
同样,与场景A(IGH=0.19IG_H=0.19IGH=0.19)相比,熵准则对场景B(IG_H=0.31)更为有利,计算过程如下:
IH(Dp)=−(0.5log2(0.5)+0.5log2(0.5))=1A:IH(Dleft)=−(34log2(34)+14log2(14))=0.81A:IH(Dright)=−(14log2(14)+34log2(34))=0.81A:IGH=1−480.81−480.81=0.19B:IH(Dleft)=−(26log2(26)+46log2(46))=0.92B:IH(Dright)=0B:IGH=1−680.92−0=0.31I_H(D_p)=-(0.5log_2(0.5)+0.5log_2(0.5))=1\\ A: I_H(D_{left})=-(\frac{3}{4}log_2(\frac{3}{4})+\frac{1}{4}log_2(\frac{1}{4}))=0.81\\ A: I_H(D_{right})=-(\frac{1}{4}log_2(\frac{1}{4})+\frac{3}{4}log_2(\frac{3}{4}))=0.81\\ A: IG_H=1-\frac{4}{8}0.81-\frac{4}{8}0.81=0.19\\ B: I_H(D_{left})=-(\frac{2}{6}log_2(\frac{2}{6})+\frac{4}{6}log_2(\frac{4}{6}))=0.92\\ B: I_H(D_{right})=0\\ B: IG_H=1-\frac{6}{8}0.92-0=0.31IH(Dp)=−(0.5log2(0.5)+0.5log2(0.5))=1A:IH(Dleft)=−(43log2(43)+41log2(41))=0.81A:IH(Dright)=−(41log2(41)+43log2(43))=0.81A:IGH=1−840.81−840.81=0.19B:IH(Dleft)=−(62log2(62)+64log2(64))=0.92B:IH(Dright)=0B:IGH=1−860.92−0=0.31
为了能更好地比较这三种不同的杂质标准,我们将把分类标签为1、概率在[0, 1]之间的杂质情况画在图上展示出来。在这里,我们添加一个小比例样本的熵(熵/2)来观察介于熵和分类误差中间的度量标准即基尼杂质,具体的代码示例如下:
import matplotlib.pyplot as plt
import numpy as np
def gini(p):
return p * (1 - p) + (1 - p) * (1 - (1 - p))
def entropy(p):
return - p * np.log2(p) - (1 - p) * np.log2((1 - p))
def error(p):
return 1 - np.max([p, 1 - p])
x = np.arange(0.0, 1.0, 0.01)
ent = [entropy(p) if p != 0 else None for p in x]
sc_ent = [e * 0.5 if e else None for e in ent]
err = [error(i) for i in x]
fig = plt.figure()
ax = plt.subplot(111)
for i, lab, ls, c, in zip([ent, sc_ent, gini(x), err],
['Entropy', 'Entropy (scaled)',
'Gini impurity', 'Misclassification error'],
['-', '-', '--', '-.'],
['black', 'lightgray', 'red', 'green', 'cyan']):
line = ax.plot(x, i, label=lab, linestyle=ls, lw=2, color=c)
ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),
ncol=5, fancybox=True, shadow=False)
ax.axhline(y=0.5, linewidth=1, color='k', linestyle='--')
ax.axhline(y=1.0, linewidth=1, color='k', linestyle='--')
plt.ylim([0, 1.1])
plt.xlabel('p(i=1)')
plt.ylabel('impurity index')
plt.show()
执行代码得到下图:
3.6.2 构建决策树
通过将特征空间划分成不同的举行,决策树可以构建复杂的决策边界。但需要注意决策树越深,决策边界就越复杂,越容易导致过拟合。假设最大深度为4,以熵作为杂质度量标准,用scikit-learn来训练决策树模型,代码如下(注意,此处为了可视化进行了特征缩放,此特征缩放并非决策树算法的要求):
from sklearn.tree import DecisionTreeClassifier
tree_model = DecisionTreeClassifier(criterion='gini',
max_depth=4,
random_state=1)
tree_model.fit(X_train, y_train)
X_combined = np.vstack((X_train, X_test))
y_combined = np.hstack((y_train, y_test))
plot_decision_regions(X_combined, y_combined,
classifier=tree_model,
test_idx=range(105, 150))
plt.xlabel('petal length [cm]')
plt.ylabel('petal width [cm]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
执行代码得到决策边界如下图:
非常好的是,scikit-learn允许我们在训练完成后实现决策树模型可视化,实现代码如下:
from sklearn import tree
tree.plot_tree(tree_model)
plt.show()
执行后轻松得到:
我们也可以使用PyDotPlus库来把.dot数据文件转换为决策树图像,画出的图像将更好看,在执行代码前,需要先在电脑命令窗口安装相关库:
pip install graphviz
pip install pyparsing
pip install pydotplus
接下来执行以下代码以在本地目录创建决策树的图像文件:
from pydotplus import graph_from_dot_data
from sklearn.tree import export_graphviz
dot_data = export_graphviz(tree_model,
filled=True,
rounded=True,
class_names=['Setosa''Versicolor''Virginica'],
feature_names=['petal length','petal width'],
out_file=None)
graph = graph_from_dot_data(dot_data)
graph.write_png('tree.png')
通过设置out_fit=None,可以直接把数据赋予dot_data变量,而不是在磁盘上产生中间文件tree.dot。参数filled、rounded、class_names和feature_names可以设置颜色、边框和边缘圆角,并在每个节点上显示大多数分类标签以及分裂标准的特征。执行代码得到下图:
在这个决策图上我们可以很好地追溯训练数据集的分裂过程。从105个样本的根节点开始,以花瓣宽度0.75厘米作为截至条件,先把样本数据分裂成大小分别为35和70的两个子节点,可以看到左边子节点的纯度已经很高,它只包含Itis-setosa类,而右边子节点则进一步分裂成Iris-versicolor和Iris-virginica两类。
在树和决策区域图上,我们可以看到决策树在花朵分类上的表现不错,但scikit-learn目前还不提供手工修剪决策树的函数。
3.6.3 多个决策树的随机森林组合
基于决策树的随机森林算法以其良好的可拓展性和易用性而闻名,可以把随机森林视为决策树的集成。随机森林的原理是对受较大方差影响的多个决策树分别求平均值,然后构建一个具有更好泛化性能和不易过拟合的强大模型。随即森林可以概括为以下四个步骤:
- 1)随即提取一个规模为n的bootstrap样本(从训练数据集中有放回地随机选取n个样本);
- 2)基于bootstrap样本生成随机数同时:a.不放回地随机选择d个特征;b.根据目标函数的要求,例如信息增益最大化,使用选定的最佳特征来分裂节点;
- 3)把以上两个步骤重复k次;
- 4)聚合每棵树的预测结果,并以多数票机制确定标签的分类。
Bootstrap:其基本思想是对现有的数据,不断再随机取小的样本,对每个小样处理数据,得到需要估计的统计量.从而来了解估计统计量的分布。它是一种重新取样的方法,从现有的样本数据中独立取样,并替换相同数量的样本,在这些重新取样的数据中进行推断。
需要注意的是,在步骤2中训练单个决策树的时候,不需要评估所有的特征来确定节点的最佳分裂方案,而是仅仅考虑其中的一个随机子集。
虽然随机森林的可解释性不如决策树,但是它的显著优势在于不必担心超参数值的选择。通常不需要修剪随机森林,因为集成模型对来自单个决策树的噪声有较强的抵抗力。唯一需要关心的参数是在步骤3中如何选择随机森林中的树的数量k,一般而言,树越多,随机森林分类器的性能就越好,当然计算成本也会随之增大。
事实上随机森林分类器的其他超参数也可以被优化,可以通过bootstrap样本规模n来控制随机森林的偏差-方差权衡。因为特定训练数据被包含在bootstrap样本中的概率比较低,所以减小bootstrap样本的规模会增加在单个树之间的多样性。因此,缩小bootstrap样本的规模可能会增加随机森林的随机性,这有助于降低过拟合。然而,较小规模的bootstrap样本通常也会导致随机森林的总体性能较差,训练样本和测试样本之间的性能差别很小,但总体测试性能较低。相反,扩大bootstrap样本的规模可能会增加过拟合。由于bootstrap样本的存在,各个决策树之间会变得更相似,它们可以通过学习更加紧密地拟合原始数据的训练数据集。
包括scikit-learn中的RandomForestClassifier实现在内,在大多数的随机森林实现中,bootstrap样本的规模一般与原始训练数据集样本的规模保持一致,这通常会有一个好的偏置-方差权衡。对于每轮分裂的特征数d,我们希望选择的值小于训练数据集的特征总数。在scikit-learn以及其他实现中,合理的默认值为d=md=\sqrt{m}d=m,这里m为训练数据集的特征总数。
在scikit-learn当中已经实现了一个分类器供我们使用,代码如下:
from sklearn.ensemble import RandomForestClassifier
forest = RandomForestClassifier(criterion='gini',n_estimators=25,random_state=1,n_jobs=2)
forest.fit(X_train, y_train)
plot_decision_regions(X_combined, y_combined, classifier=forest, test_idx=range(105,150))
plt.xlabel('petal length[cm]')
plt.ylabel('petal width[cm]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
执行代码后得到由随机森林中的树继集成产生的决策区域,如下图所示:
3.7 K近邻——一种惰性学习算法
K近邻(KNN)分类器和我们讨论过的其他分类器都有所不同,它是惰性学习中的典型例子。所谓的惰性,并不是说它看上去很简单,而是在于它不是从训练数据中学习判别函数,而是靠记忆训练过的数据集来完成任务。
参数模型和非参数模型
机器学习算法可以分为参数模型和非参数模型。
- 参数模型指我们从训练数据集估计参数来学习,一个能分类新数据点的函数,而不再需要原始训练数据集。参数模型的典型例子是感知器、逻辑回归和线性支持向量机。
- 非参数模型就不能用一组固定的参数来描述,参数的个数随着训练数据的增加而增加,非参数模型的例子就是决策树分类器、随机森林和核支持向量机。
KNN算法是基于实例的学习,属于非参数模型。基于实例学习的模型以记忆训练数据集为特征,惰性学习是基于实例学习的一种特殊情况,这与它在学习过程中付出的零代价有关。
KNN算法可以总结为以下几个步骤:
- 选择k个数和一个距离度量;
- 找到要分类样本的k近邻;
- 以多数票机制确定分类标签。

上图显示了在五个近邻中,如何基于多数票机制为新数据点(?)分配三角形标签。
基于所选择的距离度量,KNN算法从训练数据集中找到最接近要分类新数据点的k个样本,新数据点的分类标签由最靠近该点的k个数据点的多数票决定。
基于记忆方法的主要优点是当新的训练数据出现时,分类器可以立即适应。但缺点是新样本分类的计算复杂度与最坏情况下训练数集的样本数呈线性关系,除非数据集只有很少的维度(特征),而且算法实现采用有效的数据结构(如KD树)。此外,由于没有训练步骤,因此不能丢弃训练样本。因此如果要处理大型数据集,存储空间将面临挑战。
执行下列代码可以用scikit-learn的欧氏距离度量来实现KNN模型:
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=5, p=2, metric='minkowski')
knn.fit(X_train, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=knn, test_idx=range(105,150))
plt.xlabel('petal length[standardized]')
plt.ylabel('petal width[standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
在这里使用KNN模型指定5个邻居,得到以下决策边界为:
在争执不下的情况下,KNN算法更喜欢近距离样本的邻居。如果邻居有相似距离,那么该算法将在训练数据集中选择最先出现的分类标签。
正确选择k值对在过拟合和欠拟合之间找到恰当的平衡至关重要,我们还必须确保选择的距离度量适合数据集中的特征。通常用简单的欧氏距离来度量。但是如果使用欧氏距离的话,就要确保对数据进行了正确的标准化处理。
在我们上面的代码中,使用的是minkowski距离,它是欧氏距离和曼哈顿距离的推广,公式为
d(x(i),x(j))=∑k∣xk(i)xk(j)∣ppd(x^{(i)},x^{(j)})=\sqrt[p]{\sum_k|x_k^{(i)}\quad\quad x_k^{(j)}|^p}d(x(i),x(j))=pk∑∣xk(i)xk(j)∣p
在参数p中,如果p=2就是指欧氏距离,p=1就是指曼哈顿距离。
附全文代码:
#把鸢尾花数据存入特征矩阵
from sklearn import datasets
import numpy as np
iris = datasets.load_iris()
X = iris.data[:, [2, 3]]
y = iris.target
print('Class labels:',np.unique(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.3, random_state=1, stratify=y)
#定义stratify=y获得内置分层支持
print('Labels counts in y:', np.bincount(y))
print('Lbels counts in y_train:', np.bincount(y_train))
print('Lbels counts in y_test:', np.bincount(y_test))
#标准化特征
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)
#使用OvR训练感知器
from sklearn.linear_model import Perceptron
ppn = Perceptron(eta0=0.1, random_state=1)
ppn.fit(X_train_std, y_train)
#预测
y_pred = ppn.predict(X_test_std)
print('Misclassified examples: %d' % (y_test != y_pred).sum())
#调用metrics实现分类准确率
from sklearn.metrics import accuracy_score
print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))
#调用predict和accuracy_score实现分类准确率
print('Accuracy: %.3f' % ppn.score(X_test_std, y_test))
#绘制决策区
from matplotlib.colors import ListedColormap
import matplotlib.pyplot as plt
def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):
#setup marker generator and color map
markers = ('s', 'x', 'o', '^', 'v')
colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
cmap = ListedColormap(colors[:len(np.unique(y))])
#plot the decision surface
x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
np.arange(x2_min, x2_max, resolution))
Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
Z = Z.reshape(xx1.shape)
plt.contourf(xx1, xx2, Z, alpha=0.3, cmap=cmap)
plt.xlim(xx1.min(), xx1.max())
plt.ylim(xx2.min(), xx2.max())
for idx, cl in enumerate(np.unique(y)):
plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=colors[idx], marker=markers[idx], label=cl, edgecolor='black')
#highlight test examples
if test_idx:
#plot all examples
X_test, y_test = X[test_idx, :], y[test_idx]
plt.scatter(X_test[:, 0], X_test[:, 1], facecolors='none', edgecolor='black', alpha=1, linewidth=1, marker='o', s=100, label='test set')
#修改plot_decision_regions函数
X_combined_std = np.vstack((X_train_std, X_test_std))
y_combined = np.hstack((y_train, y_test))
plot_decision_regions(X=X_combined_std, y=y_combined, classifier=ppn, test_idx=range(105, 150))
plt.xlabel('petal length[standardized]')
plt.ylabel('petal width[standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#Logistic
import numpy as np
import matplotlib.pyplot as plt
# 定义Logistic激活函数
def logistic(z):
return 1.0 / (1.0 + np.exp(-z))
# 生成输入值范围
z = np.arange(-7, 7, 0.1)
phi_z = logistic(z)
plt.plot(z, phi_z)
plt.axvline(0, color='k')
plt.ylim(-0.1, 1.1)
plt.xlabel('z')
plt.ylabel('$\phi(z)$')
plt.yticks([0, 0.5, 1])
ax = plt.gca()
ax.yaxis.grid(True)
plt.tight_layout()
plt.show()
#分类训练样本的代价
def cost_1(z):
return - np.log(logistic(z))
def cost_0(z):
return - np.log(1 - logistic(z))
z = np.arange(-10, 10, 0.1)
phi_z = logistic(z)
c1 = [cost_1(x) for x in z]
plt.plot(phi_z, c1, label='J(w) if y=1')
c0 = [cost_0(x) for x in z]
plt.plot(phi_z, c0, linestyle='--', label='J(w) if y=0')
plt.ylim(0, 5.1)
plt.xlim([0,1])
plt.xlabel('$\phi(z)$')
plt.ylabel('J(w)')
plt.legend(loc='best')
plt.tight_layout()
plt.show()
#Adaline转换为逻辑回归
class LogisticRegressionGD(object):
"""Logistic Regression Classifier using gradient descent.
Parameters
------------
eta : float
Learning rate (between 0.0 and 1.0)
n_iter : int
Passes over the training dataset.
random_state : int
Random number generator seed for random weight
initialization.
Attributes
-----------
w_ : 1d-array
Weights after fitting.
cost_ : list
Logistic cost function value in each epoch.
"""
def __init__(self, eta=0.05, n_iter=100, random_state=1):
self.eta = eta
self.n_iter = n_iter
self.random_state = random_state
def fit(self, X, y):
""" Fit training data.
Parameters
----------
X : {array-like}, shape = [n_examples, n_features]
Training vectors, where n_examples is the number of examples and
n_features is the number of features.
y : array-like, shape = [n_examples]
Target values.
Returns
-------
self : object
"""
rgen = np.random.RandomState(self.random_state)
self.w_ = rgen.normal(loc=0.0, scale=0.01, size=1 + X.shape[1])
self.cost_ = []
for i in range(self.n_iter):
net_input = self.net_input(X)
output = self.activation(net_input)
errors = (y - output)
self.w_[1:] += self.eta * X.T.dot(errors)
self.w_[0] += self.eta * errors.sum()
# note that we compute the logistic `cost` now
# instead of the sum of squared errors cost
cost = -y.dot(np.log(output)) - ((1 - y).dot(np.log(1 - output)))
self.cost_.append(cost)
return self
def net_input(self, X):
"""Calculate net input"""
return np.dot(X, self.w_[1:]) + self.w_[0]
def activation(self, z):
"""Compute logistic sigmoid activation"""
return 1. / (1. + np.exp(-np.clip(z, -250, 250)))
def predict(self, X):
"""Return class label after unit step"""
return np.where(self.net_input(X) >= 0.0, 1, 0)
#只考虑两种花,验证逻辑回归的有效性
X_train_01_subset = X_train[(y_train == 0) | (y_train == 1)]
y_train_01_subset = y_train[(y_train == 1) | (y_train == 1)]
lrgd = LogisticRegressionGD(eta=0.05, n_iter=1000, random_state=1)
lrgd.fit(X_train_01_subset, y_train_01_subset)
plot_decision_regions(X=X_train_01_subset, y=y_train_01_subset, classifier=lrgd)
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#sklearn实现多元分类
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(C=100, random_state=1, solver='lbfgs', multi_class='ovr')
lr.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=lr, test_idx=range(105, 150))
plt.xlabel('petal legngth [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#预测数据集中前三类的概率:
lr.predict_proba(X_test_std[:3, :1])
#找到最大值
lr.predict_proba(X_test_std[:3, :1].argmax(axis=1))
#sklearn-predict方法找最大值
lr.predict(X_test_std[:3, :1])
#只预测一种分类标签
lr.predict(X_test_std[0, :].reshape(1, -1))
#正则化两个权重系数
weights, params = [], []
for c in np.arange(-5,5):
lr = LogisticRegression(C=10.**c, random_state=1, solver='lbfgs', multi_class='ovr')
lr.fit(X_train_std, y_train)
weights.append(lr.coef_[1])
params.append(10.**c)
weights = np.array(weights)
plt.plot(params, weights[:, 0], label='petal length')
plt.plot(params, weights[:, 1], linestyle='--', label='petal width')
plt.ylabel('weight coefficient')
plt.xlabel('C')
plt.legend(loc='upper left')
plt.xscale('log')
plt.show()
#支持向量机实现鸢尾花分类
from sklearn.svm import SVC
svm = SVC(kernel='linear', C=1, random_state=1)
svm.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=svm, test_idx=range(105, 150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#创建异或操作后的数据集
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(1)
X_xor = np.random.randn(200, 2)
y_xor = np.logical_xor(X_xor[:, 0] > 0, X_xor[:, 1] > 0)
y_xor = np.where(y_xor, 1, -1)
plt.scatter(X_xor[y_xor == 1, 0], X_xor[y_xor == 1, 1], c='b', marker='x', label='1')
plt.scatter(X_xor[y_xor == -1, 0], X_xor[y_xor == -1, 1], c='r', marker='s', label='-1')
plt.xlim([-3, 3])
plt.ylim([-3, 3])
plt.legend(loc='best')
plt.tight_layout()
plt.show()
#训练核支持向量机
svm = SVC(kernel='rbf', random_state=1, gamma=0.1, C=10)
svm.fit(X_xor, y_xor)
plot_decision_regions(X_xor, y_xor, classifier=svm)
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#使用核支持向量机到鸢尾花数据集
svm = SVC(kernel='rbf', random_state=1, gamma=0.2, C=1)
svm.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=svm, test_idx=range(105, 150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#加大gamma的值
svm = SVC(kernel='rbf', random_state=1, gamma=100, C=1)
svm.fit(X_train_std, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=svm, test_idx=range(105, 150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#比较三种杂质标准
import matplotlib.pyplot as plt
import numpy as np
def gini(p):
return p * (1 - p) + (1 - p) * (1 - (1 - p))
def entropy(p):
return - p * np.log2(p) - (1 - p) * np.log2((1 - p))
def error(p):
return 1 - np.max([p, 1 - p])
x = np.arange(0.0, 1.0, 0.01)
ent = [entropy(p) if p != 0 else None for p in x]
sc_ent = [e * 0.5 if e else None for e in ent]
err = [error(i) for i in x]
fig = plt.figure()
ax = plt.subplot(111)
for i, lab, ls, c, in zip([ent, sc_ent, gini(x), err],
['Entropy', 'Entropy (scaled)',
'Gini impurity', 'Misclassification error'],
['-', '-', '--', '-.'],
['black', 'lightgray', 'red', 'green', 'cyan']):
line = ax.plot(x, i, label=lab, linestyle=ls, lw=2, color=c)
ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15),
ncol=5, fancybox=True, shadow=False)
ax.axhline(y=0.5, linewidth=1, color='k', linestyle='--')
ax.axhline(y=1.0, linewidth=1, color='k', linestyle='--')
plt.ylim([0, 1.1])
plt.xlabel('p(i=1)')
plt.ylabel('impurity index')
plt.show()
#构建决策树
from sklearn.tree import DecisionTreeClassifier
tree_model = DecisionTreeClassifier(criterion='gini',
max_depth=4,
random_state=1)
tree_model.fit(X_train, y_train)
X_combined = np.vstack((X_train, X_test))
y_combined = np.hstack((y_train, y_test))
plot_decision_regions(X_combined, y_combined,
classifier=tree_model,
test_idx=range(105, 150))
plt.xlabel('petal length [cm]')
plt.ylabel('petal width [cm]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#决策树可视化
from sklearn import tree
tree.plot_tree(tree_model)
plt.show()
#使用pydotplus绘制决策树
from pydotplus import graph_from_dot_data
from sklearn.tree import export_graphviz
dot_data = export_graphviz(tree_model,
filled=True,
rounded=True,
class_names=['Setosa',
'Versicolor',
'Virginica'],
feature_names=['petal length',
'petal width'],
out_file=None)
graph = graph_from_dot_data(dot_data)
graph.write_png('tree.png')
#sklearn中的分类器
from sklearn.ensemble import RandomForestClassifier
forest = RandomForestClassifier(criterion='gini',n_estimators=25,random_state=1,n_jobs=2)
forest.fit(X_train, y_train)
plot_decision_regions(X_combined, y_combined, classifier=forest, test_idx=range(105,150))
plt.xlabel('petal length[cm]')
plt.ylabel('petal width[cm]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
#实现KNN模型
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=5, p=2, metric='minkowski')
knn.fit(X_train, y_train)
plot_decision_regions(X_combined_std, y_combined, classifier=knn, test_idx=range(105,150))
plt.xlabel('petal length[standardized]')
plt.ylabel('petal width[standardized]')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
更多推荐




所有评论(0)