Interl校企实训 基于机器学习的预测淡水质量

一.预测淡水质量

     1.1 背景描述

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

     1.2 预期解决方案

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

     1.3 数据集

               数据集

     1.4 要求

               需要使用英特尔® oneAPI AI分析工具包中数据分析与机器学习相关的优化库，基于参考资料的方案，进行淡水质量检测的实现。

二.数据预处理

     2.1加载所需包和库

import os
import xgboost
from xgboost import XGBClassifier
import time
import warnings
import modin.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

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

     2.2 运用组件来加速

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

              2.2.1 daal4py组件

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

              2.2.2 modin

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

              2.2.3 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. 数据探索及特征选择

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

数据读取与分析：

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

查看前部分的数据：

上图可以看到，部分字段偏差值较高，比如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()函数显示图形。

# 相关性分析
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'这些列相关度很低，需要删除不相关的列。

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

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

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

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

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

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

from imblearn.under_sampling import RandomUnderSampler
import datetime

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

四.XGBoost算法

     4.1建立模型

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

from sklearn.metrics import make_scorer, precision_score, recall_score, accuracy_score, f1_score, roc_auc_score  
from xgboost import XGBClassifier   
param_grid = {
    'max_depth': [10, 15, 20],
    "gamma": [0, 1, 2], # -> 0
    "subsample": [0.9, 1], # -> 1
    "colsample_bytree": [0.3, 0.5, 1], 
    'min_child_weight': [4, 6, 8], 
    "n_estimators": [10,50, 80, 100], 
    "alpha": [3, 4, 5]
}

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

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

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

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

start = time.time()
y_pred = rd_search.best_estimator_.predict(X_test)
from sklearn.metrics import confusion_matrix
# confusion matrix on the test data.
end = time.time()
print('模型拟合时间:{:.2f}'.format(end-start))
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']))
end = time.time()
print('模型拟合时间:{:.2f}'.format(end-start))

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

输出模型的准确率、查全率、AUC 值、F1值等各项性能指标：

print("accuracy:",accuracy_score(y_test, y_pred))
print("recall:",recall_score(y_test, y_pred))

print("roc:",roc_auc_score(y_test, rd_search.best_estimator_.predict_proba(X_test)[:, 1]))
print("f1:",f1_score(y_test, y_pred))

测试数据集：

#测试
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)
test_df = pd.read_csv('./test_data.csv')
print("Data shape: {}\n".format(df.shape))
display(df.head())
 
df.info()
 
display(test_df.head())
factor = pd.factorize(df['Color'])
print(factor)
test_df.Color = factor[0]
factor = pd.factorize(test_df['Source'])
print(factor)
test_df.Source = factor[0]
factor = pd.factorize(test_df['Month'])
print(factor)
test_df.Month = factor[0]
test_df = test_df.drop(
   columns=['Index', 'Day', 'Time of Day', 'Month', 'Water Temperature', 'Source', 'Conductivity', 'Air Temperature'])
 
test_df = test_df.fillna(0)
test_df = df.drop_duplicates()
 
test_df = test_df.drop(columns=['Lead'])
test_df = test_df.drop(columns=['Target'])
 
X_test = scaler.transform(test_df)
 
import pandas as pd
 
# 读取 test_data.csv 文件
test_df = pd.read_csv('test_data.csv')
 
# 提取 Target 列的值
y_true = test_df['Target']
 
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)
 
# 计算模型在测试集上的f1分数并输出
f1 = f1_score(y_true, y_pred)
print("模型在测试集上的f1分数:", f1)

结果：

五.结果可视化

 结果可视化：

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

使用Scikit-learn库：

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

输出结果：

六.学习总结

总结：

代码首先读取了一个数据集，并对数据进行了一系列的数据预处理和特征工程操作，包括计算偏度和峰度、因子化处理英文字段、相关性分析、删除不相关的列、处理缺失值和重复值、特征选择等操作。在数据准备好之后，进行了数据切分和标准化处理。本模型使用XGBoost，并使用随机网格搜索进行优化,经过10次训练 , f1分数在0.83左右 , 准确性较高.通过使用Intel oneAPI组件中的daal4py进行加速,使得推理时间在0.24s左右.最后评估了模型在测试集上的性能表现，包括 ROC-AUC、准确率、召回率和 F1 分数等指标，所需时间较短,准确性较高.可以看出训练的效果较好。