【Intel校企合作】淡水质量预测

1. 预测淡水质量

1.1 问题描述

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

1.2 预期解决方案

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

1.3 数据集

数据集点此处下载
数据集特征中包含：
Index：编号
pH：水资源的pH值
Iron：含铁量
Nitrate：硝酸盐含量
Chloride：氯化物含量
Lead：铅含量
Zinc：锌含量
Color：淡水颜色
Turbidity：浑浊度
Fluoride：氟化物含量
Copper：铜含量
Odor：气味
Sulfate Conductivity：硫酸导电率
Chlorine Manganese ：氯锰含量
Total Dissolved Solids ：总溶解固体
Source：来源
Temperature ：温度
Air Temperature ：空气温度
Month ：取水的月份
Day ：取水的日期
Time of Day ：取水的时间
Target：是否为可使用淡水

1.4 要求

需要使用 英特尔® ONEAPI AI分析工具包。

1.5 参考资料

欢迎参考英特尔提供的类似实现方案。
也可以参考课程实践案例。

2. 导入相关包

import datetime
import time
from sklearn.metrics import f1_score, accuracy_score, roc_auc_score, roc_curve, auc
import numpy as np
# import pandas as pd
import sklearn
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split, RandomizedSearchCV
from sklearn.preprocessing import StandardScaler
from xgboost import *
import daal4py as d4p

# oneapi优化(modin,sklearnex)
import modin.pandas as pd
from modin.config import Engine
Engine.put("dask")
 
from sklearnex import patch_sklearn
 
patch_sklearn()

运行结果截图：

Intel 加速组件介绍：

• Modin
  Modin是一个强大的工具，可以帮助用户加速大规模数据分析的处理速度。它提供了与pandas兼容的API，使得用户可以无缝地在Modin和pandas之间切换，同时还提供了并行计算、可扩展性和易于使用等特性。

• Daal4py
  daal4py 是 Intel 提供的一个 Python 库，用于实现快速的数据分析和机器学习算法。

• Scikit-learn
  Scikit-learn组件为使用Scikit-learn进行机器学习的用户提供了显著的性能提升。通过利用Intel的硬件优化技术和跨设备兼容性，这个组件使得处理大规模数据集和运行复杂算法变得更加高效和便捷。

• XGBoost
  XGBoost组件为使用XGBoost进行机器学习的用户提供了更高的性能。通过利用Intel的硬件优化技术和改进数据加载处理方式，这个组件使得处理大规模数据集和运行复杂算法变得更加高效和快速。

使用modin.pandas库时，将默认的引擎更改为dask。modin.pandas是一个用于处理大型数据集的库，它提供了类似于Pandas的功能，但可以更好地利用分布式计算资源。通过将引擎设置为dask，可以利用Dask库提供的并行计算功能来加速数据处理任务。

3. 数据预处理

3.1 读取并查看数据集

df = pd.read_csv('dataset.csv')

print(df.info(), df.describe())

3.2 正负样本数据比例查看

这里用环状图可以直观地看出数据集中正负样本数据比例大致为2:1，后续需要采取相关取样方法来使得样本平衡。

def plot_target(target_col):
    tmp=df[target_col].value_counts(normalize=True)
    target = tmp.rename(index={1:'Target 1',0:'Target 0'})
    wedgeprops = {'width':0.5, 'linewidth':10}
    
    plt.figure(figsize=(6,6))
    plt.pie(list(tmp), labels=target.index,
            startangle=90, autopct='%1.1f%%',wedgeprops=wedgeprops)
    plt.title('Label Distribution', fontsize=16)
    plt.show() 
    
plot_target(target_col='Target')

3.3 数据类型转换

在3.1中查看数据集的时候可以发现，有三个object类型的数据特征（Color、Source、Month）不利于后面进行数据分析，所以我们要先将其进行数据映射转换为数值型数据。

# 定义函数将分类变量映射为数值
def trans(l):
    data = {}
    for i in range(len(l)):
        data[l[i]] = i + 1
    return data

# 将分类变量转换为数值型
d1 = trans(df['Source'].unique())
d2 = trans(df['Month'].unique())
d3 = trans(df['Color'].unique())
df['Source'] = df['Source'].map(d1)
df['Month'] = df['Month'].map(d2)
df['Color'] = df['Color'].map(d3)

# 输出转换后的数据信息
print(df.info(), df.describe())

从图中我们可以看出刚刚object类型的数据特征（Color、Source、Month）已经成功转为float64数值型数据。

3.4 相关性分析

为了方便后面对特征数据的选取，这里用热力图进行相关性分析。

plt.figure(figsize=(30, 20))
sns.heatmap(df.corr(), annot=True, cmap='Blues', xticklabels=1, yticklabels=1)

可以发现，有的数据和Target的相关性很差，所以我们排除小于等于0.1的：

# 筛选相关性大于等于0.1的特征
col = df.corr()['Target']
data = col[col >= 0.1].index
data = df[data]

3.5 缺失值、重复值处理

查看缺失值和重复值

display(df.isna().sum())
missing = df.isna().sum().sum()
duplicates = df.duplicated().sum()
print("\n共{:,.0f} 缺失值 in the data.".format(missing))
print("共{:,.0f} 重复值 in the data.".format(duplicates))

从图中可以看出缺失值还是很多的，所以不可以简单删除。这里我们就采用上下值的平均值来填充缺失值，没有重复值则不管。

# 填充缺失值
data = data.fillna(data.interpolate())
data = data.drop_duplicates()
print(data.isna().sum(),data.isna().sum(),data.duplicated().sum())
print("\n共{:,.0f} 缺失值 in the data.".format(data.isna().sum().sum()))
print("共{:,.0f} 重复值 in the data.".format(data.duplicated().sum()))

结果如下，可以看出数据已经处理好了。

3.6 偏差值处理

通过以下代码对偏差值进行处理，通过检查皮尔逊相关系数与缺失值比例是否大于20%，数值型特征方差是否小于0.1来进行删除特征操作。

# 去除低方差特征和与目标变量不相关的特征
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))

结果如下：

3.7 数据转换

通过直方图的形式查看连续性数据的分布，只有类似于正态分布的数据分布才是符合现实规律的数据，并且这有利于模型的训练。
通过以下代码查看直方图：

# 绘制直方图
col_name = df.columns
df[col_name].hist(bins=60,figsize=(16,12))

结果如下图所示：

既有符合正态分布的，也有比较不符合的，然后对不符合正态分布的特征数据进行数据转换：

# 针对不规则分布的变量进行非线性变换，一般进行log
log_col = ['Iron', 'Turbidity', 'Copper', 'Manganese']
show_col = []
for col in log_col:
    df[col + '_log'] = np.log(df[col])
    show_col.append(col + '_log')

df[show_col].hist(bins=60, figsize=(16, 12))

变换后的结果如下：

4. 模型训练

模型采用XGBoost，并使用随机网格搜索（RandomizedSearchCV）进行优化。

4.1 平衡数据及数据切分

由上面3.2的环状图图可知数据集比例有点失衡，所以在这里，我们使用imblearn库中的RandomUnderSampler()函数进行下采样，使得正负样本比例接近1:1。这里使用了随机数种子（random_state=21）来确保每次采样的结果相同。

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

from imblearn.under_sampling import RandomUnderSampler
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

df_processed = df.drop(['Target', 'Iron', 'Turbidity', 'Copper', 'Manganese', 'Color', 'Chlorine'], axis=1)
X = df_processed.values
y = df['Target'].values

# 下采样
under_sampler = RandomUnderSampler(random_state=21)
X_resampled, y_resampled = under_sampler.fit_resample(X, y)

X_train, X_test, y_train, y_test = train_test_split(X_resampled, y_resampled, 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))

输出训练集和测试集的形状：

4.2 网格搜索模型参数

# 定义超参数网格
param_grid = {
    'max_depth': [10, 15, 20],
    "gamma": [0, 1, 2],
    "subsample": [0.9, 1],
    "colsample_bytree": [0.3, 0.5, 1],
    'min_child_weight': [4, 6, 8],
    "n_estimators": [10,50, 80, 100],
    "alpha": [3, 4, 5]
}

其中网络参数包括

树的最大深度(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，∞]

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

4.3 模型训练

定义XGBoost分类器。

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

# 使用随机搜索（RandomizedSearchCV）来搜索XGBoost模型的最佳超参数组合，并基于F1分数（F1 score）来进行评估和选择最佳模型
# RandomizedSearchCV 是用于超参数搜索和模型选择的类。
# xgb 是 XGBoost 模型对象。
# param_grid 包含了要搜索的超参数范围。
# n_iter=10 指定了进行随机搜索的迭代次数。
# cv=3 指定了交叉验证的折数。
# refit=refit_score 指定了要使用的评分指标。
# scoring=sklearn.metrics.make_scorer(f1_score) 指定了评分函数，这里使用了 F1 分数。
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=sklearn.metrics.make_scorer(f1_score), verbose=10, return_train_score=True)
start = time.time()
rd_search.fit(X_train, y_train)
end = time.time()
print(rd_search.best_params_)
print(rd_search.best_score_)
print(rd_search.best_estimator_)
print(datetime.datetime.now() - start_time)
print('模型拟合时间:{:.2f}'.format(end - start))

通过搜索之后得到的相关最优参数以及模拟拟合时间如下：

1.未使用modin的模拟拟合时间：

2. 使用modin优化后的模拟拟合时间：

由上面两个图可以看出，未使用modin的模拟拟合时间是430.90s，而使用modin优化后的模拟拟合时间为246.37s，其速度整整提升了接近两倍！

然后再用搜索出来的最佳参数重新进行训练，并将模型转化为更高效的daal模型：

xgb = rd_search.best_estimator_
 
daal_model = d4p.get_gbt_model_from_xgboost(xgb.get_booster())

最后输出f1分数：

## Calculate predictions ##
daal_prob = d4p.gbt_classification_prediction(nClasses=2,
    resultsToEvaluate="computeClassLabels|computeClassProbabilities",
    fptype='float').compute(X_test, daal_model).probabilities # or .predictions
start1 = time.time()
xgb_pred = pd.Series(np.where(daal_prob[:,1]>.5, 1, 0), name='Target')
end1 = time.time()

xgb_auc = roc_auc_score(y_test, daal_prob[:,1])
xgb_f1 = f1_score(y_test, xgb_pred)
xgb_accuracy = accuracy_score(y_test, y_pred)
print('推理时间:{:.6f}'.format(end1 - start1))    
## Plot model results ##
print("\nTest F1 Scores: {:.2f}%\nAUC: {:.5f}\nAccuracy: {:.2f}".format(xgb_f1*100, xgb_auc, xgb_accuracy)) 
# plot_model_res(model_name='XGBoost', y_test=y_test, y_prob=daal_prob[:,1])

从结果来看，时间挺快的，F1分数也还可以。

4.4 结果可视化

当评估分类模型性能时，常用的指标包括 AUC（Area Under the ROC Curve）和 AUPRC（Area Under the Precision-Recall Curve）。这两个指标用于衡量模型在不同阈值下的性能表现。

绘制 AUC 曲线时，横轴为 FPR（False Positive Rate），纵轴为 TPR（True Positive Rate），即ROC曲线。AUC 衡量的是ROC曲线下的面积，取值范围在0到1之间，数值越接近1表示模型性能越好，0.5表示随机猜测。

绘制 AUPRC 曲线时，横轴为召回率（Recall），纵轴为精确率（Precision），即P-R曲线。AUPRC 衡量的是P-R曲线下的面积，同样取值范围在0到1之间，数值越高表示模型性能越好。

在绘制这两种曲线时，通常会通过模型预测值和真实标签计算出各个阈值下的 FPR、TPR、Precision 和 Recall 值，然后根据这些值绘制出 ROC 曲线和 P-R 曲线，最后计算出对应的 AUC 和 AUPRC 值来评估模型的性能优劣。

# 绘制 AUC 图像
fpr, tpr, thresholds = roc_curve(y_test, daal_prob[:,1])
roc_auc = auc(fpr, tpr)

plt.figure(figsize=(8, 6))
plt.plot(fpr, tpr, color='darkorange', lw=2, label='ROC curve (AUC = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver Operating Characteristic (ROC) Curve')
plt.legend(loc="lower right")
plt.grid(True)
plt.show()

# 绘制 AUPRC 图像
precision, recall, _ = sklearn.metrics.precision_recall_curve(y_test, daal_prob[:,1])
auprc = auc(recall, precision)

plt.figure(figsize=(8, 6))
plt.plot(recall, precision, color='darkorange', lw=2, label='AUPRC curve (AUPRC = %0.2f)' % auprc)
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.title('Area Under Precision-Recall Curve (AUPRC)')
plt.legend(loc="lower right")
plt.grid(True)
plt.show()

结果如下图所示：

5. 对给定测试集test_data进行预测

5.1 读取并查看测试集

df1 = pd.read_csv('test_data.csv')

print(df1.info(), df1.describe())

5.2 处理测试集数据

1. 数据类型转换

# 定义函数将分类变量映射为数值
def trans(l):
    data = {}
    for i in range(len(l)):
        data[l[i]] = i + 1
    return data

# 将分类变量转换为数值型
d1 = trans(df1['Source'].unique())
d2 = trans(df1['Month'].unique())
d3 = trans(df1['Color'].unique())
df1['Source'] = df1['Source'].map(d1)
df1['Month'] = df1['Month'].map(d2)
df1['Color'] = df1['Color'].map(d3)

# 输出转换后的数据信息
print(df1.info(), df1.describe())
# print(df.info)

2. 缺失值和重复值处理

display(df1.isna().sum())
missing = df1.isna().sum().sum()
duplicates = df1.duplicated().sum()
print("\n共{:,.0f} 缺失值 in the data.".format(missing))
print("共{:,.0f} 重复值 in the data.".format(duplicates))

# 填充缺失值
data = data.fillna(data.interpolate())
data = data.drop_duplicates()
print(data.isna().sum(),data.duplicated().sum())
print("\n共{:,.0f} 缺失值 in the data.".format(data.isna().sum().sum()))
# print("共{:,.0f} 重复值 in the data.".format(data.duplicated().sum()))

3. 数据集准备

# 特征工程和准备数据集
X = df.drop(['Target', 'Iron', 'Turbidity', 'Copper', 'Manganese', 'Color', 'Chlorine'], axis=1)
y = df['Target']

# 将整个数据集作为测试集
X_test = X
y_test = y

scaler = StandardScaler()
X_test_scaled = scaler.fit_transform(X_test)

print("Test Shape: {}".format(X_test_scaled.shape))

5.3 结果预测

测试集处理及预测是直接在上面训练集处理的后面直接写的，所以未保存模型，但是测试集预测是直接运用的上面训练集训练的模型。

daal_prob = d4p.gbt_classification_prediction(nClasses=2,
    resultsToEvaluate="computeClassLabels|computeClassProbabilities",
    fptype='float').compute(X_test, daal_model).probabilities # or .predictions
start1 = time.time()
xgb_pred = pd.Series(np.where(daal_prob[:,1] > 0.5, 1, 0), name='Target')
end1 = time.time()

xgb_auc = roc_auc_score(y_test, daal_prob[:,1])
xgb_f1 = f1_score(y_test, xgb_pred)
xgb_accuracy = accuracy_score(y_test, xgb_pred)  # 使用xgb_pred替换y_pred
print('推理时间:{:.6f}'.format(end1 - start1))    
## Plot model results ##
print("\nTest F1 Scores: {:.2f}%\nAUC: {:.5f}\nAccuracy: {:.2f}".format(xgb_f1*100, xgb_auc, xgb_accuracy)) 
# plot_model_res(model_name='XGBoost', y_test=y_test, y_prob=daal_prob[:,1])

6. 总结

通过该项目，实现了一个完整的分类模型开发流程，从数据准备、模型训练到模型评估和预测。借助XGBoost和daal4py等库，能够高效地构建和部署机器学习模型。通过可视化和指标评估，能够全面地了解模型的性能和预测效果。在处理新数据时，能够灵活地使用保存的模型进行预测，并得到所需的结果。

根据上面得到的结果也可以看出运用XGBoost和daal4py等库，能明显提高数据处理的速度，大大缩短训练时间。