Python-数据科学秘籍-四-

Python 数据科学秘籍（四）

原文：annas-archive.org/md5/a4f348a4e11e27ea41410c793e63daff

译者：飞龙

协议：CC BY-NC-SA 4.0

第七章 机器学习 2

在本章中，我们将介绍以下几种方法：

• 使用回归预测实值

• 学习带 L2 收缩的回归——岭回归

• 学习带 L1 收缩的回归——LASSO

• 使用交叉验证迭代器与 L1 和 L2 收缩

介绍

在本章中，我们将介绍回归技术及其在 Python 中的实现方法。接着，我们将讨论回归方法固有的一些缺点，并探讨如何通过收缩方法来解决这些问题。收缩方法中需要设置一些参数。我们将讨论交叉验证技术，以找到收缩方法的最佳参数值。

在上一章中，我们讨论了分类问题。在本章中，让我们将注意力转向回归问题。在分类问题中，响应变量Y要么是二元的，要么是一组离散值（在多类和多标签问题中）。而在回归中，响应变量是一个实值数字。

回归可以看作是一种函数逼近。回归的任务是找到一个函数，使得当X（一组随机变量）作为输入传递给该函数时，它应该返回Y，即响应变量。X也被称为自变量，而Y则被称为因变量。

我们将利用上一章中学到的技术，将数据集分为训练集、开发集和测试集，在训练集上迭代地构建模型，并在开发集上进行验证。最后，我们将使用测试集来全面评估模型的优劣。

本章将从使用最小二乘估计的简单线性回归方法开始。在第一种方法的开头，我们将简要介绍回归的框架，这是理解本章其他方法所必需的基础背景信息。尽管非常强大，简单回归框架也存在一个缺点。由于无法控制线性回归系数的上下限，它们往往会过拟合给定的数据。（线性回归的代价函数是无约束的，我们将在第一种方法中进一步讨论）。输出的回归模型可能在未见过的数据集上表现不佳。为了应对这个问题，使用了收缩方法。收缩方法也被称为正则化方法。在接下来的两种方法中，我们将介绍两种不同的收缩方法——LASSO 和岭回归。在最后一种方法中，我们将介绍交叉验证的概念，并展示如何利用它来估计传递给岭回归（作为一种收缩方法）的参数 alpha。

使用回归预测实值

在我们深入探讨此方案之前，先快速了解回归一般是如何运作的。这个介绍对于理解当前方案以及后续方案非常重要。

回归是一种特殊形式的函数逼近。以下是预测变量集：

每个实例xi具有m个属性：

回归的任务是找到一个函数，使得当X作为输入提供给该函数时，应该返回一个Y响应变量。Y是一个实值向量：

我们将使用波士顿房价数据集来解释回归框架。

以下链接提供了波士顿房价数据集的良好介绍：

archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names。

在本例中，响应变量Y是波士顿地区自有住房的中位数价值。有 13 个预测变量。前面的网页链接提供了所有预测变量的良好描述。回归问题的定义是找到一个函数F，使得如果我们给这个函数一个以前未见过的预测值，它应该能够给出中位数房价。

函数F(X)是我们线性回归模型的输出，它是输入X的线性组合，因此得名线性回归：

wi变量是前述方程中的未知值。建模过程的关键在于发现wi变量。通过使用我们的训练数据，我们将找到wi的值；wi被称为回归模型的系数。

线性回归建模问题可以表述为：使用训练数据来找到系数：

使得：

上述公式值尽可能小。

该方程的值越低（在优化术语中称为代价函数），线性回归模型越好。因此，优化问题就是最小化前述方程，即找到wi系数的值，使其最小化该方程。我们不会深入探讨优化算法的细节，但了解这个目标函数是有必要的，因为接下来的两个方案需要你理解它。

现在，问题是我们如何知道使用训练数据构建的模型，即我们新找到的系数w1, w2,..wm，是否足够准确地预测未见过的记录？我们将再次利用之前定义的成本函数。当我们将模型应用于开发集或测试集时，我们找到实际值和预测值之间差异平方的平均值，如下所示：

上述方程被称为均方误差（mean squared error）——这是我们用来判断回归模型是否值得使用的标准。我们希望得到一个输出模型，使得实际值和预测值之间的差异平方的平均值尽可能低。寻找系数的这个方法被称为最小二乘估计。

我们将使用 scikit-learn 的LinearRegression类。然而，它在内部使用scipy.linalg.lstsq方法。最小二乘法为我们提供了回归问题的闭式解。有关最小二乘法及其推导的更多信息，请参考以下链接：

en.wikipedia.org/wiki/Least_squares。

en.wikipedia.org/wiki/Linear_least_squares_(mathematics)。

我们对回归做了一个非常简单的介绍。感兴趣的读者可以参考以下书籍：www.amazon.com/exec/obidos/ASIN/0387952845/trevorhastie-20。

www.amazon.com/Neural-Networks-Learning-Machines-Edition/dp/0131471392

准备就绪

波士顿数据集有 13 个属性和 506 个实例。目标变量是一个实数，房屋的中位数值在千美元左右。

参考以下 UCI 链接了解有关波士顿数据集的更多信息：archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names。

我们将提供这些预测变量和响应变量的名称，如下所示：

如何实现…

我们将从加载所有必要的库开始。接下来，我们将定义第一个函数get_data()。在这个函数中，我们将读取波士顿数据集，并将其作为预测变量x和响应变量y返回：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

在我们的build_model函数中，我们将使用给定的数据构建线性回归模型。以下两个函数，view_model和model_worth，用于检查我们所构建的模型：

def build_model(x,y):
    """
    Build a linear regression model
    """
    model = LinearRegression(normalize=True,fit_intercept=True)
    model.fit(x,y)
    return model    

def view_model(model):
    """
    Look at model coeffiecients
    """
    print "\n Model coeffiecients"
    print "======================\n"
    for i,coef in enumerate(model.coef_):
        print "\tCoefficient %d  %0.3f"%(i+1,coef)

    print "\n\tIntercept %0.3f"%(model.intercept_)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\tMean squared error = %0.2f"%(mean_squared_error(true_y,predicted_y))

plot_residual函数用于绘制我们回归模型中的误差：

def plot_residual(y,predicted_y):
    """
    Plot residuals
    """
    plt.cla()
    plt.xlabel("Predicted Y")
    plt.ylabel("Residual")
    plt.title("Residual Plot")
    plt.figure(1)
    diff = y - predicted_y
    plt.plot(predicted_y,diff,'go')
    plt.show()

最后，我们将编写我们的main函数，用于调用之前的所有函数：

if __name__ == "__main__":

    x,y = get_data()

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    # Build the model
    model = build_model(x_train,y_train)
    predicted_y = model.predict(x_train)

    # Plot the residual
    plot_residual(y_train,predicted_y)
    # View model coeffiecients    
    view_model(model)

    print "\n Model Performance in Training set\n"
    model_worth(y_train,predicted_y)  

    # Apply the model on dev set
    predicted_y = model.predict(x_dev)
    print "\n Model Performance in Dev set\n"
    model_worth(y_dev,predicted_y)  

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(2)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_dev_poly   = poly_features.transform(x_dev)

    # Build model with polynomial features
    model_poly = build_model(x_train_poly,y_train)
    predicted_y = model_poly.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial features)\n"
    model_worth(y_train,predicted_y)  

    # Apply the model on dev set
    predicted_y = model_poly.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

    # Apply the model on Test set
    x_test_poly = poly_features.transform(x_test)
    predicted_y = model_poly.predict(x_test_poly)

    print "\n Model Performance in Test set  (Polynomial features)\n"
    model_worth(y_test,predicted_y)  

    predicted_y = model.predict(x_test)
    print "\n Model Performance in Test set  (Regular features)\n"
    model_worth(y_test,predicted_y)  

它是如何工作的…

让我们从主模块开始，跟着代码走。我们将使用 get_data 函数加载预测变量 x 和响应变量 y：

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

该函数调用了 scikit-learn 提供的便捷 load_boston() 函数，以便将波士顿房价数据集作为 NumPy 数组进行检索。

我们将使用 Scikit 库中的 train_test_split 函数，将数据划分为训练集和测试集。我们将保留数据集的 30% 用于测试：

x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)

我们将在下一行中提取开发集：

x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

在下一行，我们将通过调用 build_model 方法使用训练数据集来构建我们的模型。该模型创建一个 LinearRegression 类型的对象。LinearRegression 类封装了 SciPy 的最小二乘法：

    model = LinearRegression(normalize=True,fit_intercept=True)

让我们看看初始化该类时传递的参数。

fit_intercept 参数设置为 True。这告诉线性回归类进行数据中心化。通过数据中心化，我们将每个预测变量的均值设置为零。线性回归方法要求数据根据其均值进行中心化，以便更好地解释截距。除了按均值对每个特征进行中心化外，我们还将通过其标准差对每个特征进行归一化。我们将使用 normalize 参数并将其设置为 True 来实现这一点。有关如何按列进行归一化的更多信息，请参考第三章，缩放与数据标准化 的相关内容。通过 fit_intercept 参数，我们将指示算法包含一个截距，以适应响应变量中的任何常数偏移。最后，我们将通过调用 fit 函数并使用响应变量 y 和预测变量 x 来拟合模型。

注意

有关线性回归方法的更多信息，请参考 Trevor Hastie 等人所著的《统计学习元素》一书。

检查我们构建的模型是一种良好的做法，这样我们可以更好地理解模型，以便进一步改进或提升可解释性。

现在，让我们绘制残差图（预测的 y 与实际的 y 之间的差异），并将预测的 y 值作为散点图展示。我们将调用 plot_residual 方法来实现：

    # Plot the residual
    plot_residual(y_train,predicted_y)

让我们看看下面的图表：

我们可以使用这个散点图来验证数据集中的回归假设。我们没有看到任何模式，点均匀地分布在零残差值附近。

注意

有关使用残差图来验证线性回归假设的更多信息，请参考 Daniel. T. Larose 所著的《数据挖掘方法与模型》一书。

然后，我们将使用view_model方法检查我们的模型。在这个方法中，我们将打印截距和系数值。线性回归对象有两个属性，一个叫做coef_，它提供了系数的数组，另一个叫做intercept_，它提供截距值：

让我们来看一下coefficient 6，即房屋中可居住的房间数量。系数值的解释是：每增加一个房间，价格上升三倍。

最后，我们将通过调用model_worth函数来评估我们的模型好坏，该函数使用我们预测的响应值和实际响应值，二者均来自训练集和开发集。

这个函数打印出均方误差值，即实际值与预测值之间差值的平方的平均值：

我们在开发集中的值较低，这表明我们的模型表现良好。让我们看看是否可以改善均方误差。如果我们为模型提供更多特征会怎么样？让我们从现有的属性中创建一些特征。我们将使用 scikit-learn 中的PolynomialFeatures类来创建二阶多项式：

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(2)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_dev_poly   = poly_features.transform(x_dev)

我们将2作为参数传递给PolynomialFeatures，表示我们需要二阶多项式。如果类初始化为空，2也是默认值：

快速查看新x的形状，我们现在有 105 个属性，而原本只有 13 个。让我们使用新的多项式特征来构建模型，并检查模型的准确性：

    # Build model with polynomial features
    model_poly = build_model(x_train_poly,y_train)
    predicted_y = model_poly.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial           features)\n"
    model_worth(y_train,predicted_y)  

    # Apply the model on dev set
    predicted_y = model_poly.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

我们的模型已经很好地拟合了训练数据集。在开发集和训练集中，我们的多项式特征表现优于原始特征。

最后，让我们看看使用多项式特征的模型和使用常规特征的模型在测试集上的表现：

    # Apply the model on Test set
    x_test_poly = poly_features.transform(x_test)
    predicted_y = model_poly.predict(x_test_poly)

    print "\n Model Performance in Test set  (Polynomial features)\n"
    model_worth(y_test,predicted_y)  

predicted_y = model.predict(x_test)
    print "\n Model Performance in Test set  (Regular features)\n"
    model_worth(y_test,predicted_y)  

我们可以看到，在测试数据集上，使用多项式特征的表现优于我们原始特征的表现。

这就是你需要知道的如何在 Python 中进行线性回归。我们了解了线性回归的工作原理，以及如何构建模型来预测实数值。

还有更多…

在继续之前，我们将查看PolynomialFeatures类中的另一个参数设置，叫做interaction_only：

poly_features = PolynomialFeatures(interaction_only=True)

通过将interaction_only设置为true-with x1和x2 attributes-only，只会创建x1*x2属性。x1和x2的平方不会创建，假设阶数为二。

我们的测试集结果与开发集的结果不如预期，尤其是在普通特征和多项式特征上。这是线性回归的一个已知问题。线性回归对于处理方差的能力较弱。我们面临的问题是高方差和低偏差。随着模型复杂度的增加，也就是呈现给模型的特征数量增多，模型往往能够很好地拟合训练数据——因此偏差较低——但开始在测试数据上出现下降的结果。对此问题，有几种技术可以应对。

我们来看看一个叫做递归特征选择的方法。所需特征的数量作为参数传递给该方法。它递归地筛选特征。在第 i 次运行中，会对数据拟合一个线性模型，并根据系数的值来筛选特征；那些权重较低的特征会被剔除。如此一来，迭代就继续进行，直到我们得到所需数量的特征时，迭代才会停止。接下来，让我们看看一个代码示例：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
from itertools import combinations
from sklearn.feature_selection import RFE

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

def build_model(x,y,no_features):
    """
    Build a linear regression model
    """
    model = LinearRegression(normalize=True,fit_intercept=True)
    rfe_model = RFE(estimator=model,n_features_to_select=no_features)
    rfe_model.fit(x,y)
    return rfe_model    

def view_model(model):
    """
    Look at model coeffiecients
    """
    print "\n Model coeffiecients"
    print "======================\n"
    for i,coef in enumerate(model.coef_):
        print "\tCoefficient %d  %0.3f"%(i+1,coef)

    print "\n\tIntercept %0.3f"%(model.intercept_)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\tMean squared error = %0.2f"%(mean_squared_error(true_y,predicted_y))
    return mean_squared_error(true_y,predicted_y)

def plot_residual(y,predicted_y):
    """
    Plot residuals
    """
    plt.cla()
    plt.xlabel("Predicted Y")
    plt.ylabel("Residual")
    plt.title("Residual Plot")
    plt.figure(1)
    diff = y - predicted_y
    plt.plot(predicted_y,diff,'go')
    plt.show()

def subset_selection(x,y):
    """
    subset selection method
    """
    # declare variables to track
    # the model and attributes which produces
    # lowest mean square error
    choosen_subset = None
    low_mse = 1e100
    choosen_model = None
    # k value ranges from 1 to the number of 
    # attributes in x
    for k in range(1,x.shape[1]+1):
        print "k= %d "%(k)
        # evaluate all attribute combinations
        # of size k+1
        subsets = combinations(range(0,x.shape[1]),k+1)
        for subset in subsets:
            x_subset = x[:,subset]
            model = build_model(x_subset,y)
            predicted_y = model.predict(x_subset)
            current_mse = mean_squared_error(y,predicted_y)
            if current_mse < low_mse:
                low_mse = current_mse
                choosen_subset = subset
                choosen_model = model

    return choosen_model, choosen_subset,low_mse    

if __name__ == "__main__":

    x,y = get_data()

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(interaction_only=True)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_dev_poly   = poly_features.transform(x_dev)

    #choosen_model,choosen_subset,low_mse = subset_selection(x_train_poly,y_train)    
    choosen_model = build_model(x_train_poly,y_train,20)
    #print choosen_subse
    predicted_y = choosen_model.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial features)\n"
    mse  = model_worth(y_train,predicted_y)  

    # Apply the model on dev set
    predicted_y = choosen_model.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

    # Apply the model on Test set
    x_test_poly = poly_features.transform(x_test)
    predicted_y = choosen_model.predict(x_test_poly)

    print "\n Model Performance in Test set  (Polynomial features)\n"
    model_worth(y_test,predicted_y)  

这段代码与之前的线性回归代码非常相似，唯一不同的是build_model方法：

def build_model(x,y,no_features):
    """
    Build a linear regression model
    """
    model = LinearRegression(normalize=True,fit_intercept=True)
    rfe_model = RFE(estimator=model,n_features_to_select=no_features)
    rfe_model.fit(x,y)
    return rfe_model    

除了预测变量x和响应变量y，build_model方法还接受一个参数，即要保留的特征数量no_features。在这个例子中，我们传递了一个值 20，要求递归特征消除方法只保留 20 个重要特征。如你所见，我们首先创建了一个线性回归对象。这个对象被传递给RFE类。RFE 代表递归特征消除，这是 scikit-learn 提供的一个类，用于实现递归特征消除。现在，我们来评估一下我们的模型在训练集、开发集和测试集上的表现：

测试数据集的均方误差为 13.20，几乎是之前的一半。因此，我们能够有效地使用递归特征消除方法进行特征选择，从而提升我们的模型。

另请参见

• 第三章中的数据标准化方法，数据分析 – 探索与整理

• 第三章中的数据标准化方法，数据分析 – 探索与整理

• 第六章中的为模型构建准备数据方法，机器学习 I

使用 L2 收缩学习回归 – 岭回归

让我们将之前讨论的回归技术扩展，以包括正则化。在训练线性回归模型时，一些系数可能会取非常高的值，导致模型的不稳定。正则化或收缩是一种控制系数权重的方法，目的是使它们不会取非常大的值。让我们再次看看线性回归的成本函数，以理解回归固有的问题，以及我们通过控制系数权重所指的内容：

线性回归试图找到系数w0…wm，使得它最小化上述方程。线性回归存在一些问题。

如果数据集包含许多相关的预测变量，数据的非常小的变化可能导致模型不稳定。此外，我们还会面临解释模型结果的问题。例如，如果我们有两个负相关的变量，这些变量将对响应变量产生相反的影响。我们可以手动查看这些相关性，移除其中一个有影响的变量，然后继续构建模型。然而，如果我们能够自动处理这些场景，那将更有帮助。

我们在之前的例子中介绍了一种方法，叫做递归特征消除，用来保留最具信息性的属性并丢弃其余的属性。然而，在这种方法中，我们要么保留一个变量，要么不保留它；我们的决策是二元的。在这一部分，我们将看到一种方法，能够以控制变量权重的方式，强烈惩罚不必要的变量，使它们的权重降到极低。

我们将修改线性回归的成本函数，以包含系数。如你所知，成本函数的值应该在最佳模型时达到最小。通过将系数包含在成本函数中，我们可以对取非常高值的系数进行严重惩罚。一般来说，这些技术被称为收缩方法，因为它们试图收缩系数的值。在这个例子中，我们将看到 L2 收缩，通常称为岭回归。让我们看看岭回归的成本函数：

如你所见，系数的平方和被加入到成本函数中。这样，当优化程序试图最小化上述函数时，它必须大幅减少系数的值才能达到目标。α参数决定了收缩的程度。α值越大，收缩越大。系数值会被减小至零。

在了解了这些数学背景后，让我们进入实际操作，看看岭回归的应用。

准备工作

再次使用波士顿数据集来演示岭回归。波士顿数据集有 13 个属性和 506 个实例。目标变量是一个实数，房屋的中位数价格在几千美元之间。有关波士顿数据集的更多信息，请参考以下 UCI 链接：

archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names

我们打算生成二次的多项式特征，并仅考虑交互效应。在本教程的最后，我们将看到系数如何受到惩罚。

如何操作……

我们将首先加载所有必要的库。接下来，我们将定义第一个函数get_data()。在这个函数中，我们将读取波士顿数据集，并将其返回为预测变量x和响应变量y：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    x    = x - np.mean(x,axis=0)

    return x,y    

在我们的下一个build_model函数中，我们将使用给定的数据构建岭回归模型。以下两个函数，view_model和model_worth，用于检查我们构建的模型：

def build_model(x,y):
    """
    Build a Ridge regression model
    """
    model = Ridge(normalize=True,alpha=0.015)
    model.fit(x,y)
    # Track the scores- Mean squared residual for plot
    return model    

def view_model(model):
    """
    Look at model coeffiecients
    """
    print "\n Model coeffiecients"
    print "======================\n"
    for i,coef in enumerate(model.coef_):
        print "\tCoefficient %d  %0.3f"%(i+1,coef)

    print "\n\tIntercept %0.3f"%(model.intercept_)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\tMean squared error = %0.2f"%(mean_squared_error(true_y,predicted_y))
    return mean_squared_error(true_y,predicted_y)

最后，我们将编写我们的main函数，该函数用于调用所有前面的函数：

if __name__ == "__main__":

    x,y = get_data()

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(interaction_only=True)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_dev_poly   = poly_features.transform(x_dev)
    x_test_poly = poly_features.transform(x_test)

    #choosen_model,choosen_subset,low_mse = subset_selection(x_train_poly,y_train)    
    choosen_model = build_model(x_train_poly,y_train)

    predicted_y = choosen_model.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial features)\n"
    mse  = model_worth(y_train,predicted_y)  
    view_model(choosen_model)

    # Apply the model on dev set
    predicted_y = choosen_model.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

    # Apply the model on Test set
    predicted_y = choosen_model.predict(x_test_poly)

    print "\n Model Performance in Test set  (Polynomial features)\n"
    model_worth(y_test,predicted_y)  

它是如何工作的……

让我们从主模块开始，并跟随代码。我们使用get_data函数加载了预测变量x和响应变量y。该函数调用了 scikit-learn 的便捷函数load_boston()，将波士顿房价数据集作为 NumPy 数组获取。

接下来，我们将使用 scikit-learn 库中的train_test_split函数将数据集划分为训练集和测试集。我们将保留 30%的数据集用于测试。在下一行，我们将从中提取开发集。

然后，我们将构建多项式特征：

poly_features = PolynomialFeatures(interaction_only=True)
poly_features.fit(x_train)

如你所见，我们将interaction_only设置为true。通过将interaction_only设置为true，对于x1和x2属性，仅创建x1*x2属性，而不创建x1和x2的平方，假设度数为 2。默认的度数为 2：

x_train_poly = poly_features.transform(x_train)
x_dev_poly = poly_features.transform(x_dev)
x_test_poly = poly_features.transform(x_test)

使用transform函数，我们将转换我们的训练集、开发集和测试集，以包含多项式特征。

在下一行，我们将通过调用build_model方法使用训练数据集构建我们的岭回归模型：

model = Ridge(normalize=True,alpha=0.015)
model.fit(x,y)

数据集中的属性通过其均值进行中心化，并使用标准差进行标准化，方法是使用normalize参数并将其设置为true。Alpha控制收缩的程度，其值设置为0.015。我们不是凭空得出这个数字的，而是通过多次运行模型得出的。稍后在本章中，我们将看到如何通过经验得出此参数的正确值。我们还将使用fit_intercept参数来拟合该模型的截距。然而，默认情况下，fit_intercept参数设置为true，因此我们不会显式指定它。

现在让我们看看模型在训练集中的表现。我们将调用 model_worth 方法来获取均方误差。此方法需要预测的响应变量和实际的响应变量，以返回均方误差：

predicted_y = choosen_model.predict(x_train_poly)
print "\n Model Performance in Training set (Polynomial features)\n"
mse = model_worth(y_train,predicted_y) 

我们的输出如下所示：

在将模型应用于测试集之前，让我们先看看系数的权重。我们将调用一个名为 view_model 的函数来查看系数的权重：

view_model(choosen_model)

我们没有展示所有的系数，共有 92 个。然而，从一些系数来看，收缩效应应该是显而易见的。例如，系数 1 几乎为 0（请记住，它是一个非常小的值，这里仅显示了前三位小数）。

让我们继续看看模型在开发集中的表现：

predicted_y = choosen_model.predict(x_dev_poly)
print "\n Model Performance in Dev set (Polynomial features)\n"
model_worth(y_dev,predicted_y) 

不错，我们已经达到了比训练误差更低的均方误差。最后，让我们来看一下我们模型在测试集上的表现：

与前面配方中的线性回归模型相比，我们在测试集上的表现更好。

还有更多…

我们之前提到过，线性回归模型对数据集中的任何小变化都非常敏感。让我们看一个小例子，来演示这一点：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    x    = x - np.mean(x,axis=0)

    return x,y    

在这段代码中，我们将使用 build_model 函数对原始数据进行线性回归和岭回归模型的拟合：

lin_model,ridg_model = build_model(x,y)

我们将在原始数据中引入少量噪声，具体如下：

# Add some noise to the dataset
noise = np.random.normal(0,1,(x.shape))
x = x + noise

我们将再次在噪声数据集上拟合模型。最后，我们将比较系数的权重：

在添加了少量噪声后，当我们尝试使用线性回归拟合模型时，分配的权重与前一个模型分配的权重非常不同。现在，让我们看看岭回归的表现：

权重在第一个和第二个模型之间没有发生剧烈变化。希望这可以展示岭回归在噪声数据条件下的稳定性。

选择合适的 alpha 值总是很棘手。一个粗暴的方法是通过多个值来运行，并跟踪系数的路径。通过路径，选择一个系数变化不剧烈的 alpha 值。我们将使用 coeff_path 函数绘制系数的权重。

我们来看看 coeff_path 函数。它首先生成一组 alpha 值：

alpha_range = np.linspace(10,100.2,300)

在这种情况下，我们生成了从 10 到 100 之间均匀分布的 300 个数字。对于这些 alpha 值中的每一个，我们将构建一个模型并保存其系数：

for alpha in alpha_range:
    model = Ridge(normalize=True,alpha=alpha)
    model.fit(x,y)
    coeffs.append(model.coef_)

最后，我们将绘制这些系数权重与 alpha 值的关系图：

如您所见，值在 alpha 值为 100 附近趋于稳定。您可以进一步缩放到接近 100 的范围，并寻找一个理想的值。

另见

• 在第七章中的使用回归预测实数值方法，机器学习 II

• 在第三章中的数据缩放方法，数据分析 – 探索与整理

• 在第三章中的标准化数据方法，数据分析 – 探索与整理

• 在第六章中的准备数据以构建模型方法，机器学习 I

使用 L1 收缩学习回归 – LASSO

最小绝对收缩与选择算子（LASSO）是另一种常用的回归问题收缩方法。与岭回归相比，LASSO 能得到稀疏解。所谓稀疏解，指的是大部分系数被缩减为零。在 LASSO 中，很多系数被设为零。对于相关变量，LASSO 只选择其中一个，而岭回归则为两个变量的系数分配相等的权重。因此，LASSO 的这一特性可以用来进行变量选择。在这个方法中，让我们看看如何利用 LASSO 进行变量选择。

让我们看一下 LASSO 回归的代价函数。如果你已经跟随前两个方法，你应该能很快识别出区别：

系数会受到系数绝对值之和的惩罚。再次强调，alpha 控制惩罚的程度。我们来尝试理解为什么 L1 收缩会导致稀疏解的直观原因。

我们可以将前面的方程重写为一个无约束的代价函数和一个约束，如下所示：

最小化：

受约束条件的影响：

记住这个方程，我们在系数空间中绘制w0和w1的代价函数值：

蓝色线条表示不同w0和w1值下代价函数（无约束）的等高线。绿色区域表示由 eta 值决定的约束形状。两个区域交汇的优化值是w0设为 0 时的情况。我们展示了一个二维空间，在这个空间中，通过将w0设为 0，我们的解变得稀疏。在多维空间中，我们将有一个菱形区域，LASSO 通过将许多系数缩减为零，给出一个稀疏解。

准备开始

我们将再次使用波士顿数据集来演示 LASSO 回归。波士顿数据集有 13 个属性和 506 个实例。目标变量是一个实数，房屋的中位数价值在千元范围内。

有关更多波士顿数据集的信息，请参考以下 UCI 链接：

archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names。

我们将看到如何使用 LASSO 进行变量选择。

如何操作……

我们将首先加载所有必要的库。接下来我们定义我们的第一个函数get_data()。在此函数中，我们将读取波士顿数据集，并将其作为预测变量x和响应变量y返回：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import train_test_split
from sklearn.linear_model import Lasso, LinearRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from sklearn.preprocessing import PolynomialFeatures
import numpy as np

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

在我们接下来的build_model函数中，我们将使用给定的数据构建 LASSO 回归模型。接下来的两个函数，view_model和model_worth，用于检查我们构建的模型：

def build_models(x,y):
    """
    Build a Lasso regression model
    """
    # Alpha values uniformly
    # spaced between 0.01 and 0.02
    alpha_range = np.linspace(0,0.5,200)
    model = Lasso(normalize=True)
    coeffiecients = []
    # Fit a model for each alpha value
    for alpha in alpha_range:
        model.set_params(alpha=alpha)
        model.fit(x,y)
        # Track the coeffiecients for plot
        coeffiecients.append(model.coef_)
    # Plot coeffients weight decay vs alpha value
    # Plot model RMSE vs alpha value
    coeff_path(alpha_range,coeffiecients)
    # View coeffiecient value
    #view_model(model)

def view_model(model):
    """
    Look at model coeffiecients
    """
    print "\n Model coeffiecients"
    print "======================\n"
    for i,coef in enumerate(model.coef_):
        print "\tCoefficient %d  %0.3f"%(i+1,coef)

    print "\n\tIntercept %0.3f"%(model.intercept_)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\t Mean squared error = %0.2f\n"%(mean_squared_error(true_y,predicted_y))

我们将定义两个函数，coeff_path和get_coeff，来检查我们的模型系数。coeff_path函数由build_model函数调用，用于绘制不同 alpha 值下的系数权重。get_coeff函数则由主函数调用：

def coeff_path(alpha_range,coeffiecients):
    """
    Plot residuals
    """
    plt.close('all')
    plt.cla()

    plt.figure(1)
    plt.xlabel("Alpha Values")
    plt.ylabel("Coeffiecient Weight")
    plt.title("Coeffiecient weights for different alpha values")
    plt.plot(alpha_range,coeffiecients)
    plt.axis('tight')

    plt.show()

def get_coeff(x,y,alpha):
    model = Lasso(normalize=True,alpha=alpha)
    model.fit(x,y)
    coefs = model.coef_
    indices = [i for i,coef in enumerate(coefs) if abs(coef) > 0.0]
    return indices

最后，我们将编写我们的main函数，用于调用之前所有的函数：

if __name__ == "__main__":

    x,y = get_data()
    # Build multiple models for different alpha values
    # and plot them    
    build_models(x,y)
    print "\nPredicting using all the variables"
    full_model = LinearRegression(normalize=True)
    full_model.fit(x,y)
    predicted_y = full_model.predict(x)
    model_worth(y,predicted_y)    

    print "\nModels at different alpha values\n"
    alpa_values = [0.22,0.08,0.01]
    for alpha in alpa_values:

        indices = get_coeff(x,y,alpha)   
        print "\t alpah =%0.2f Number of variables selected = %d "%(alpha,len(indices))
        print "\t attributes include ", indices
        x_new = x[:,indices]
        model = LinearRegression(normalize=True)
        model.fit(x_new,y)
        predicted_y = model.predict(x_new)
        model_worth(y,predicted_y)

它是如何工作的……

让我们从主模块开始，跟随代码进行。我们将使用get_data函数加载预测变量x和响应变量y。该函数调用了 scikit-learn 提供的便捷load_boston()函数，从而将波士顿房价数据集作为 NumPy 数组载入。

我们将继续调用build_models。在build_models中，我们将为不同的alpha值构建多个模型：

alpha_range = np.linspace(0,0.5,200)
model = Lasso(normalize=True)
coeffiecients = []
# Fit a model for each alpha value
for alpha in alpha_range:
model.set_params(alpha=alpha)
model.fit(x,y)
# Track the coeffiecients for plot
coeffiecients.append(model.coef_)

正如你所看到的，在 for 循环中，我们还将不同 alpha 值的系数值存储在一个列表中。

让我们通过调用coeff_path函数来绘制不同 alpha 值下的系数值：

plt.close('all')
plt.cla()

plt.figure(1)
plt.xlabel("Alpha Values")
plt.ylabel("Coeffiecient Weight")
plt.title("Coeffiecient weights for different alpha values")
plt.plot(alpha_range,coeffiecients)
plt.axis('tight')
plt.show()

在x轴上，你可以看到我们有 alpha 值，而在y轴上，我们将为给定的 alpha 值绘制所有系数。让我们看看输出的图表：

不同颜色的线条代表不同的系数值。正如你所看到的，随着 alpha 值的增大，系数权重趋向于零。从这个图表中，我们可以选择 alpha 的值。

作为参考，我们先拟合一个简单的线性回归模型：

print "\nPredicting using all the variables"
full_model = LinearRegression(normalize=True)
full_model.fit(x,y)
predicted_y = full_model.predict(x)
model_worth(y,predicted_y) 

让我们看看当我们尝试使用新构建的模型进行预测时的均方误差：

让我们继续根据 LASSO 来选择系数：

print "\nModels at different alpha values\n"
alpa_values = [0.22,0.08,0.01]
for alpha in alpa_values:
indices = get_coeff(x,y,alpha) 

根据之前的图表，我们选择了0.22、0.08和0.01作为 alpha 值。在循环中，我们将调用get_coeff方法。该方法使用给定的 alpha 值拟合 LASSO 模型，并仅返回非零系数的索引：

model = Lasso(normalize=True,alpha=alpha)
model.fit(x,y)
coefs = model.coef_

indices = [i for i,coef in enumerate(coefs) if abs(coef) > 0.0]

本质上，我们只选择那些系数值非零的属性——特征选择。让我们回到我们的for循环，在那里我们将使用减少后的系数拟合线性回归模型：

print "\t alpah =%0.2f Number of variables selected = %d "%(alpha,len(indices))
print "\t attributes include ", indices
x_new = x[:,indices]
model = LinearRegression(normalize=True)
model.fit(x_new,y)
predicted_y = model.predict(x_new)
model_worth(y,predicted_y)

我们想要了解的是，如果我们使用减少后的属性集来预测模型，与使用整个数据集最初构建的模型相比，模型的效果如何：

看看我们的第一次尝试，alpha 值为0.22时。只有两个系数的非零值，分别为5和12。均方误差为30.51，仅比使用所有变量拟合的模型多了9。

类似地，对于0.08的 alpha 值，存在三个非零系数。我们可以看到均方误差有所改善。最后，对于0.01的 alpha 值，13 个属性中有 9 个被选择，且均方误差非常接近使用所有属性构建的模型。

如你所见，我们并没有使用所有属性来拟合模型。我们能够使用 LASSO 自动选择一部分属性子集。因此，我们已经看到了 LASSO 如何用于变量选择。

还有更多内容……

通过保留最重要的变量，LASSO 避免了过拟合。然而，如你所见，均方误差值并不理想。我们可以看到，LASSO 导致了预测能力的损失。

如前所述，对于相关变量，LASSO 只选择其中一个，而岭回归则对两个变量的系数赋予相等的权重。因此，与 LASSO 相比，岭回归具有更高的预测能力。然而，LASSO 能够进行变量选择，而岭回归则不能。

注意

请参考Trevor Hastie 等人的《稀疏统计学习：Lasso 与泛化》一书，了解更多关于 LASSO 和岭回归方法的信息。

另见：

• 在第三章中，数据标准化的内容，数据分析 – 探索与清洗

• 在第三章中，数据标准化的内容，数据分析 – 探索与清洗

• 在第六章中，模型构建数据准备的内容，机器学习 I

• 在第七章中，L2 收缩回归 – 岭回归的内容，机器学习 II

使用交叉验证迭代器与 L1 和 L2 收缩

在上一章中，我们介绍了将数据划分为训练集和测试集的方法。在接下来的步骤中，我们再次对测试数据集进行划分，以得到开发数据集。其目的是将测试集从模型构建周期中分离出来。然而，由于我们需要不断改进模型，因此我们使用开发集来测试每次迭代中的模型准确性。虽然这是一个好的方法，但如果数据集不是很大，实施起来会比较困难。我们希望尽可能多地提供数据来训练模型，但仍然需要保留部分数据用于评估和最终测试。在许多实际场景中，获取一个非常大的数据集是非常罕见的。

在这个示例中，我们将看到一种叫做交叉验证的方法，帮助我们解决这个问题。这个方法通常被称为 k 折交叉验证。训练集被分成 k 个折叠。模型在 K-1（K 减 1）个折叠上进行训练，剩下的折叠用来测试。这样，我们就不需要单独的开发数据集了。

让我们看看 scikit-learn 库提供的一些迭代器，以便有效地执行 k 折交叉验证。掌握了交叉验证的知识后，我们将进一步了解如何利用交叉验证选择收缩方法中的 alpha 值。

准备工作

我们将使用鸢尾花数据集来演示各种交叉验证迭代器的概念。接着，我们会回到波士顿住房数据集，演示如何利用交叉验证成功找到收缩方法中的理想 alpha 值。

如何实现…

让我们看看如何使用交叉验证迭代器：

from sklearn.datasets import load_iris
from sklearn.cross_validation import KFold,StratifiedKFold

def get_data():
    data = load_iris()
    x = data['data']
    y = data['target']
    return x,y

def class_distribution(y):
        class_dist = {}
        total = 0
        for entry in y:
            try:
                class_dist[entry]+=1
            except KeyError:
                class_dist[entry]=1
            total+=1

        for k,v in class_dist.items():
            print "\tclass %d percentage =%0.2f"%(k,v/(1.0*total))

if __name__ == "__main__":
    x,y = get_data()
    # K Fold
    # 3 folds
    kfolds = KFold(n=y.shape[0],n_folds=3)
    fold_count  =1
    print
    for train,test in kfolds:
        print "Fold %d x train shape"%(fold_count),x[train].shape,\
        " x test shape",x[test].shape
        fold_count==1
    print
    #Stratified KFold
    skfolds = StratifiedKFold(y,n_folds=3)
    fold_count  =1
    for train,test in skfolds:
        print "\nFold %d x train shape"%(fold_count),x[train].shape,\
        " x test shape",x[test].shape
        y_train = y[train]
        y_test  = y[test]
        print "Train Class Distribution"
        class_distribution(y_train)
        print "Test Class Distribution"
        class_distribution(y_test)

        fold_count+=1

    print

在我们的主函数中，我们将调用get_data函数来加载鸢尾花数据集。然后我们将演示一个简单的 k 折和分层 k 折交叉验证。

通过掌握 k 折交叉验证的知识，让我们写一个食谱，利用这些新获得的知识来提升岭回归：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import KFold,train_test_split
from sklearn.linear_model import Ridge
from sklearn.grid_search import GridSearchCV
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures
import numpy as np

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

我们将首先加载所有必要的库。接下来，我们将定义我们的第一个函数get_data()。在这个函数中，我们将读取波士顿数据集，并将其返回为预测变量x和响应变量y。

在我们的下一个build_model函数中，我们将使用给定的数据构建岭回归模型，并利用 k 折交叉验证。

以下两个函数，view_model和model_worth，用于 introspect（内省）我们构建的模型。

最后，我们将编写display_param_results函数，以查看每个折叠中的模型误差：

def build_model(x,y):
    """
    Build a Ridge regression model
    """
    kfold = KFold(y.shape[0],5)
    model = Ridge(normalize=True)

    alpha_range = np.linspace(0.0015,0.0017,30)
    grid_param = {"alpha":alpha_range}
    grid = GridSearchCV(estimator=model,param_grid=grid_param,cv=kfold,scoring='mean_squared_error')
    grid.fit(x,y)
    display_param_results(grid.grid_scores_)
    print grid.best_params_
    # Track the scores- Mean squared residual for plot
    return grid.best_estimator_

def view_model(model):
    """
    Look at model coeffiecients
    """
    #print "\n Estimated Alpha = %0.3f"%model.alpha_
    print "\n Model coeffiecients"
    print "======================\n"
    for i,coef in enumerate(model.coef_):
        print "\tCoefficient %d  %0.3f"%(i+1,coef)

    print "\n\tIntercept %0.3f"%(model.intercept_)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\tMean squared error = %0.2f"%(mean_squared_error(true_y,predicted_y))
    return mean_squared_error(true_y,predicted_y)

def display_param_results(param_results):
    fold = 1
    for param_result in param_results:
        print "Fold %d Mean squared error %0.2f"%(fold,abs(param_result[1])),param_result[0]
        fold+=1

最后，我们将编写我们的main函数，用来调用所有前面的函数：

if __name__ == "__main__":

    x,y = get_data()

    # Divide the data into Train and test    
    x_train,x_test,y_train,y_test = train_test_split(x,y,test_size = 0.3,random_state=9)

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(interaction_only=True)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_test_poly  = poly_features.transform(x_test)

    choosen_model = build_model(x_train_poly,y_train)
    predicted_y = choosen_model.predict(x_train_poly)
    model_worth(y_train,predicted_y)

    view_model(choosen_model)

    predicted_y = choosen_model.predict(x_test_poly)
    model_worth(y_test,predicted_y)

工作原理…

让我们从主方法开始。我们将从KFold类开始。这个迭代器类是通过实例化数据集中实例的数量和我们所需的折数来创建的：

kfolds = KFold(n=y.shape[0],n_folds=3)

现在，我们可以按以下方式遍历这些折叠：

fold_count =1
print
for train,test in kfolds:
print "Fold %d x train shape"%(fold_count),x[train].shape,\
" x test shape",x[test].shape
fold_count==1

让我们看看打印语句的输出：

我们可以看到数据被分成三部分，每部分有 100 个实例用于训练，50 个实例用于测试。

接下来，我们将转到StratifiedKFold。回想一下我们在上一章中讨论的关于训练集和测试集中的类分布均匀的问题。StratifiedKFold在三个折叠中实现了均匀的类分布。

它是这样调用的：

skfolds = StratifiedKFold(y,n_folds=3)

由于它需要知道数据集中类标签的分布，因此这个迭代器对象将响应变量y作为其参数之一。另一个参数是请求的折叠数。

让我们打印出这三折数据集中训练集和测试集的形状，并查看它们的类分布。我们将使用class_distribution函数打印每一折中的类分布：

fold_count =1
for train,test in skfolds:
print "\nFold %d x train shape"%(fold_count),x[train].shape,\
" x test shape",x[test].shape
y_train = y[train]
y_test = y[test]
print "Train Class Distribution"
class_distribution(y_train)
print "Test Class Distribution"
class_distribution(y_test)

fold_count+=1

你可以看到类的分布是均匀的。

假设你构建了一个五折交叉验证的数据集，拟合了五个不同的模型，并得到了五个不同的准确度分数。现在，你可以取这些分数的平均值来评估模型的好坏。如果你不满意，你可以继续调整参数，重新构建模型，再在五折数据上运行并查看平均准确度分数。这样，你可以通过仅使用训练数据集找到合适的参数值，不断改进模型。

掌握了这些知识后，我们来重新审视我们之前的岭回归问题。

让我们从main模块开始，跟着代码走。我们将使用get_data函数加载预测变量x和响应变量y。该函数调用 scikit-learn 便捷的load_boston()函数，将波士顿房价数据集作为 NumPy 数组提取出来。

接下来，我们将使用 scikit-learn 库中的train_test_split函数将数据分割为训练集和测试集。我们将保留 30%的数据集用于测试。

然后我们将构建多项式特征：

poly_features = PolynomialFeatures(interaction_only=True)
poly_features.fit(x_train)

如你所见，我们将interaction_only设置为true。通过将interaction_only设置为true——并配合x1和x2属性——只会创建x1*x2属性。x1和x2的平方不会被创建，假设度数为二。默认度数为二：

x_train_poly = poly_features.transform(x_train)
x_test_poly = poly_features.transform(x_test)

使用transform函数，我们将转换训练集和测试集，以包含多项式特征。我们来调用build_model函数。在build_model函数中，我们首先注意到的是 k 折交叉验证的声明。我们将在此应用交叉验证的知识，并创建一个五折数据集：

kfold = KFold(y.shape[0],5)

然后我们将创建我们的岭回归对象：

model = Ridge(normalize=True)

现在让我们看看如何利用 k 折交叉验证来找出岭回归的理想 alpha 值。在接下来的行中，我们将使用GridSearchCV创建一个对象：

grid = GridSearchCV(estimator=model,param_grid=grid_param,cv=kfold,scoring='mean_squared_error')

GridSearchCV是 scikit-learn 中的一个便捷函数，帮助我们使用一系列参数训练模型。在这种情况下，我们希望找到理想的 alpha 值，因此我们想用不同的 alpha 值训练我们的模型。让我们看看传递给GridSearchCV的参数：

估算器：这是应该使用给定参数和数据运行的模型类型。在我们的例子中，我们希望运行岭回归。因此，我们将创建一个岭回归对象，并将其传递给GridSearchCV。

参数网格：这是我们希望用来评估模型的参数字典。让我们详细处理这个问题。我们将首先声明要用于构建模型的 alpha 值范围：

alpha_range = np.linspace(0.0015,0.0017,30)

这给我们一个 NumPy 数组，包含 30 个均匀间隔的元素，起始于 0.0015，结束于 0.0017。我们希望为每个这些值构建一个模型。我们将创建一个名为grid_param的字典对象，并在 alpha 键下添加生成的 NumPy 数组作为 alpha 值：

grid_param = {"alpha":alpha_range}

我们将把这个字典作为参数传递给GridSearchCV。看一下条目param_grid=grid_param。

Cv：这定义了我们感兴趣的交叉验证类型。我们将传递之前创建的 k 折（五折）迭代器作为 cv 参数。

最后，我们需要定义一个评分函数。在我们的案例中，我们关注的是求得平方误差。这是我们用来评估模型的指标。

因此，内部的GridSearchCV将为我们的每个参数值构建五个模型，并在排除的折中测试时返回平均分数。在我们的案例中，我们有五个测试折，因此会返回这五个测试折中分数的平均值。

解释完这些后，我们将开始拟合我们的模型，也就是说，启动网格搜索活动。

最后，我们想查看不同参数设置下的输出。我们将使用display_param_results函数来显示不同折中的平均均方误差：

输出的每一行显示了参数 alpha 值和来自测试折的平均均方误差。我们可以看到，当我们深入到 0.0016 的范围时，均方误差在增加。因此，我们决定停在 0.0015。我们可以查询网格对象以获取最佳参数和估计器：

print grid.best_params_
return grid.best_estimator_

这不是我们测试的第一组 alpha 值。我们最初的 alpha 值如下：

alpha_range = np.linspace(0.01,1.0,30)

以下是我们的输出：

当我们的 alpha 值超过 0.01 时，均方误差急剧上升。因此，我们又给出了一个新的范围：

alpha_range = np.linspace(0.001,0.1,30)

我们的输出如下：

这样，我们通过迭代的方式得到了一个范围，从 0.0015 开始，到 0.0017 结束。

然后，我们将从我们的网格搜索中获取最佳估计器，并将其应用于我们的训练和测试数据：

choosen_model = build_model(x_train_poly,y_train)
predicted_y = choosen_model.predict(x_train_poly)
model_worth(y_train,predicted_y)

我们的model_worth函数打印出我们训练数据集中的均方误差值：

让我们查看我们的系数权重：

我们没有显示所有的系数，但当你运行代码时，可以查看所有的值。

最后，让我们将模型应用到我们的测试数据集：

因此，我们通过交叉验证和网格搜索成功地得出了岭回归的 alpha 值。与岭回归配方中的值相比，我们的模型实现了较低的均方误差。

还有更多内容……

scikit-learn 提供了其他交叉验证迭代器。在这种情况下，特别值得关注的是留一法迭代器。你可以在scikit-learn.org/stable/modules/cross_validation.html#leave-one-out-loo上阅读更多关于它的信息。

在这种方法中，给定折数，它会留出一条记录进行测试，并将其余的记录用于训练。例如，如果你的输入数据有 100 个实例，并且我们要求五折交叉验证，那么每一折中我们将有 99 个实例用于训练，一个实例用于测试。

在我们之前使用的网格搜索方法中，如果没有为交叉验证（cv）参数提供自定义迭代器，它将默认使用留一法交叉验证方法：

grid = GridSearchCV(estimator=model,param_grid=grid_param,cv=None,scoring='mean_squared_error')

另见

• 第三章中的数据缩放方法，数据分析 – 探索与整理

• 第三章中的标准化数据方法，数据分析 – 探索与整理

• 第六章中的为模型构建准备数据方法，机器学习 I

• 第七章中的L2 缩减回归 – 岭回归方法，机器学习 II

• 第七章中的L2 缩减回归 – Lasso方法，机器学习 II

第八章 集成方法

在本章中，我们将讨论以下内容：

• 理解集成方法——袋装法（Bagging）

• 理解集成方法——提升法（Boosting），AdaBoost

• 理解集成方法——梯度提升（Gradient Boosting）

介绍

在本章中，我们将讨论涉及集成方法的内容。当我们在现实生活中面临不确定性决策时，通常会向多个朋友寻求意见。我们根据从这些朋友那里获得的集体知识做出决策。集成方法在机器学习中的概念类似。前几章中，我们为我们的数据集构建了单一模型，并使用该模型对未见过的测试数据进行预测。如果我们在这些数据上构建多个模型，并根据所有这些模型的预测结果做出最终预测，会怎么样呢？这就是集成方法的核心思想。使用集成方法解决某个问题时，我们会构建多个模型，并利用它们所有的预测结果来对未见过的数据集做出最终预测。对于回归问题，最终输出可能是所有模型预测值的平均值。在分类问题中，则通过多数投票来决定输出类别。

基本思想是拥有多个模型，每个模型在训练数据集上产生略有不同的结果。有些模型比其他模型更好地学习数据的某些方面。其信念是，所有这些模型的最终输出应该优于任何单一模型的输出。

如前所述，集成方法的核心是将多个模型结合在一起。这些模型可以是相同类型的，也可以是不同类型的。例如，我们可以将神经网络模型的输出与贝叶斯模型的输出结合在一起。本章中我们将限制讨论使用相同类型的模型进行集成。通过像袋装法（Bagging）和提升法（Boosting）这样的技术，数据科学社区广泛使用相同类型模型的集成方法。

自助聚合（Bootstrap aggregation），通常称为袋装法（Bagging），是一种优雅的技术，用于生成大量模型并将它们的输出结合起来以做出最终预测。袋装法中的每个模型只使用部分训练数据。袋装法的核心思想是减少数据的过拟合。如前所述，我们希望每个模型与其他模型略有不同。因此，我们对每个模型的训练数据进行有放回的抽样，从而引入了变异性。引入模型变异的另一种方式是对属性进行抽样。我们不会将所有属性提供给模型，而是让不同的模型获得不同的属性集。袋装法可以轻松并行化。基于现有的并行处理框架，可以在不同的训练数据子集上并行构建模型。袋装法不适用于线性预测器，如线性回归。

Boosting 是一种集成技术，能够产生一系列逐渐复杂的模型。它是按顺序工作的，通过基于前一个模型的错误来训练新的模型。每一个训练出来的模型都会有一个与之关联的权重，这个权重是根据该模型在给定数据上表现的好坏来计算的。当最终做出预测时，这些权重决定了某个特定模型对最终结果的影响程度。与 Bagging 不同，Boosting 不太容易进行并行化处理。由于模型是按序列构建的，它无法并行化处理。序列中较早的分类器所犯的错误被视为难以分类的实例。该框架的设计是让序列中后续的模型拾取前一个预测器所犯的错误或误分类，并尝试改进它们。Boosting 通常使用非常弱的分类器，例如决策树桩，即一个只有一个分裂节点和两个叶子的决策树，作为集成中的一部分。关于 Boosting 的一个非常著名的成功案例是 Viola-Jones 人脸检测算法，其中多个弱分类器（决策树桩）被用来找出有效特征。你可以在以下网站上查看更多关于该成功案例的内容：

en.wikipedia.org/wiki/Viola%E2%80%93Jones_object_detection_framework

在本章中，我们将详细研究 Bagging 和 Boosting 方法。我们将在最后的配方中扩展讨论一种特殊类型的 Boosting，叫做梯度 Boosting。我们还将探讨回归和分类问题，并看看如何通过集成学习来解决这些问题。

理解集成学习—Bagging 方法

集成方法属于一种叫做基于委员会学习的方法家族。与其将分类或回归的决策交给单一模型，不如通过一组模型来共同决策。Bagging 是一种著名且广泛使用的集成方法。

Bagging 也叫做自助聚合。只有当我们能够在基础模型中引入变异性时，Bagging 才能发挥作用。也就是说，如果我们能够成功地在基础数据集中引入变异性，就会导致具有轻微变化的模型。

我们利用自助采样为这些模型引入数据集的变异性。自助采样是通过随机抽取给定数据集中的实例，指定次数地进行采样，且可以选择是否有放回。 在 Bagging 中，我们通过自助采样生成 m 个不同的数据集，并为每个数据集构建一个模型。最后，我们对所有模型的输出进行平均，来产生回归问题中的最终预测结果。

假设我们对数据进行 m 次自助采样，我们将得到 m 个模型，也就是 y m 个值，最终的预测结果如下：

对于分类问题，最终的输出是通过投票决定的。假设我们的集成模型中有一百个模型，并且我们有一个二类分类问题，类别标签为{+1,-1}。如果超过 50 个模型预测输出为+1，我们就将预测结果判定为+1。

随机化是另一种可以在模型构建过程中引入变化的技术。例如，可以随机选择每个集成模型的一个属性子集。这样，不同的模型将拥有不同的属性集。这种技术叫做随机子空间方法。

对于非常稳定的模型，Bagging 可能无法取得很好的效果。如果基础分类器对数据的微小变化非常敏感，Bagging 则会非常有效。例如，决策树非常不稳定，未剪枝的决策树是 Bagging 的良好候选模型。但是，如果是一个非常稳定的模型，比如最近邻分类器 K，我们可以利用随机子空间方法，为最近邻方法引入一些不稳定性。

在接下来的方案中，你将学习如何在 K 最近邻算法中应用 Bagging 和随机子空间方法。我们将解决一个分类问题，最终的预测将基于多数投票。

准备工作…

我们将使用 Scikit learn 中的KNeighborsClassifier进行分类，并使用BaggingClassifier来应用 Bagging 原理。我们将通过make_classification便捷函数生成本方案的数据。

如何实现

让我们导入必要的库，并编写一个get_data()函数，为我们提供数据集，以便我们可以处理这个方案：

from sklearn.datasets import make_classification
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import BaggingClassifier
from sklearn.metrics import classification_report
from sklearn.cross_validation import train_test_split

def get_data():
    """
    Make a sample classification dataset
    Returns : Independent variable y, dependent variable x
    """
    no_features = 30
    redundant_features = int(0.1*no_features)
    informative_features = int(0.6*no_features)
    repeated_features = int(0.1*no_features)
    print no_features,redundant_features,informative_features,repeated_features
    x,y = make_classification(n_samples=500,n_features=no_features,flip_y=0.03,\
            n_informative = informative_features, n_redundant = redundant_features \
            ,n_repeated = repeated_features,random_state=7)
    return x,y

让我们开始编写三个函数：

函数build_single_model用给定数据构建一个简单的 K 最近邻模型。

函数build_bagging_model实现了 Bagging 过程。

函数view_model用于检查我们已经构建的模型：

def build_single_model(x,y):
    model = KNeighborsClassifier()
    model.fit(x,y)
    return model

def build_bagging_model(x,y):
	bagging = BaggingClassifier(KNeighborsClassifier(),n_estimators=100,random_state=9 \
             ,max_samples=1.0,max_features=0.7,bootstrap=True,bootstrap_features=True)
	bagging.fit(x,y)
	return bagging

def view_model(model):
    print "\n Sampled attributes in top 10 estimators\n"
    for i,feature_set in  enumerate(model.estimators_features_[0:10]):
        print "estimator %d"%(i+1),feature_set

最后，我们将编写主函数，调用其他函数：

if __name__ == "__main__":
    x,y = get_data()    

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    # Build a single model    
    model = build_single_model(x_train,y_train)
    predicted_y = model.predict(x_train)
    print "\n Single Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    # Build a bag of models
    bagging = build_bagging_model(x_train,y_train)
    predicted_y = bagging.predict(x_train)
    print "\n Bagging Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
	view_model(bagging)

    # Look at the dev set
    predicted_y = model.predict(x_dev)
    print "\n Single Model Accuracy on Dev data\n"
    print classification_report(y_dev,predicted_y)

    print "\n Bagging Model Accuracy on Dev data\n"
    predicted_y = bagging.predict(x_dev)
    print classification_report(y_dev,predicted_y)

它是如何工作的…

我们从主方法开始。首先调用get_data函数，返回一个包含预测变量的矩阵 x 和响应变量的向量 y。让我们来看一下get_data函数：

    no_features = 30
    redundant_features = int(0.1*no_features)
    informative_features = int(0.6*no_features)
    repeated_features = int(0.1*no_features)
 x,y =make_classification(n_samples=500,n_features=no_features,flip_y=0.03,\
n_informative = informative_features, n_redundant = redundant_features \
            ,n_repeated = repeated_features,random_state=7)

看一下传递给make_classification方法的参数。第一个参数是所需的实例数；在这种情况下，我们说需要 500 个实例。第二个参数是每个实例所需的属性数。我们说我们需要30个属性，这由变量no_features定义。第三个参数，flip_y，会随机交换 3%的实例。这是为了在我们的数据中引入一些噪声。下一个参数指定了这30个特征中应该有多少个是足够信息量大的，可以用于我们的分类。我们指定 60%的特征，也就是 30 个中的 18 个应该是有信息量的。下一个参数是关于冗余特征的。这些特征是由信息性特征的线性组合生成的，用于在特征之间引入相关性。最后，重复特征是从信息性特征和冗余特征中随机抽取的重复特征。

让我们使用train_test_split将数据划分为训练集和测试集。我们保留 30%的数据用于测试：

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)

我们再次利用train_test_split将我们的测试数据分为开发集和测试集。

    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

在划分好数据以构建、评估和测试模型后，我们继续构建模型。我们将首先通过调用以下代码构建一个单一模型，使用KNeighborsClassifier：

model = build_single_model(x_train,y_train)

在这个函数内部，我们创建了一个KNeighborsClassifier类型的对象，并将我们的数据进行拟合，如下所示：

def build_single_model(x,y):
    model = KNeighborsClassifier()
    model.fit(x,y)
    return model

如前一节所述，KNearestNeighbor是一个非常稳定的算法。让我们看看这个模型的表现。我们在训练数据上执行预测，并查看我们的模型指标：

    predicted_y = model.predict(x_train)
    print "\n Single Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)

classification_report是 Scikit learn 模块中的一个方便函数。它给出一个包含precision、recall和f1-score的表格：

在350个实例中，我们的精确度为 87%。有了这个结果，让我们继续构建我们的集成模型：

    bagging = build_bagging_model(x_train,y_train)

我们使用训练数据调用build_bagging_model函数，构建一个分类器集合，如下所示：

def build_bagging_model(x,y):
bagging =             BaggingClassifier(KNeighborsClassifier(),n_estimators=100,random_state=9 \
           ,max_samples=1.0,max_features=0.7,bootstrap=True,bootstrap_features=True)
bagging.fit(x,y)
return bagging

在这个方法中，我们调用了BaggingClassifier类。让我们看一下我们传递给这个类的参数，以便初始化它。

第一个参数是底层估计器或模型。通过传递KNeighborClassifier，我们告诉集成分类器我们想要构建一个由KNearestNeighbor分类器组成的集合。下一个参数指定了我们将构建的估计器的数量。在这种情况下，我们说需要100个估计器。random_state参数是随机数生成器使用的种子。为了在不同的运行中保持一致性，我们将其设置为一个整数值。

我们的下一个参数是 max_samples，指定在从输入数据集进行自助采样时，每个估计器选择的实例数。在这种情况下，我们要求集成程序选择所有实例。

接下来，参数 max_features 指定在为估计器进行自助采样时要包含的属性数量。我们假设只选择 70%的属性。因此，在集成中的每个估计器/模型将使用不同的属性子集来构建模型。这就是我们在上一节中介绍的随机空间方法。该函数继续拟合模型并将模型返回给调用函数。

    bagging = build_bagging_model(x_train,y_train)
    predicted_y = bagging.predict(x_train)
    print "\n Bagging Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)

让我们看看模型的准确性：

你可以看到模型指标有了大幅提升。

在测试模型之前，让我们先查看通过调用 view_model 函数分配给不同模型的属性：

    view_model(bagging)

我们打印出前十个模型选择的属性，如下所示：

def view_model(model):
    print "\n Sampled attributes in top 10 estimators\n"
    for i,feature_set in  enumerate(model.estimators_features_[0:10]):
        print "estimator %d"%(i+1),feature_set

从结果中可以看出，我们已经几乎随机地为每个估计器分配了属性。通过这种方式，我们为每个估计器引入了变异性。

让我们继续检查我们单一分类器和估计器集合在开发集中的表现：

    # Look at the dev set
    predicted_y = model.predict(x_dev)
    print "\n Single Model Accuracy on Dev data\n"
    print classification_report(y_dev,predicted_y)

    print "\n Bagging Model Accuracy on Dev data\n"
    predicted_y = bagging.predict(x_dev)
    print classification_report(y_dev,predicted_y)

正如预期的那样，与单一分类器相比，我们的估计器集合在开发集中的表现更好。

还有更多……

正如我们之前所说，对于分类问题，获得最多票数的标签将被视为最终预测。除了投票方案，我们还可以要求组成模型输出标签的预测概率。最终，可以取这些概率的平均值来决定最终的输出标签。在 Scikit 的情况下，API 文档提供了关于如何执行最终预测的详细信息：

'输入样本的预测类别是通过选择具有最高平均预测概率的类别来计算的。如果基础估计器没有实现predict phobia方法，则会使用投票。'

scikit-learn.org/stable/modules/generated/sklearn.ensemble.BaggingClassifier.html

在上一章中，我们讨论了交叉验证。虽然交叉验证看起来与 Bagging 非常相似，但它们在实际使用中是不同的。在交叉验证中，我们创建 K 折，并根据这些折的模型输出来选择模型的参数，就像我们为岭回归选择 alpha 值一样。这样做主要是为了避免在模型构建过程中暴露我们的测试数据。交叉验证可以用于 Bagging，以确定我们需要向 Bagging 模块添加的估计器数量。

然而，Bagging 的一个缺点是我们失去了模型的可解释性。考虑一个经过剪枝的简单决策树，它很容易解释。但一旦我们有了 100 个这样的模型集合，它就变成了一个黑箱。为了提高准确性，我们牺牲了可解释性。

有关 Bagging 的更多信息，请参阅 Leo Breiman 的以下论文：

*Leo Breiman. 1996. Bagging predictors.*Mach. Learn.24, 2 (August 1996), 123-140. DOI=10.1023/A:1018054314350 dx.doi.org/10.1023/A:1…

另请参见

• 使用交叉验证迭代器，请参阅第七章，《机器学习 2》

• 构建决策树解决多类问题，请参阅第六章，《机器学习 1》

理解集成 - 提升方法

Boosting 是一种强大的集成技术。它几乎被用在大多数数据科学应用中。事实上，它是数据科学家工具包中最重要的工具之一。Boosting 技术利用了类似于 Bagging 的一组估计器。但在这里相似性就结束了。在我们深入研究我们的方法之前，让我们快速看一下 Boosting 如何作为一种非常有效的集成技术。

让我们来看一个熟悉的二类分类问题，其中输入是一组预测变量（X），输出是一个响应变量（Y），它可以取0或1作为值。分类问题的输入表示如下：

分类器的任务是找到一个可以近似的函数：

分类器的误分类率定义为：

假设我们构建了一个非常弱的分类器，其错误率略好于随机猜测。在 Boosting 中，我们在略微修改的数据集上构建了一系列弱分类器。我们为每个分类器轻微修改数据，并最终得到 M 个分类器：

最后，它们所有的预测通过加权多数投票进行组合：

这种方法称为 AdaBoost。

权重 alpha 和模型构建的顺序方式是 Boosting 与 Bagging 不同的地方。正如前面提到的，Boosting 在每个分类器上构建了一系列略微修改的数据集上的弱分类器。让我们看看这个微小的数据修改指的是什么。正是从这个修改中我们得出了我们的权重 alpha。

最初对于第一个分类器，m=1，我们将每个实例的权重设置为 1/N，也就是说，如果有一百条记录，每条记录的权重为 0.001。让我们用 w 来表示权重-现在我们有一百个这样的权重：

所有记录现在都有相等的机会被分类器选择。我们构建分类器，并在训练数据上测试它，以获得误分类率。请参考本节前面给出的误分类率公式。我们将稍微修改它，加入权重，如下所示：

其中 abs 表示结果的绝对值。根据这个误差率，我们计算我们的 alpha（模型权重）如下：

其中 epsilon 是一个非常小的值。

假设我们的模型 1 的误差率是 0.3，也就是说，模型能够正确分类 70%的记录。因此，该模型的权重大约是 0.8，这是一个不错的权重。基于此，我们将返回并设置个别记录的权重，如下所示：

正如你所看到的，所有被错误分类的属性的权重都会增加。这增加了被下一个分类器选择的错误分类记录的概率。因此，下一个分类器在序列中选择权重更大的实例并尝试拟合它。通过这种方式，所有未来的分类器将开始集中处理之前分类器错误分类的记录。

这就是提升方法的威力。它能够将多个弱分类器转化为一个强大的集成模型。

让我们看看提升方法如何应用。在我们编写代码的过程中，我们还将看到 AdaBoost 的一个小变种，称为 SAMME。

开始使用…

我们将利用 scikit-learn 的DecisionTreeClassifier类进行分类，并使用AdaBoostClassifier应用提升原理。我们将使用make_classification便捷函数生成数据。

如何实现

让我们导入必要的库，并编写一个函数get_data()，为我们提供一个数据集来实现这个方案。

from sklearn.datasets import make_classification
from sklearn.ensemble import AdaBoostClassifier
from sklearn.metrics import classification_report,zero_one_loss
from sklearn.cross_validation import train_test_split
from sklearn.tree import DecisionTreeClassifier
import numpy as np
import matplotlib.pyplot as plt
import itertools

def get_data():
    """
    Make a sample classification dataset
    Returns : Independent variable y, dependent variable x
    """
    no_features = 30
    redundant_features = int(0.1*no_features)
    informative_features = int(0.6*no_features)
    repeated_features = int(0.1*no_features)
    print no_features,redundant_features,informative_features,repeated_features
    x,y = make_classification(n_samples=500,n_features=no_features,flip_y=0.03,\
            n_informative = informative_features, n_redundant = redundant_features \
            ,n_repeated = repeated_features,random_state=7)
    return x,y

def build_single_model(x,y):
    model = DecisionTreeClassifier()
    model.fit(x,y)
    return model

def build_boosting_model(x,y,no_estimators=20):
    boosting = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1,min_samples_leaf=1),random_state=9 \
    ,n_estimators=no_estimators,algorithm="SAMME")
    boosting.fit(x,y)
    return boosting

def view_model(model):
    print "\n Estimator Weights and Error\n"
    for i,weight in  enumerate(model.estimator_weights_):
        print "estimator %d weight = %0.4f error = %0.4f"%(i+1,weight,model.estimator_errors_[i])

    plt.figure(1)
    plt.title("Model weight vs error")
    plt.xlabel("Weight")
    plt.ylabel("Error")
    plt.plot(model.estimator_weights_,model.estimator_errors_)

def number_estimators_vs_err_rate(x,y,x_dev,y_dev):
    no_estimators = range(20,120,10)
    misclassy_rate = []
    misclassy_rate_dev = []

    for no_estimator in no_estimators:
        boosting = build_boosting_model(x,y,no_estimators=no_estimator)
        predicted_y = boosting.predict(x)
        predicted_y_dev = boosting.predict(x_dev)        
        misclassy_rate.append(zero_one_loss(y,predicted_y))
        misclassy_rate_dev.append(zero_one_loss(y_dev,predicted_y_dev))

    plt.figure(2)
    plt.title("No estimators vs Mis-classification rate")
    plt.xlabel("No of estimators")
    plt.ylabel("Mis-classification rate")
    plt.plot(no_estimators,misclassy_rate,label='Train')
    plt.plot(no_estimators,misclassy_rate_dev,label='Dev')

    plt.show() 

if __name__ == "__main__":
    x,y = get_data()    
    plot_data(x,y)

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    # Build a single model    
    model = build_single_model(x_train,y_train)
    predicted_y = model.predict(x_train)
    print "\n Single Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_train,predicted_y)*100),"%"

    # Build a bag of models
    boosting = build_boosting_model(x_train,y_train, no_estimators=85)
    predicted_y = boosting.predict(x_train)
    print "\n Boosting Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_train,predicted_y)*100),"%"

    view_model(boosting)

    # Look at the dev set
    predicted_y = model.predict(x_dev)
    print "\n Single Model Accuracy on Dev data\n"
    print classification_report(y_dev,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_dev,predicted_y)*100),"%"

    print "\n Boosting Model Accuracy on Dev data\n"
    predicted_y = boosting.predict(x_dev)
    print classification_report(y_dev,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_dev,predicted_y)*100),"%"

    number_estimators_vs_err_rate(x_train,y_train,x_dev,y_dev)

让我们继续并编写以下三个函数：

函数 build_single_model 用于使用给定数据构建一个简单的决策树模型。

函数 build_boosting_model，它实现了提升算法。

函数 view_model，用于检查我们构建的模型。

def build_single_model(x,y):
    model = DecisionTreeClassifier()
    model.fit(x,y)
    return model

def build_boosting_model(x,y,no_estimators=20):
    boosting = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1,min_samples_leaf=1),random_state=9 \
    ,n_estimators=no_estimators,algorithm="SAMME")
    boosting.fit(x,y)
    return boosting

def view_model(model):
    print "\n Estimator Weights and Error\n"
    for i,weight in  enumerate(model.estimator_weights_):
        print "estimator %d weight = %0.4f error = %0.4f"%(i+1,weight,model.estimator_errors_[i])

    plt.figure(1)
    plt.title("Model weight vs error")
    plt.xlabel("Weight")
    plt.ylabel("Error")
    plt.plot(model.estimator_weights_,model.estimator_errors_)

然后我们编写一个名为 number_estimators_vs_err_rate 的函数。我们使用此函数来查看随着集成模型中模型数量的变化，我们的误差率是如何变化的。

def number_estimators_vs_err_rate(x,y,x_dev,y_dev):
    no_estimators = range(20,120,10)
    misclassy_rate = []
    misclassy_rate_dev = []

    for no_estimator in no_estimators:
        boosting = build_boosting_model(x,y,no_estimators=no_estimator)
        predicted_y = boosting.predict(x)
        predicted_y_dev = boosting.predict(x_dev)        
        misclassy_rate.append(zero_one_loss(y,predicted_y))
        misclassy_rate_dev.append(zero_one_loss(y_dev,predicted_y_dev))

    plt.figure(2)
    plt.title("No estimators vs Mis-classification rate")
    plt.xlabel("No of estimators")
    plt.ylabel("Mis-classification rate")
    plt.plot(no_estimators,misclassy_rate,label='Train')
    plt.plot(no_estimators,misclassy_rate_dev,label='Dev')

    plt.show()

最后，我们将编写主函数，它将调用其他函数。

if __name__ == "__main__":
    x,y = get_data()    
    plot_data(x,y)

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    # Build a single model    
    model = build_single_model(x_train,y_train)
    predicted_y = model.predict(x_train)
    print "\n Single Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_train,predicted_y)*100),"%"

    # Build a bag of models
    boosting = build_boosting_model(x_train,y_train, no_estimators=85)
    predicted_y = boosting.predict(x_train)
    print "\n Boosting Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_train,predicted_y)*100),"%"

    view_model(boosting)

    # Look at the dev set
    predicted_y = model.predict(x_dev)
    print "\n Single Model Accuracy on Dev data\n"
    print classification_report(y_dev,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_dev,predicted_y)*100),"%"

    print "\n Boosting Model Accuracy on Dev data\n"
    predicted_y = boosting.predict(x_dev)
    print classification_report(y_dev,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_dev,predicted_y)*100),"%"

    number_estimators_vs_err_rate(x_train,y_train,x_dev,y_dev)

它是如何工作的…

让我们从主方法开始。我们首先调用get_data函数返回数据集，其中 x 为预测变量矩阵，y 为响应变量向量。让我们深入了解get_data函数：

    no_features = 30
    redundant_features = int(0.1*no_features)
    informative_features = int(0.6*no_features)
    repeated_features = int(0.1*no_features)
 x,y =make_classification(n_samples=500,n_features=no_features,flip_y=0.03,\
n_informative = informative_features, n_redundant = redundant_features \
            ,n_repeated = repeated_features,random_state=7)

查看传递给make_classification方法的参数。第一个参数是所需的实例数量；在这种情况下，我们说我们需要 500 个实例。第二个参数给出了每个实例所需的属性数量。我们说我们需要 30 个属性，正如no_features变量所定义的那样。第三个参数flip_y，随机交换 3%的实例。这是为了在数据中引入一些噪音。接下来的参数指定了这些 30 个特征中应该有多少个是足够有用的，可以用于分类。我们指定 60%的特征，即 30 个特征中的 18 个应该是有信息量的。接下来的参数是冗余特征。这些特征是通过有信息特征的线性组合生成的，用来引入特征之间的相关性。最后，重复特征是从有信息特征和冗余特征中随机选取的重复特征。

让我们使用train_test_split将数据分为训练集和测试集。我们将 30%的数据保留用于测试。

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)

再次，我们利用train_test_split将我们的测试数据分为开发集和测试集。

    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

在将数据分为用于构建、评估和测试模型的部分后，我们开始构建我们的模型。

让我们先拟合一棵单一的决策树，并查看该树在训练集上的表现：

    # Build a single model    
    model = build_single_model(x_train,y_train)

我们通过调用build_single_model函数，传入预测变量和响应变量来构建模型。在这个过程中，我们拟合一棵单一的决策树，并将树返回给调用函数。

def build_single_model(x,y):
    model = DecisionTreeClassifier()
    model.fit(x,y)
    return model

让我们使用classification_report评估模型的好坏，这是来自 Scikit learn 的一个实用函数，它显示一组指标，包括precision（精确度）、recall（召回率）和f1-score；我们还会显示误分类率。

    predicted_y = model.predict(x_train)
    print "\n Single Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    print "Fraction of misclassfication =     
           %0.2f"%(zero_one_loss(y_train,predicted_y)*100),"%"

如你所见，我们的决策树模型完美地拟合了数据——我们的误分类率为 0。在我们在开发集上测试该模型之前，让我们构建我们的集成模型：

    # Build a bag of models
    boosting = build_boosting_model(x_train,y_train, no_estimators=85)

使用build_boosting_model方法，我们按如下方式构建我们的集成模型：

    boosting = AdaBoostClassifier(DecisionTreeClassifier(max_depth=1,min_samples_leaf=1),random_state=9 \
    ,n_estimators=no_estimators,algorithm="SAMME")
    boosting.fit(x,y)

我们利用 Scikit learn 中的AdaBoostClassifier构建我们的提升集成。我们使用以下参数实例化该类：

估计器——在我们的案例中，我们说我们想要构建一个决策树的集成。因此，我们传入DecisionTreeClassifier对象。

max_depth——我们不希望在集成中使用完全生长的树木。我们只需要树桩——只有两个叶子节点和一个分割节点的树。因此，我们将max_depth参数设置为 1。

使用n_estimators参数，我们指定要生成的树木数量；在此案例中，我们将生成 86 棵树。

最后，我们有一个参数叫做 algorithm，它被设置为SAMME。SAMME代表逐阶段加法建模，使用多类指数损失函数。SAMME是对 AdaBoost 算法的改进。它试图将更多的权重分配给误分类的记录。模型权重α是SAMME与 AdaBoost 的区别所在。

我们在前面的公式中忽略了常数 0.5。让我们来看一下新的添加项：log(K-1)。如果 K=2，那么前面的公式就简化为 AdaBoost。在这里，K 是响应变量中的类别数。对于二分类问题，SAMME 就会简化为 AdaBoost，正如前面所述。

让我们拟合模型，并将其返回给调用函数。我们在训练数据集上运行该模型，再次查看模型的表现：

    predicted_y = boosting.predict(x_train)
    print "\n Boosting Model Accuracy on training data\n"
    print classification_report(y_train,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_train,predicted_y)*100),"%"

结果与我们原始模型的表现没有太大不同。我们已经正确分类了几乎 98%的记录。

在对开发集进行测试之前，让我们首先查看我们构建的 Boosting 集成模型：

    view_model(boosting)

在 view_model 内部，我们首先打印出分配给每个分类器的权重：

    print "\n Estimator Weights and Error\n"
    for i,weight in  enumerate(model.estimator_weights_):
        print "estimator %d weight = %0.4f error = %0.4f"%(i+1,weight,model.estimator_errors_[i])

这里我们展示了前 20 个集成的权重。根据它们的误分类率，我们为这些估计器分配了不同的权重。

让我们继续绘制一个图表，显示估计器权重与每个估计器所产生错误之间的关系：

    plt.figure(1)
    plt.title("Model weight vs error")
    plt.xlabel("Weight")
    plt.ylabel("Error")
    plt.plot(model.estimator_weights_,model.estimator_errors_)

正如你所看到的，正确分类的模型被分配了比错误率较高的模型更多的权重。

现在让我们看看单一决策树和集成决策树在开发数据上的表现：

    # Look at the dev set
    predicted_y = model.predict(x_dev)
    print "\n Single Model Accuracy on Dev data\n"
    print classification_report(y_dev,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_dev,predicted_y)*100),"%"

    print "\n Boosting Model Accuracy on Dev data\n"
    predicted_y = boosting.predict(x_dev)
    print classification_report(y_dev,predicted_y)
    print "Fraction of misclassfication = %0.2f"%(zero_one_loss(y_dev,predicted_y)*100),"%"

就像我们对训练数据做的那样，我们打印出分类报告和误分类率：

正如你所看到的，单一决策树表现不佳。尽管它在训练数据上显示了 100%的准确率，但在开发数据上却误分类了近 40%的记录——这是过拟合的迹象。相比之下，Boosting 模型能够更好地拟合开发数据。

我们如何改进 Boosting 模型？其中一种方法是测试训练集中的错误率与我们想要在集成中包含的分类器数量之间的关系。

    number_estimators_vs_err_rate(x_train,y_train,x_dev,y_dev)

以下函数会根据集成的数量递增并绘制错误率：

def number_estimators_vs_err_rate(x,y,x_dev,y_dev):
    no_estimators = range(20,120,10)
    misclassy_rate = []
    misclassy_rate_dev = []

    for no_estimator in no_estimators:
        boosting = build_boosting_model(x,y,no_estimators=no_estimator)
        predicted_y = boosting.predict(x)
        predicted_y_dev = boosting.predict(x_dev)        
        misclassy_rate.append(zero_one_loss(y,predicted_y))
        misclassy_rate_dev.append(zero_one_loss(y_dev,predicted_y_dev))

    plt.figure(2)
    plt.title("No estimators vs Mis-classification rate")
    plt.xlabel("No of estimators")
    plt.ylabel("Mis-classification rate")
    plt.plot(no_estimators,misclassy_rate,label='Train')
    plt.plot(no_estimators,misclassy_rate_dev,label='Dev')

    plt.show()

如你所见，我们声明了一个列表，起始值为 20，结束值为 120，步长为 10。在for循环中，我们将列表中的每个元素作为估计器参数传递给build_boosting_model，然后继续访问模型的错误率。接着我们检查开发集中的错误率。现在我们有两个列表—一个包含训练数据的所有错误率，另一个包含开发数据的错误率。我们将它们一起绘制，x轴是估计器的数量，y轴是开发集和训练集中的错误分类率。

上述图表给出了一个线索，在大约 30 到 40 个估计器时，开发集中的错误率非常低。我们可以进一步实验树模型参数，以获得一个良好的模型。

还有更多内容…

Boosting 方法首次在以下开创性论文中提出：

Freund, Y. & Schapire, R. (1997), 'A decision theoretic generalization of on-line learning and an application to boosting', Journal of Computer and System Sciences 55(1), 119–139.

最初，大多数 Boosting 方法将多类问题简化为二类问题和多个二类问题。以下论文将 AdaBoost 扩展到多类问题：

Multi-class AdaBoost Statistics and Its Interface, Vol. 2, No. 3 (2009), pp. 349-360, doi:10.4310/sii.2009.v2.n3.a8 by Trevor Hastie, Saharon Rosset, Ji Zhu, Hui Zou

本文还介绍了 SAMME，这是我们在配方中使用的方法。

另请参见

• 在第六章中，构建决策树解决多类问题的配方，机器学习 I

• 在第七章中，使用交叉验证迭代器的配方，机器学习 II

• 在第八章中，理解集成方法 – 自助法的配方，模型选择与评估

理解集成方法 – 梯度 Boosting

让我们回顾一下前面配方中解释的 Boosting 算法。在 Boosting 中，我们以逐步的方式拟合加法模型。我们顺序地构建分类器。构建每个分类器后，我们估计分类器的权重/重要性。根据权重/重要性，我们调整训练集中实例的权重。错误分类的实例比正确分类的实例权重更高。我们希望下一个模型能够选择那些错误分类的实例并在其上进行训练。训练集中的那些没有正确拟合的实例会通过这些权重被识别出来。换句话说，这些记录是前一个模型的不足之处。下一个模型试图克服这些不足。

梯度 Boosting 使用梯度而非权重来识别这些不足之处。让我们快速了解如何使用梯度来改进模型。

让我们以一个简单的回归问题为例，假设我们已知所需的预测变量X和响应变量Y，其中Y是一个实数。

梯度提升过程如下：

它从一个非常简单的模型开始，比如均值模型。

预测值仅仅是响应变量的均值。

然后它会继续拟合残差。残差是实际值 y 与预测值 y_hat 之间的差异。

接下来的分类器是在如下数据集上进行训练的：

随后模型会在前一个模型的残差上进行训练，因此算法会继续构建所需数量的模型，最终形成集成模型。

让我们尝试理解为什么我们要在残差上进行训练。现在应该清楚，Boosting 方法构建的是加性模型。假设我们建立两个模型F1(X)和F2(X)来预测Y1。根据加性原理，我们可以将这两个模型结合起来，如下所示：

也就是说，我们结合两个模型的预测结果来预测 Y_1。

等价地，我们可以这样说：

残差是模型未能很好拟合的部分，或者简单来说，残差就是前一个模型的不足之处。因此，我们利用残差来改进模型，即改进前一个模型的不足。基于这个讨论，你可能会好奇为什么这种方法叫做梯度提升（Gradient Boosting）而不是残差提升（Residual Boosting）。

给定一个可微的函数，梯度表示该函数在某些值处的一阶导数。在回归问题中，目标函数为：

其中，F(xi)是我们的回归模型。

线性回归问题是通过最小化前述函数来解决的。我们在F(xi)的值处求该函数的一阶导数，如果用该导数值的负值来更新权重系数，我们将朝着搜索空间中的最小解前进。前述成本函数关于F(xi)的一阶导数是F(xi ) – yi。请参阅以下链接了解推导过程：

zh.wikipedia.org/wiki/%E6%A0%B9%E5%87%BB%E4%BD%8D%E9%9A%8F

F(xi ) – yi，即梯度，是我们残差yi – F(xi)的负值，因此得名梯度提升（Gradient Boosting）。

有了这个理论，我们可以进入梯度提升的具体步骤。

开始使用……

我们将使用波士顿数据集来演示梯度提升。波士顿数据集有 13 个特征和 506 个实例。目标变量是一个实数，即房屋的中位数价格（单位为千美元）。波士顿数据集可以从 UCI 链接下载：

archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.names

我们打算生成 2 次方的多项式特征，并仅考虑交互效应。

如何操作

让我们导入必要的库，并编写一个get_data()函数来提供我们需要使用的数据集，以便完成这个任务：

# Load libraries
from sklearn.datasets import load_boston
from sklearn.cross_validation import train_test_split
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import PolynomialFeatures
import numpy as np
import matplotlib.pyplot as plt

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

def build_model(x,y,n_estimators=500):
    """
    Build a Gradient Boost regression model
    """
    model = GradientBoostingRegressor(n_estimators=n_estimators,verbose=10,\
            subsample=0.7, learning_rate= 0.15,max_depth=3,random_state=77)
    model.fit(x,y)
    return model    

def view_model(model):
    """
    """
    print "\n Training scores"
    print "======================\n"
    for i,score in enumerate(model.train_score_):
        print "\tEstimator %d score %0.3f"%(i+1,score)

    plt.cla()
    plt.figure(1)
    plt.plot(range(1,model.estimators_.shape[0]+1),model.train_score_)
    plt.xlabel("Model Sequence")
    plt.ylabel("Model Score")
    plt.show()

    print "\n Feature Importance"
    print "======================\n"
    for i,score in enumerate(model.feature_importances_):
        print "\tFeature %d Importance %0.3f"%(i+1,score)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\tMean squared error = %0.2f"%(mean_squared_error(true_y,predicted_y))

if __name__ == "__main__":

    x,y = get_data()

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(2,interaction_only=True)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_dev_poly   = poly_features.transform(x_dev)

    # Build model with polynomial features
    model_poly = build_model(x_train_poly,y_train)
    predicted_y = model_poly.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial features)\n"
    model_worth(y_train,predicted_y)  

    # View model details
    view_model(model_poly)

    # Apply the model on dev set
    predicted_y = model_poly.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

    # Apply the model on Test set
    x_test_poly = poly_features.transform(x_test)
    predicted_y = model_poly.predict(x_test_poly)

    print "\n Model Performance in Test set  (Polynomial features)\n"
    model_worth(y_test,predicted_y)  

让我们编写以下三个函数。

build_model 函数实现了梯度提升算法。

view_model 和 model_worth 函数，用于检查我们构建的模型：

def build_model(x,y,n_estimators=500):
    """
    Build a Gradient Boost regression model
    """
    model = GradientBoostingRegressor(n_estimators=n_estimators,verbose=10,\
            subsample=0.7, learning_rate= 0.15,max_depth=3,random_state=77)
    model.fit(x,y)
    return model    

def view_model(model):
    """
    """
    print "\n Training scores"
    print "======================\n"
    for i,score in enumerate(model.train_score_):
        print "\tEstimator %d score %0.3f"%(i+1,score)

    plt.cla()
    plt.figure(1)
    plt.plot(range(1,model.estimators_.shape[0]+1),model.train_score_)
    plt.xlabel("Model Sequence")
    plt.ylabel("Model Score")
    plt.show()

    print "\n Feature Importance"
    print "======================\n"
    for i,score in enumerate(model.feature_importances_):
        print "\tFeature %d Importance %0.3f"%(i+1,score)

def model_worth(true_y,predicted_y):
    """
    Evaluate the model
    """
    print "\tMean squared error = %0.2f"%(mean_squared_error(true_y,predicted_y))

最后，我们编写主函数，该函数将调用其他函数：

if __name__ == "__main__":

    x,y = get_data()

    # Divide the data into Train, dev and test    
    x_train,x_test_all,y_train,y_test_all = train_test_split(x,y,test_size = 0.3,random_state=9)
    x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

    #Prepare some polynomial features
    poly_features = PolynomialFeatures(2,interaction_only=True)
    poly_features.fit(x_train)
    x_train_poly = poly_features.transform(x_train)
    x_dev_poly   = poly_features.transform(x_dev)

    # Build model with polynomial features
    model_poly = build_model(x_train_poly,y_train)
    predicted_y = model_poly.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial features)\n"
    model_worth(y_train,predicted_y)  

    # View model details
    view_model(model_poly)

    # Apply the model on dev set
    predicted_y = model_poly.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

    # Apply the model on Test set
    x_test_poly = poly_features.transform(x_test)
    predicted_y = model_poly.predict(x_test_poly)

    print "\n Model Performance in Test set  (Polynomial features)\n"
    model_worth(y_test,predicted_y)  

如何运作…

让我们从主模块开始，跟随代码进行操作。我们使用 get_data 函数加载预测变量 x 和响应变量 y：

def get_data():
    """
    Return boston dataset
    as x - predictor and
    y - response variable
    """
    data = load_boston()
    x    = data['data']
    y    = data['target']
    return x,y    

该函数调用 Scikit learn 的便捷函数load_boston()来获取波士顿房价数据集，并将其作为 numpy 数组加载。

我们继续使用 Scikit 库中的 train_test_split 函数将数据划分为训练集和测试集。我们保留 30%的数据集用于测试。

x_train,x_test_all,y_train,y_test_all = 
train_test_split(x,y,test_size = 0.3,random_state=9)

从这 30%的数据中，我们在下一行再次提取开发集：

x_dev,x_test,y_dev,y_test = train_test_split(x_test_all,y_test_all,test_size=0.3,random_state=9)

我们接下来构建多项式特征如下：

poly_features = PolynomialFeatures(interaction_only=True)
poly_features.fit(x_train)

如你所见，我们已将 interaction_only 设置为 True。将 interaction_only 设置为 True 时，给定 x1 和 x2 属性时，只会创建 x1*x2 的交互项，而不会创建 x1 的平方和 x2 的平方，假设多项式的度数为 2。默认的度数是 2。

x_train_poly = poly_features.transform(x_train)
x_dev_poly = poly_features.transform(x_dev)
x_test_poly = poly_features.transform(x_test)

使用 transform 函数，我们将训练集、开发集和测试集转换为包含多项式特征的数据集：

让我们继续构建我们的模型：

    # Build model with polynomial features
    model_poly = build_model(x_train_poly,y_train)

在 build_model 函数内部，我们按如下方式实例化 GradientBoostingRegressor 类：

    model = GradientBoostingRegressor(n_estimators=n_estimators,verbose=10,\
            subsample=0.7, learning_rate= 0.15,max_depth=3,random_state=77)

让我们看看这些参数。第一个参数是集成模型的数量。第二个参数是 verbose——当设置为大于 1 的数字时，它会在每个模型（在此为树）的构建过程中打印进度。下一个参数是 subsample，它决定了模型将使用的训练数据百分比。在本例中，0.7 表示我们将使用 70%的训练数据集。下一个参数是学习率。它是一个收缩参数，用于控制每棵树的贡献。接下来的参数 max_depth 决定了构建树的大小。random_state 参数是随机数生成器使用的种子。为了在不同的运行中保持一致性，我们将其设置为一个整数值。

由于我们将 verbose 参数设置为大于 1，在拟合模型时，我们会在每次模型迭代过程中看到以下结果：

如你所见，训练损失随着每次迭代而减少。第四列是袋外改进得分。在子采样中，我们仅选择了数据集的 70%；袋外得分是用剩下的 30%计算的。与前一个模型相比，损失有所改善。例如，在第二次迭代中，相较于第一次迭代构建的模型，我们有 10.32 的改进。

让我们继续检查集成模型在训练数据上的表现：

    predicted_y = model_poly.predict(x_train_poly)
    print "\n Model Performance in Training set (Polynomial features)\n"
    model_worth(y_train,predicted_y)  

如你所见，我们的提升集成模型已经完美地拟合了训练数据。

model_worth 函数打印出模型的更多细节，具体如下：

在详细输出中看到的每个不同模型的得分，都作为属性存储在模型对象中，可以通过以下方式检索：

print "\n Training scores"
print "======================\n"
for i,score in enumerate(model.train_score_):
print "\tEstimator %d score %0.3f"%(i+1,score)

让我们在图表中展示这个结果：

plt.cla()
plt.figure(1)
plt.plot(range(1,model.estimators_.shape[0]+1),model.train_score_)
plt.xlabel("Model Sequence")
plt.ylabel("Model Score")
plt.show()

x 轴表示模型编号，y 轴显示训练得分。记住，提升是一个顺序过程，每个模型都是对前一个模型的改进。

如图所示，均方误差（即模型得分）随着每一个后续模型的增加而减小。

最后，我们还可以看到与每个特征相关的重要性：

    print "\n Feature Importance"
    print "======================\n"
    for i,score in enumerate(model.feature_importances_):
        print "\tFeature %d Importance %0.3f"%(i+1,score)

让我们看看各个特征之间的堆叠情况。

梯度提升将特征选择和模型构建统一为一个操作。它可以自然地发现特征之间的非线性关系。请参考以下论文，了解如何将梯度提升用于特征选择：

Zhixiang Xu, Gao Huang, Kilian Q. Weinberger, 和 Alice X. Zheng. 2014. 梯度提升特征选择。在 第 20 届 ACM SIGKDD 国际知识发现与数据挖掘大会论文集(KDD '14)。ACM, 纽约, NY, USA, 522-531。

让我们将开发数据应用到模型中并查看其表现：

    # Apply the model on dev set
    predicted_y = model_poly.predict(x_dev_poly)
    print "\n Model Performance in Dev set  (Polynomial features)\n"
    model_worth(y_dev,predicted_y)  

最后，我们来看一下测试集上的表现。

如你所见，与开发集相比，我们的集成模型在测试集上的表现极为出色。

还有更多内容…

欲了解更多关于梯度提升的信息，请参考以下论文：

Friedman, J. H. (2001). 贪婪函数逼近：一种梯度提升机。统计年鉴，第 1189–1232 页。

在这份报告中，我们用平方损失函数解释了梯度提升。然而，梯度提升应该被视为一个框架，而不是一个方法。任何可微分的损失函数都可以在这个框架中使用。用户可以选择任何学习方法和可微损失函数，并将其应用于梯度提升框架。

Scikit Learn 还提供了一种用于分类的梯度提升方法，称为 GradientBoostingClassifier。

scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html

另见

• 理解集成方法, Bagging 方法 配方见 第八章，模型选择与评估

• 理解集成方法，AdaBoost 增强方法 配方见 第八章，模型选择与评估

• 使用回归预测实值数 配方见 第七章，机器学习 II

• 使用 LASSO 回归进行变量选择 配方见 第七章，机器学习 II

• 使用交叉验证迭代器 配方见 第七章，机器学习 II