【预测淡水质量】

1.问题描述：

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

2.预期解决方案：

通过参考英特尔的类似实现方案，预测淡水是否可以安全饮用和被依赖淡水的生态系统所使用，从而可以帮助全球水安全和环境可持续性发展。这里分类准确度和推理时间将作为评分的主要依据。

3.数据预处理:

数据来自2023年春季英特尔oneAPI校园黑客松竞赛，首先引入 Python 相关模块包，包括 scikit-learn、xgboost, pandas、numpy、matplotlib 等核心组件。
其中，sklearn.preprocessing 组件主要是对数据预处理相关,例如可使用 LabelBinarizer对标签进行数值化;应用 StandardScaler 对样本特征进行归一化处理等;pandas 主要是对各种格式数据进行高效读取、检查、特征增删，以及实现数据格式转换等;numpy 是科学计算的库,主要用于高维向量的数学函数运算;matplotlib 是对数据进行可视化分析的常见工具之一。

import os
import xgboost
from xgboost import XGBClassifier
import time
import warnings
import pandas as pd
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
import zipfile
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,避免其对列数进行隐藏,这样利于后续查看各样本列的情况。

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

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

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

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

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

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

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

3.1 数据探索及特征选择

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

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

运行后结果：

可以看到此 CSV 文件中含有 5956842 条数据,有24 列,浮点型(float64)占 19 列,整型(int64)字段占 2 列,对象型(object)占 3 列，总共占用内存空间大约 1.1GB。

为了方便对前述各项指标进行质量检查,包括对数据的分布情况、数据类型、最大值、最小值、均值、标准差等,另外重点查看变量的有效记

df = df.apply(pd.to_numeric, errors='coerce')
desc = df.describe()

skewness = df.skew()
kurtosis = df.kurtosis()

desc.loc['skewness'] = skewness
desc.loc['kurtosis'] = kurtosis

idx = pd.Index(['count', 'mean', 'std', 'min', '25%', '50%', '75%', 'max', 'skewness', 'kurtosis'], name='Summary Statistic')
desc = desc.set_index(idx)

desc_styled = desc.style.format('{:,.3f}').background_gradient(subset=pd.IndexSlice['mean':'kurtosis', :], cmap='OrRd')

display(desc_styled)

代码运行后,将数字型相关字段进行统计,其输出结果如下

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

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

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]

3.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()函数显示图形。

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

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

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

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

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))

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

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

下面进行偏差值处理，删除偏差值比较大的列。
首先，使用scipy.stats库中的pearsonr()函数计算数据框df中各个列与最后一列的皮尔逊相关系数，并检查每个特征与最后一列的相关性是否显著（p值小于0.05），如果不显著，则将该特征从数据框中删除。
然后，检查每个特征的缺失值比例是否超过20%，如果超过，则将该特征从数据框中删除。
接下来，检查每个数值型特征的方差是否小于等于0.1，如果小于等于0.1，则将该特征从数据框中删除。
最后，输出处理后的数据框df中剩余的特征及其数量。

from scipy.stats import pearsonr
 
variables = df.columns
df = df
 
var = df.var()
numeric = df.columns
df = df.fillna(df.interpolate())
for i in range(0, len(var) - 1):
    if var[i] <= 0.1:  # 方差大于10％
        print(variables[i])
        df = df.drop(numeric[i],axis=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],axis=1)
 
variables = df.columns
print(variables)
print(len(variables))
print(df.Target.value_counts())

3.4 数据不平衡处理

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()

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

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

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

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

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

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

4. XGBoost算法

4.1数据建模

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

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}

其中网络参数包括

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

然后，定义XGBoost分类器。

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)

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)进行验证，并将结果输出

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']))

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

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

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。

计算推理时间

#测试
from datetime import datetime
# 记录开始时间
inference_start_time = datetime.now()
 
# 模型推理代码
y_pred = rd_search.best_estimator_.predict(X_test)
 
# 计算模型推理时间
inference_time = datetime.now() - inference_start_time
print("模型推理时间:", inference_time)

4.2结果分析与可视化

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

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()

5. 总结

本项目的目标是通过机器学习模型对淡水质量进行预测，并实现高准确度和高泛化性能，以提供准确的水质评估和监测。对数据进行了预处理，包括缺失值处理、异常值处理、特征选择和数据归一化等。确保数据的准确性和一致性。
在多个机器学习算法中，选择了适合分类问题的模型，如决策树、随机森林或支持向量机等。分割数据集为训练集和测试集，并使用交叉验证来评估模型的性能。
迭代优化模型参数，并探索不同特征工程的策略，以提高模型性能。使用准确度和推理时间作为评估指标，准确度反映模型整体分类准确性，推理时间反映了模型在给定的硬件和软件环境下能够实时地进行预测。

在测试集上，模型达到了88.10% 的准确度和0.18 的推理时间，表明模型具备较好的分类能力和实际应用中进行快速推理能力。