2023年春季英特尔oneAPI校园黑客松竞赛: 淡水质量预测

淡水质量预测

淡水是我们最重要和最稀缺的自然资源之一，仅占地球总水量的 3%。它几乎触及我们日常生活的方方面面，从饮用、游泳和沐浴到生产食物、电力和我们每天使用的产品。获得安全卫生的供水不仅对人类生活至关重要，而且对正在遭受干旱、污染和气温升高影响的周边生态系统的生存也至关重要。因此，利用机器学习技术预测淡水质量是一项至关重要的任务，关乎人类的生存安全和生态系统的健康。

1. 业务理解

机器学习技术将被用于分析历史数据，预测未来淡水需求和供应的趋势，以便人们能够更好地规划淡水的利用和分配，确保淡水的供应能够满足人类和生态系统的需求。预测结果将帮助人们制定更有效的水资源管理政策，缓解淡水资源短缺的问题，同时确保淡水的供应符合卫生和安全标准。

2. 数据预处理

数据来自2023年春季英特尔oneAPI校园黑客松竞赛，首先引入 Python 相关模块包，包括 scikit-learn、xgboost,  pandas、numpy、matplotlib 等核心组件,代码如下:

其中，sklearn.preprocessing 组件主要是对数据预处理相关,例如可使用 LabelBinarizer对标签进行数值化;应用 StandardScaler 对样本特征进行归一化处理等;pandas 主要是对各种格式数据进行高效读取、检查、特征增删，以及实现数据格式转换等;numpy 是科学计算的库,主要用于高维向量的数学函数运算;matplotlib 是对数据进行可视化分析的常见工具之一。

SQL import os import xgboost from xgboost import XGBClassifier import time import warnings import numpy as np import seaborn as sns import matplotlib.pyplot as plt import plotly.io as pio import plotly.graph_objects as go from sklearn.utils import resample from sklearn.model_selection import train_test_split, StratifiedKFold, GridSearchCV, RandomizedSearchCV from sklearn.preprocessing import RobustScaler from sklearn.metrics import roc_auc_score, roc_curve, auc, accuracy_score, f1_score from sklearn.preprocessing import StandardScaler import sklearn from sklearn.metrics import precision_recall_curve, average_precision_score

为了方便查出输出结果，对运行环境进行设置,使其支持中文显示，并配置 pandas 最大列数为 100,避免其对列数进行隐藏,这样利于后续查看各样本列的情况,代码如下:

SQL os.environ['NLSLANG']='SIMPLIFIED CHINESE_CHINA.UTF8' pd.set_option( 'display.max_columns', 100)

我们使用intel加速组件，来加速训练和预测速度。

1. daal4py组件

可将训练好的xgboost转换为 daal4py 模型，以便进一步改进预测时间性能, 利用底层的英特尔® 高级矢量扩展指令集（英特尔® AVX-512）硬件，最大限度地提高英特尔® 至强® 处理器上的梯度提升性能。

1. modin

Modin 使用Ray、Dask或Unidist提供一种轻松的方式来加速您pandas 、脚本和库。与其他分布式 DataFrame 库不同，Modin 提供与现有 pandas 代码的无缝集成和兼容性。即使使用 DataFrame 构造函数也是相同的。

1. sklearnex

借助面向 Scikit-learn* 的英特尔® 扩展，加速 Scikit-learn ，并且仍然完全符合所有 Scikit-Learn API 和算法。英特尔® Extension for Scikit-learn* 是一款免费软件 AI 加速器，可为各种应用程序带来超过 10-100 倍的加速。

面向 Scikit-learn* 的英特尔® 扩展为您提供了一种加速现有 scikit-learn 代码的方法。加速是通过打补丁实现的：用扩展提供的优化版本替换现有的 scikit-learn 算法。

代码如下：

SQL import daal4py as d4p import modin.pandas as pd from modin.config import Engine Engine.put("dask") from sklearnex import patch_sklearn patch_sklearn()

2.1 数据探索及特征选择

下载2023年春季英特尔oneAPI校园黑客松竞赛的数据。

SQL wget https://filerepo.idzcn.com/hack2023/datasetab75fb3.zip

使用pandas将用户相关的特征读取进来，转化为DataFrame对象，并将样本数量字段数量、数据类型，非空值数量等详细信息进行显示

Python # Read data df = pd.read_csv('./dataset.csv') print("Data shape: {}\n".format(df.shape)) display(df.head()) df.info()

运行后，其结果输出如下：

Python <class 'pandas.core.frame.DataFrame'> RangeIndex: 5956842 entries, 0 to 5956841 Data columns (total 24 columns):  #   Column                  Dtype   ---  ------                  -----    0   Index                   int64    1   pH                      float64  2   Iron                    float64  3   Nitrate                 float64  4   Chloride                float64  5   Lead                    float64  6   Zinc                    float64  7   Color                   object  8   Turbidity               float64  9   Fluoride                float64  10  Copper                  float64  11  Odor                    float64  12  Sulfate                 float64  13  Conductivity            float64  14  Chlorine                float64  15  Manganese               float64  16  Total Dissolved Solids  float64  17  Source                  object  18  Water Temperature       float64  19  Air Temperature         float64  20  Month                   object  21  Day                     float64  22  Time of Day             float64  23  Target                  int64   dtypes: float64(19), int64(2), object(3) memory usage: 1.1+ GB

可以看到此 CSV 文件中含有 5956842 条数据,有24 列,浮点型(float64)占 19 列,整型(int64)字段占 2 列,对象型(object)占 3 列，总共占用内存空间大约 1.1GB。为了方便对前述各项指标进行质量检查,包括对数据的分布情况、数据类型、最大值、最小值、均值、标准差等,另外重点查看变量的有效记录数的比例、完整百分比、有效记录数等作为字段筛选的依据,代码如下:

Python desc=df.describe() skewness=df.skew() kurtosis=df.kurtosis() desc=desc.append([skewness, kurtosis],ignore_index=True) idx=pd.Series(['count','mean','std','min','25%','50%','75%','max','skewness','kurtosis'],name='Summary Statistic') desc=pd.concat([desc,idx],1).set_index('Summary Statistic') display(desc.style.format('{:,.3f}').background_gradient(subset=(desc.index[1:],desc.columns[:]),cmap='OrRd'))

代码运行后,将数字型相关(21 列)字段进行统计,其输出结果下图所示。

上图可以看到，部分字段偏差值较高，比如Iron，Lead等，需要进一步做处理。

将Color， Source，Month等英文字段进行因子化（Farcotize）处理为数字型变量，代码如下：

Python display(df.head()) factor = pd.factorize(df['Color']) print(factor) df.Color = factor[0] factor = pd.factorize(df['Source']) print(factor) df.Source = factor[0] factor = pd.factorize(df['Month']) print(factor) df.Month = factor[0] df

2.2 相关性分析

首先，使用pandas库中的corr()函数计算数据框df中各个列（除了Target列）与Target列的相关系数，并使用abs()函数将其转换为绝对值，然后使用sort_values()函数按照相关系数的大小进行排序，得到一个Series类型的bar变量，其中每个元素代表一个特征与Target的相关系数的绝对值。

然后，使用matplotlib库中的bar()函数画出一个条形图，其中x轴代表各个特征，y轴代表它们与Target的相关系数的绝对值。条形图中每个条形的高度表示对应特征与Target的相关系数的绝对值大小。

接下来，使用numpy库中的arange()函数生成一个长度为数据框df中列数的等差数列，然后将其作为x轴上各个条形的位置坐标。

然后，使用rcParams.update()函数更新matplotlib库的参数，其中'figure.figsize'参数用于设置x轴间距。

接下来，使用xticks()函数设置x轴上各个条形的标签，并设置旋转角度和字体大小。

最后，使用show()函数显示图形。

Python # 相关性分析 bar = df.corr()['Target'].abs().sort_values(ascending=False)[1:] plt.bar(bar.index, bar, width=0.5) # 设置figsize的大小 pos_list = np.arange(len(df.columns)) params = {     'figure.figsize': '20, 10' } plt.rcParams.update(params) plt.xticks(bar.index, bar.index, rotation=-60, fontsize=10) plt.show()

 

从上图可以看出，'Index', 'Day', 'Time of Day', 'Month', 'Water Temperature', 'Source', 'Conductivity', 'Air Temperature'这些列相关度很低，需要删除不相关的列

Python  # 删除不相关的列    df = df.drop(     columns=['Index', 'Day', 'Time of Day', 'Month', 'Water Temperature', 'Source', 'Conductivity', 'Air Temperature'])

2.3 缺失值、重复值、偏差值处理

查看样本缺失值、重复值情况

Python display(df.isna().sum()) missing = df.isna().sum().sum() duplicates = df.duplicated().sum() print("\nThere are {:,.0f} missing values in the data.".format(missing)) print("There are {:,.0f} duplicate records in the data.".format(duplicates))

输出结果如下：

Python pH                        116054 Iron                       39753 Nitrate                   105725 Chloride                  175531 Lead                       26909 Zinc                      156126 Color                          0 Turbidity                  49815 Fluoride                  189156 Copper                    199402 Odor                      178891 Sulfate                   197418 Chlorine                   57825 Manganese                 109583 Total Dissolved Solids      1670 Target                         0 dtype: int64 There are 1,603,858 missing values in the data. There are 143,032 duplicate records in the data.

数据集重复值、缺失值比较多，需要进一步处理。对所有特征值为空的样本以0填充，删除重复值，代码如下：

Python df = df.fillna(0) df = df.drop_duplicates()

下面进行偏差值处理，删除偏差值比较大的列。

首先，使用scipy.stats库中的pearsonr()函数计算数据框df中各个列与最后一列的皮尔逊相关系数，并检查每个特征与最后一列的相关性是否显著（p值小于0.05），如果不显著，则将该特征从数据框中删除。

然后，检查每个特征的缺失值比例是否超过20%，如果超过，则将该特征从数据框中删除。

接下来，检查每个数值型特征的方差是否小于等于0.1，如果小于等于0.1，则将该特征从数据框中删除。

最后，输出处理后的数据框df中剩余的特征及其数量。

Python from scipy.stats import pearsonr variables = df.columns df = df var = df.var() numeric = df.columns df = df.fillna(0) for i in range(0, len(var) - 1):     if var[i] <= 0.1:  # 方差大于10％         print(variables[i])         df = df.drop(numeric[i], 1) variables = df.columns for i in range(0, len(variables)):     x = df[variables[i]]     y = df[variables[-1]]     if pearsonr(x, y)[1] > 0.05:         print(variables[i])         df = df.drop(variables[i], 1) variables = df.columns print(variables) print(len(variables))

Python Lead Index(['pH', 'Iron', 'Nitrate', 'Chloride', 'Zinc', 'Color', 'Turbidity',        'Fluoride', 'Copper', 'Odor', 'Sulfate', 'Chlorine', 'Manganese',        'Total Dissolved Solids', 'Target'],       dtype='object') 15

处理之后，删除了Lead 列。只剩余偏差值正常的15列

处理之后的结果集如下所示：

Python pH                        0 Iron                      0 Nitrate                   0 Chloride                  0 Zinc                      0 Color                     0 Turbidity                 0 Fluoride                  0 Copper                    0 Odor                      0 Sulfate                   0 Chlorine                  0 Manganese                 0 Total Dissolved Solids    0 Target                    0 dtype: int64 There are 0 missing values in the data. There are 0 duplicate records in the data.

2.4 数据不平衡处理

Python print(df.Target.value_counts()) target = df.Target.value_counts() target.rename(index={1: 'state 1', 0: 'state o'}, inplace=True) plt.pie(target, [0, 0.05], target.index, autopct='%1.1f%%') plt.show()

0    4057860

1    1754678

Name: Target, dtype: int64

使用下采样对数据进行平衡处理

首先，将数据框df中的特征和目标变量分别赋值给X和y。

然后，使用imblearn库中的RandomUnderSampler()函数进行下采样，使得正负样本比例接近1:1。这里使用了随机数种子（random_state=21）来确保每次采样的结果相同。

接下来，使用train_test_split()函数将X和y分割成训练集和测试集，其中测试集占总数据的20%。同时，使用StandardScaler()函数对训练集和测试集进行标准化处理，使得每个特征的均值为0，标准差为1，以提高模型的性能。

最后，输出训练集和测试集的形状。

Python from imblearn.under_sampling import RandomUnderSampler import datetime X = df.iloc[:, 0:len(df.columns.tolist()) - 1].values y = df.iloc[:, len(df.columns.tolist()) - 1].values # # 下采样 under_sampler = RandomUnderSampler(random_state=21) X, y = under_sampler.fit_resample(X, y) X = df.drop('Target', axis=1) y = df['Target'] X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=21) scaler = StandardScaler() X_train_scaled = scaler.fit_transform(X_train) X_test_scaled = scaler.transform(X_test) print("Train Shape: {}".format(X_train_scaled.shape)) print("Test Shape: {}".format(X_test_scaled.shape)) X_train, X_test = X_train_scaled, X_test_scaled

3. XGBoost算法

数据建模

模型采用XGBoost，并使用随机网格搜索（RandomizedSearchCV）进行优化，模型参数定义如下

Python param_grid = {     'max_depth': [10, 15, 20],     "gamma": [0, 1, 2], # -> 0     "subsample": [0.9, 1], # -> 1     "colsample_bytree": [0.3, 0.5, 1], # -> 1     'min_child_weight': [4, 6, 8], # -> 6     "n_estimators": [10、50, 80, 100], # -> 80     "alpha": [3, 4, 5] # -> 4 } scorers = {     'precision_score': make_scorer(precision_score),     'recall_score': make_scorer(recall_score),     'accuracy_score': make_scorer(accuracy_score),     'f1_score': make_scorer(f1_score),     'roc_auc_score': make_scorer(roc_auc_score), } # {'subsample': 0.5, 'n_estimators': 10, 'min_child_weight': 5, 'max_depth': 12, 'gamma': 0, 'colsample_bytree': 1, 'alpha': 4} # {'subsample': 0.6, 'n_estimators': 5, 'min_child_weight': 2, 'max_depth': 12, 'alpha': 4} # {'subsample': 1, 'n_estimators': 6, 'min_child_weight': 4, 'max_depth': 15, 'alpha': 4} # {'n_estimators': 10, 'min_child_weight': 6, 'max_depth': 15} # {'n_estimators': 20, 'min_child_weight': 6, 'max_depth': 20} # {'n_estimators': 80, 'min_child_weight': 6, 'max_depth': 100} # {'n_estimators': 80, 'min_child_weight': 6, 'max_depth': 100} # {'max_depth': 15}

其中网络参数包括

1. 树的最大深度(max_depth),其取值范围是10，15，20;
2. 评估器数量参数(n _estimators),可选取值为 10、50、80、100;
3. gamma参数 （min_split_loss），可选值为0，1，2， 默认是0，分裂节点时，损失函数减小值只有大于等于gamma节点才分裂，gamma值越大，算法越保守，越不容易过拟合，但性能就不一定能保证，需要平衡。
4. colsample_bytree：可选值为0.3, 0.5, 1，默认= 1，列采样率，也就是特征采样率。范围为（0，1]
5. subsample：可选值为0.9, 1，默认值= 1，构建每棵树对样本的采样率，如果设置成0.5，XGBoost会随机选择一半的样本作为训练集。范围：（0,1]
6. min_child_weight：可选值为4, 6, 8，默认值= 1，如果新分裂的节点的样本权重和小于min_child_weight则停止分裂 。这个可以用来减少过拟合，但是也不能太高，会导致欠拟合。范围：[0，∞]

7.alpha（reg_alpha）：可选值为3, 4, 5， 默认= 0，权重的L1正则化项。增加此值将使模型更加保守。

然后，定义XGBoost分类器。

Python xgb = XGBClassifier(     learning_rate=0.1,     n_estimators=15,     max_depth=12,     min_child_weight=6,     gamma=0,     subsample=1,     colsample_bytree=1,     objective='binary:logistic', # 二元分类的逻辑回归，输出概率     nthread=4,     alpha=4,     scale_pos_weight=1,     seed=27)

Python refit_score = "f1_score" start_time = datetime.datetime.now() print(start_time) rd_search = RandomizedSearchCV(xgb, param_grid, n_iter=10, cv=3, refit=refit_score, scoring=scorers, verbose=10, return_train_score=True) rd_search.fit(X_train, y_train) print(rd_search.best_params_) print(rd_search.best_score_) print(rd_search.best_estimator_) print(datetime.datetime.now() - start_time)

 

在网格搜索过程中,模型采用 precision_score 作为评价指标，基于 StratifiedKFold 将数据分为3份,尝试之前定义的参数列表,使用 fit 方法对模型进行训练。训练完成后,使用测试集(X_test)进行验证，并将结果输出，代码如下:

将测试集的预测结果与实际结果合起来，应用 confusion_matrix 构建混淆矩阵，结果如下:

Python y_pred = rd_search.best_estimator_.predict(X_test) # confusion matrix on the test data. print('\nConfusion matrix of Random Forest optimized for {} on the test data:'.format(refit_score)) print(pd.DataFrame(confusion_matrix(y_test, y_pred),                    columns=['pred_neg', 'pred_pos'], index=['neg', 'pos']))

Python {'subsample': 0.5, 'n_estimators': 10, 'min_child_weight': 8, 'max_depth': 15, 'gamma': 0, 'colsample_bytree': 1, 'alpha': 4}

Python Confusion matrix of Random Forest optimized for f1_score on the test data:      pred_neg  pred_pos neg    683276    128350 pos     10000    341131

使用如下代码输出模型的准确率、查全率、AUC 值、F1值等各项性能指标

Python print(accuracy_score(y_test, y_pred)) print(recall_score(y_test, y_pred)) print(roc_auc_score(y_test, rd_search.best_estimator_.predict_proba(X_test)[:, 1])) print(f1_score(y_test, y_pred))

上述代码运行之后的输入结果如下：

 

可以看到，模型的准确性为 88%，查全率为97%，AUC为0.91，F1值为0.83

结果分析与可视化

ROC曲线是一个二元分类器系统在不同阈值下的性能表现的图形化表示，它通过绘制真正例率（TPR）与假正例率（FPR）在各种阈值设置下的表现来创建。下面使用Scikit-learn库对测试集进行分层K折交叉验证。交叉验证将数据分成5个折叠，并在剩余的折叠上进行验证，同时在4个折叠上训练模型。使用的模型是随机搜索CV算法应用于机器学习模型的结果。 ROC曲线分别为每个折叠绘制，使用不同的线型，使用 `linetypes` 变量。平均ROC曲线也以蓝色绘制，阴影区域表示不同折叠之间性能的变化。为每个折叠和平均曲线计算曲线下面积（AUC）。AUC是一个衡量分类器总体性能的指标。

Python from scipy import interp params = {'legend.fontsize': 'x-large',           'figure.figsize': (12, 9),          'axes.labelsize': 'x-large',          'axes.titlesize':'x-large',          'xtick.labelsize':'x-large',          'ytick.labelsize':'x-large'} plt.rcParams.update(params) tprs = [] aucs = [] mean_fpr = np.linspace(0, 1, 100) skf = StratifiedKFold(n_splits=5) linetypes = ['--',':','-.','-','-','O'] i = 0 cv_data =skf.split(X_test, y_test) for train, test in cv_data:     probas_ = rd_search.predict_proba(X_test[test])     # Compute ROC curve and area the curve     fpr, tpr, thresholds = roc_curve(y_test[test], probas_[:, 1])     tprs.append(interp(mean_fpr, fpr, tpr))     tprs[-1][0] = 0.0     roc_auc = auc(fpr, tpr)     aucs.append(roc_auc)     plt.plot(fpr, tpr, lw=1.5,linestyle = linetypes[i], alpha=0.8,              label='ROC fold %d (AUC = %0.2f)' % (i, roc_auc))     i += 1 plt.plot([0, 1], [0, 1], linestyle='--', lw=1, color='r',          label='Chance', alpha=.6) mean_tpr = np.mean(tprs, axis=0) mean_tpr[-1] = 1.0 mean_auc = auc(mean_fpr, mean_tpr) std_auc = np.std(aucs) plt.plot(mean_fpr, mean_tpr, color='b',          label=r'Mean ROC (AUC = %0.2f $\pm$ %0.2f)' % (mean_auc, std_auc),          lw=2, alpha=.8) std_tpr = np.std(tprs, axis=0) tprs_upper = np.minimum(mean_tpr + std_tpr, 1) tprs_lower = np.maximum(mean_tpr - std_tpr, 0) plt.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2,                  label=r'$\pm$ 1 std. dev.') plt.xlim([-0.02, 1.02]) plt.ylim([-0.02, 1.02]) plt.xlabel('FPR',fontsize=20) plt.ylabel('TPR',fontsize=20) # plt.title('ROC') plt.legend(loc="lower right") plt.show()

4. 随机森林算法

数据建模

模型采用随机森林算法，并使用随机网格搜索（RandomizedSearchCV）进行优化，模型参数定义如下

Python param_grid = {     'max_depth': [10, 15, 20],     "gamma": [0, 1, 2], # -> 0     "subsample": [0.9, 1], # -> 1     "colsample_bytree": [0.3, 0.5, 1], # -> 1     'min_child_weight': [4, 6, 8], # -> 6     "n_estimators": [10、50, 80, 100], # -> 80     "alpha": [3, 4, 5] # -> 4 } scorers = {     'precision_score': make_scorer(precision_score),     'recall_score': make_scorer(recall_score),     'accuracy_score': make_scorer(accuracy_score),     'f1_score': make_scorer(f1_score),     'roc_auc_score': make_scorer(roc_auc_score), } classifier = RandomForestClassifier(criterion = 'entropy', oob_score=True,random_state = 42) refit_score = "f1_score" start_time = datetime.datetime.now() print(start_time) rd_search = RandomizedSearchCV(classifier, param_grid, n_iter=10, cv=3, refit=refit_score, scoring=scorers, verbose=10,                                return_train_score=True) rd_search.fit(X_train, y_train) print(rd_search.best_params_) print(rd_search.best_score_) print(rd_search.best_estimator_) print(datetime.datetime.now() - start_time) y_pred = randomForestClf.predict(X_test) print("roc:",roc_auc_score(y_test, randomForestClf.predict_proba(X_test)[:,1])) print("accuracy:",accuracy_score(y_test, y_pred)) print("recall:",recall_score(y_test, y_pred)) print("f1:",f1_score(y_test, y_pred))

上述代码运行之后的输入结果如下：

 

 

可以看到其准确率约为0.88，查全率约为0.97，ROC约为0.91，F1值约为0.83

5. 与其他算法比较

逻辑回归

Python  from sklearn.linear_model import LogisticRegression model = LogisticRegression() model.fit(X_train, y_train)

Python y_pred = model.predict(X_test) print("roc:", roc_auc_score(y_test, randomForestClf.predict_proba(X_test)[:, 1])) print("accuracy:", accuracy_score(y_test, y_pred)) print("recall:", recall_score(y_test, y_pred)) print("f1:", f1_score(y_test, y_pred))

 

可以看到其准确率约为0.79，查全率约为0.72，ROC约为0.91，F1值约为0.67

决策树

Python from sklearn.tree import DecisionTreeClassifier model = DecisionTreeClassifier() model.fit(X_train, y_train) y_pred = model.predict(X_test) print("roc:",roc_auc_score(y_test, model.predict_proba(X_test)[:,1])) print("accuracy:",accuracy_score(y_test, y_pred)) print("recall:",recall_score(y_test, y_pred)) print("f1:",f1_score(y_test, y_pred))

可以看到其准确率约为0.82，查全率约为0.70，AUC约为0.70，F1值约为0.70

 

总结

计算不同算法在准确率、查全率、F1值、AUC值指标，其结果如下表所示

算法	准确率	查全率	F1值	AUC值
逻辑回归	0.79	0.72	0.67	0.91
决策树	0.82	0.70	0.70	0.70
随机森林	0.88	0.97	0.83	0.92
XGBoost	0.88	0.97	0.83	0.91

从表中可以看出，XGBoost和随机森林的F1值和准确率差不多，F1值达到83%，具有较强的实用价值，所以最终采用随机森林和XGBboost作为适用模型。

参考文献

1. 《Python机器学习实战案例》（第二版）第四章 基于分类算法的学习失败预警 - 赵卫东 董亮著
2. 《机器学习》（第二版）第三章 决策树与分类算法 - 赵卫东 董亮著