使用阿里云机器学习PAI做回归预测

在DSW的Sample Notebook中，使用XGBoost回归算法预测房价。

XGBoost算法原理

XGBoost（eXtreme Gradient Boosting）是一种基于决策树集成学习算法的机器学习模型。它是由陈天奇在2016年提出的，经过不断的优化和改进，已成为了在机器学习竞赛中广泛使用的模型之一。XGBoost的特点是在GBDT（Gradient Boosting Decision Tree）算法的基础上，加入了正则化项和二阶导数信息，从而提高了模型的泛化能力和准确性。

XGBoost模型的预测函数可以表示为：

其中，K 表示决策树的数量，f_k(x_i) 表示第 k 棵树对样本 x_i 的预测值。为了使模型更加准确，XGBoost在GBDT的基础上，引入了正则化项和二阶导数信息，构建了如下的目标函数：

其中，l(y_i,\hat{y}_i) 表示预测值和真实值之间的损失函数，\Omega(f_k) 表示对决策树进行正则化的惩罚项。为了求解这个目标函数，XGBoost采用了梯度提升算法和加权最小二乘法进行优化。在每一轮迭代中，XGBoost根据当前模型的预测值和真实值之间的残差，构建一棵新的决策树，并根据正则化项和二阶导数信息对决策树进行剪枝和修剪，最终得到一个泛化能力更强的模型。

XGBoost的优点在于它既具有高准确性，又具有高效性。它能够处理大规模、高维度的数据集，并能够快速训练和预测。因此，它被广泛应用于各种机器学习任务，如分类、回归、排序、推荐等。

XGBoost的参数介绍

XGBoost的全部参数介绍比较详细，以下是常用参数的介绍：

通用超参数:

• booster [default=gbtree]: 选择哪种booster，可选gbtree, gblinear或dart。
• silent [default=0]: 是否在运行时输出运行信息。
• nthread [default to maximum number of threads available if not set]: 用于运行XGBoost的线程数量。

增强器参数:

• eta [default=0.3]: 学习率，控制每一步的权重缩减。
• gamma [default=0]: 控制是否进行节点分裂的参数，只有分裂后损失函数得到显著下降，才会进行分裂。
• max_depth [default=6]: 树的最大深度，避免过拟合。
• min_child_weight [default=1]: 定义一个节点在分裂时所需的最小样本权重和。
• max_delta_step [default=0]: 允许每个树的权重估计的最大增量步长，该参数可以帮助收敛速度，通常不需要调整。
• subsample [default=1]: 训练模型时的子采样率。
• colsample_bytree [default=1]: 在建立树时，对特征进行采样的比例。
• colsample_bylevel [default=1]: 控制每个级别（深度）的列采样比例。
• lambda [default=1]: 控制L2正则化项的权重。
• alpha [default=0]: 控制L1正则化项的权重。
• tree_method [default=auto]: 控制使用何种方法来构建树，可选exact、approx和hist。
• scale_pos_weight [default=1]: 解决样本不平衡问题，将正样本的权重提高。

学习任务参数:

• objective [default=reg:squarederror]: 目标函数，可选reg:squarederror（回归问题）、binary:logistic（二分类问题）或multi:softmax（多分类问题）等。
• num_class [default=1]: 多分类问题的类别数量。
• eval_metric [default according to objective]: 评估指标，根据目标函数自动选择。

以上是常用的参数介绍，XGBoost还有一些高级参数和调试参数，可根据实际需求进行设置。

创建DSW实例

导入三方库

import numpy as np  
import pandas as pd  
%matplotlib inline
import matplotlib.pyplot as plt   
import seaborn as sns
color = sns.color_palette()
sns.set_style('darkgrid')

import warnings
warnings.filterwarnings("ignore")


from scipy import stats
from scipy.stats import norm, skew  

pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x))

数据清洗与预处理

# 加载数据
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')

原始数据都会存在很多问题，不利于分析和训练，所以要清洗和预处理，比如去无关列，去重，缺失值处理，去极端值等。

#保存'Id'列
train_ID = train['Id']
test_ID = test['Id']

#删除'Id'列，因为它对于预测过程是不必要的。
train.drop("Id", axis = 1, inplace = True)
test.drop("Id", axis = 1, inplace = True)

# 删除GrLivArea特征值大于4000且SalePrice特征值小于300000的异常点
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

# 对SalePrice特征使用np.log1p函数进行平滑处理，使其更接近标准正态分布。
train["SalePrice"] = np.log1p(train["SalePrice"])

# 将训练集(train)和测试集(test)按行合并，并删除合并后数据集(all_data)中的SalePrice特征。
ntrain = train.shape[0]
ntest = test.shape[0]
y_train = train.SalePrice.values
all_data = pd.concat((train, test)).reset_index(drop=True)
all_data.drop(['SalePrice'], axis=1, inplace=True)

# 计算all_data中每个特征的缺失率，取出前30个缺失率较高的特征，放在missing_data中。
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})

# 使用 "None" 来填充空值的特征
none_cols = ["PoolQC", "MiscFeature", "Alley", "Fence", "FireplaceQu", "GarageType",
             "GarageFinish", "GarageQual", "GarageCond", "BsmtQual", "BsmtCond",
             "BsmtExposure", "BsmtFinType1", "BsmtFinType2", "MasVnrType"]
for feature in none_cols:
    all_data[feature].fillna('None', inplace=True)

# 使用 0 来填充空值的特征
zero_cols = ["GarageYrBlt", "GarageArea", "GarageCars", "BsmtFinSF1", "BsmtFinSF2",
             "BsmtUnfSF", "TotalBsmtSF", "BsmtFullBath", "BsmtHalfBath", "MasVnrArea"]
for feature in zero_cols:
    all_data[feature].fillna(0, inplace=True)

# 使用中位数来填充空值的特征
all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
    lambda x: x.fillna(x.median()))

# 使用众数来填充空值的特征
mode_cols = ["MSZoning", "Electrical", "KitchenQual", "Exterior1st", "Exterior2nd", "SaleType"]
for feature in mode_cols:
    all_data[feature] = all_data[feature].fillna(all_data[feature].mode()[0])

# 删除 "Utilities" 特征
all_data = all_data.drop(["Utilities"], axis=1)

# 使用 "Typ" 来填充空值的特征
all_data["Functional"] = all

# 统计 all_data 中每个特征列的缺失值比例
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
# 选出所有缺失值比例为零的特征列的索引，并在 all_data_na 中将其从 DataFrame 中删除。
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
# 根据 all_data_na 生成一个 missing_data
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})

# 类型转换
cols = ['MSSubClass', 'OverallCond', 'YrSold', 'MoSold']
all_data[cols] = all_data[cols].astype(str)

# 编码(encoding)
from sklearn.preprocessing import LabelEncoder
for c in all_data.columns:
    if all_data[c].dtype == 'object':
        lbl = LabelEncoder() 
        lbl.fit(list(all_data[c].values)) 
        all_data[c] = lbl.transform(list(all_data[c].values))

print('Shape all_data: {}'.format(all_data.shape))

# 根据相关的行业知识，创建新的feature
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']

# 一共有多少个洗漱间
all_data['TotalBath'] = all_data[['BsmtFullBath','BsmtHalfBath','FullBath','HalfBath']].sum(axis=1)

# 门廊的面积
all_data['TotalPorchSF'] = all_data[['OpenPorchSF','EnclosedPorch','3SsnPorch','ScreenPorch','WoodDeckSF']].sum(axis=1)

# 计算feature的偏度
numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index

skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
print("\n数值类型feature的偏度: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})

# 对偏度大于0，75的进行平滑处理
skewness = skewness[abs(skewness) > 0.75]
print("一共有 {} 个feature需要处理".format(skewness.shape[0]))

from scipy.special import boxcox1p
skewed_features = skewness.index
lam = 0.15
for feat in skewed_features:
    all_data[feat] = boxcox1p(all_data[feat], lam)

all_data = pd.get_dummies(all_data)

# 产生最终的数据集
train = all_data[:ntrain]
test = all_data[ntrain:]

导入建模所需的包

from sklearn.kernel_ridge import KernelRidge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler
from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_squared_error
import xgboost as xgb

编制函数计算模型在训练集上的交叉验证准确率

# 参数：model: 训练好的模型，train: 训练数据，y_train: 目标变量
# 返回值：交叉验证得分列表
def get_accuracy(model,train,y_train):
    n_folds=7
    kf1 = KFold(n_folds, shuffle=True, random_state=42)
    kf_cv_scores = cross_val_score(model,train,y_train,cv=kf1)
    return kf_cv_scores

编制函数计算均方根对数误差

def rmsle(y, y_pred):
    return np.sqrt(mean_squared_error(y, y_pred))

使用XGBoost建立回归模型

model_xgb = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, 
                             learning_rate=0.05, max_depth=3, 
                             min_child_weight=1.7817, n_estimators=2200,
                             reg_alpha=0.4640, reg_lambda=0.8571,
                             subsample=0.5213, 
                             random_state =7, nthread = -1)

训练

model_xgb.fit(train, y_train)

XGBRegressor(base_score=0.5, booster='gbtree', colsample_bylevel=1,
             colsample_bynode=1, colsample_bytree=0.4603,
             enable_categorical=False, gamma=0.0468, gpu_id=-1,
             importance_type=None, interaction_constraints='',
             learning_rate=0.05, max_delta_step=0, max_depth=3,
             min_child_weight=1.7817, missing=nan, monotone_constraints='()',
             n_estimators=2200, n_jobs=8, nthread=-1, num_parallel_tree=1,
             predictor='auto', random_state=7, reg_alpha=0.464,
             reg_lambda=0.8571, scale_pos_weight=1, subsample=0.5213,
             tree_method='exact', validate_parameters=1, verbosity=None)

kf_cv_scores = get_accuracy(model_xgb,train,y_train)
print("Average score: %.2f" % kf_cv_scores.mean())

结果

Average score: 0.92

打印其中一棵树的树结构

tree_struct = xgb.to_graphviz(model_xgb,num_trees=111)
tree_struct