初步观察完相关系数之后我们着手建模.笔者选择了决策树对数据进行分类.在建立决策树之前,为了防止过拟合,除了算法中自带的剪枝之外,还需要限制决策树的深度.结合数据集大小,我们对可能的决策树深度[1,16)进行网格搜索调参.使用sklearn包中的GridSearchCV进行网格搜索和交叉验证.使用网格搜索交叉验证出的最优决策树深度进行决策树的训练.

网格搜索结论是决策树深度为2的时候交叉验证准确率较高,为0.7030.选择决策树深度2建立决策树.在训练完模型之后from sklearn.tree import plot_tree进行决策树可视化,可视化结果如图:

多次运行下来决策树的右枝条会有变化.变化如图:

从树图中可以看出,右边枝条其实对风化情况的区分并不明显.附加的颜色条件最多只能区分1个样本.这一点也符合热力图分析.而玻璃类型对风化与否的区分较为明显.值得注意的是,在高钾玻璃中纹饰B全部风化.热力图中我们也无法拒绝纹饰B和风化之间的关系.

综上所述,我们认为风化情况主要跟玻璃类别有关.铅钡玻璃更容易风化,高钾玻璃中纹饰B更容易风化.风化与否与颜色并无较大关联,在除高钾玻璃之外,风化与否与纹饰也并没有较大关联.

第一小问(Python)所有代码:

import numpy as np
import pandas as pd
from sklearn.model_selection import GridSearchCV # 通过网格方式来搜索参数
from sklearn import tree
from matplotlib import pyplot as plt
import matplotlib
font = {#中文绘图
    'family':'SimHei',
    'weight':'bold',
    'size':12
}
matplotlib.rc("font", **font)
import seaborn as sns

data = pd.read_csv("表单1.csv")#导入数据集
features = ["文物编号","纹饰","类型","颜色"]
X = data[features]#设置待估Xy
y = data["表面风化"]
X.set_index("文物编号",inplace=True)#更改索引
features = ["纹饰","类型","颜色"]

plt.figure(figsize=[15, 10])#生成热力图,初步观察变量相关情况
sns.heatmap(pd.get_dummies(data[["纹饰","类型","颜色","表面风化"]]).corr(method='spearman'),cmap="Greens",mask=abs(pd.get_dummies(data[["纹饰","类型","颜色","表面风化"]]).corr(method='spearman').corr(method='spearman'))<0.279, linewidths=0.5, annot=True)
plt.show()

X = pd.get_dummies(X[features])#onehot

# 设置需要搜索的参数值，在这里寻找最优的决策树深度
parameters = {'max_depth':[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]}
model = tree.DecisionTreeClassifier()  # 注意：在这里不用指定参数
# GridSearchCV
clf = GridSearchCV(model, parameters, cv=5)   
clf.fit(X, y)
# 输出最好的参数以及对应的准确率
print ("best score is: %.4f"%clf.best_score_, "  best param: ",clf.best_params_,)

model = tree.DecisionTreeClassifier(max_depth=2)
model = model.fit(X, y)

fig = plt.figure(figsize=(10,10))
class_names = ["无风化","风化"]
from sklearn.tree import plot_tree # 树图
plot_tree(
    model, 
    feature_names = X.columns,  
    class_names= class_names,
    filled=True,
)

第二小问

(2)结合玻璃的类型, 分析文物样品表面有无风化化学成分含量的统计规律

在这里我仍然选择使用热力图先初步观察一下.过程同第一小问.

从图中可以发现明显的相关关系.正值代表正相关,负值代表负相关.绝对值越接近1,越接近线性相关.越接近0,关系越小.

由于该题询问的是统计规律,所以我们使用Excel进行描述统计分析.题目要求结合玻璃类型进行分析,于是我们把数据分为四组:高钾风化,高钾不风化,铅钡风化,铅钡不风化.在此笔者只对高钾风化玻璃文物的描述统计进行文字分析.其余的描述统计和图表参考文件:

<待上传的文件:四分组化学成分描述统计.xlsx>

上图展示了高钾风化玻璃数据集和其描述统计量.其中有整列为0的统计量是由化学成分缺失导致的.用到了一些不常用的统计量:偏度系数和峰度系数 .具体描述参考链接,我们直接进行描述统计:

从平均值角度,SiO2在玻璃中占比毫无疑问的最大,它是所有玻璃的主成分.别的所有微量成分的占比都在2%以下.

从峰度系数角度,二氧化硅,氧化钾,氧化镁,三氧化二铁四个组分的峰度系数为负,说明这四个组分的分布比正态分布更平缓,两侧具有较多的离散值.其余组分峰度系数为正,说明其分布比正态分布更陡峭,两侧”尾巴”比较瘦,分布更集中.

从偏度系数角度,二氧化硅,氧化钙,氧化镁,三氧化二铝,氧化铜成正偏态,指示该成分分布中有右长尾,可能存在离对应化学成分均值较远,较大的成分.其余成分分布偏左,体现出可能具有距离对应化学成分均值较远较小的化学成分.

然后对高钾玻璃的各个成分进行可视化.由相关系数可知,二氧化硅是随着风化而流失的.但风化后二氧化硅占比反而上升,说明微量成分更容易流失.

考虑到二氧化硅占比巨大,于是在微量成分分析之中去掉了二氧化硅.

对于微量元素分析发现风化前后有部分微量元素从图表中消失.他们是氧化钾,氧化镁.根据相关系数,它们可能是在风化过程中完全流失了.结合峰度系数和相关系数分析可以猜测风化进程使得随风化流失的成分分布变得极端.而氧化铜,三氧化二铝这样的成分可能不容易流失,尤其是三氧化二铝可能具有抗风化作用.所以占比变大.

在做第三小题的时候,顺便使用第三小题的岭回归模型去试探性探索了风化程度,玻璃类型和化学成分之间的定量关系.在线性模型中,负号指示对自变量有负贡献,正号指示对自变量有正贡献,系数矩阵如下:

1.高钾,铅钡和风化程度是因变量,各个化学成分占比是自变量.

回归截距
(高钾)rich_k 0.6396314
(铅钡)PbBa 0.3603686
(风化)weathering 0.4576642

高钾铅钡风化程度
SiO2 0.04857851 -0.04857851 -0.02583662
Na2O -3.66590529 3.66590529 -5.57513205
K2O 3.97230658 -3.97230658 -1.60494187
CaO 0.00000000 0.00000000 0.00000000
MgO -0.73045301 0.73045301 -1.92845812
Al2O3 -1.91151436 1.91151436 -1.82996762
Fe2O3 -0.10065749 0.10065749 -6.48401062
CuO 1.18870830 -1.18870830 1.11851119
PbO -1.00796999 1.00796999 0.75990028
BaO -1.51898888 1.51898888 -0.88100097
P2O5 -0.27666001 0.27666001 4.06530995
SrO 0.00000000 0.00000000 0.00000000
SnO2 0.00000000 0.00000000 0.00000000
SO2 0.74433109 -0.74433109 6.83425779

从该矩阵可以看出线性模型如何分出高钾和铅钡两类玻璃,并能够表现出各个化学成分对风化程度的贡献.该矩阵与相关性系数热力图互相验证.

2.高钾玻璃各个成分作为自变量,风化程度作为因变量:

回归截距
0.3591105

风化程度
SiO2 0.5155128
Na2O -1.7055526
K2O -1.5235123
CaO -1.6537503
MgO -5.6838663
Al2O3 -1.8538162
Fe2O3 -1.0543517
CuO -0.5243118
PbO -6.3371961
BaO -3.1674837
P2O5 -2.3377587
SrO -24.2124372
SnO2 -8.6608621
SO2 -18.8954178

3.铅钡玻璃各个成分作为自变量,风化程度作为因变量:

回归截距
0.5840007

风化程度
SiO2 -1.3158729
Na2O -1.8380630
K2O -16.6983000
CaO -6.9069850
MgO 5.1302300
Al2O3 2.0967995
Fe2O3 -4.1001075
CuO 0.3635241
PbO 1.4328986
BaO -0.2036333
P2O5 4.1038535
SrO -14.5897249
SnO2 36.3990629
SO2 6.4818275

从2,3两个系数矩阵可以发现,铅钡玻璃和高钾玻璃之间有明显不同的风化进程.

用到的R代码如下,需要注意的是笔者在这里使用了很多共用的变量,所以这三个模型需要依次分块执行,不能一次全部运行,否则得到的结果将是最后一个模型的结果.

library(foreign)
library(glmnet)
class_weathering_chemical <- read.csv("./class_weathering_chemical.csv")
#####完全多变量回归#####
x <- as.matrix(class_weathering_chemical[,c(-1,-16:-18)])#定义回归自变量因变量
y <- as.matrix(class_weathering_chemical[,c(16:18)])
fit <- glmnet(x,y,"mgaussian",nlambda = 1000,alpha=1)#为了查看变量影响因素而进行的回归
plot(fit, xvar="lambda", label=TRUE)
lasso_fit <- cv.glmnet(x,y,family="mgaussian",alpha=1,type.measure = "mse",nlambda=1000)#交叉验证,最佳lambda
plot(lasso_fit)
lasso_best <- glmnet(x,y,family="mgaussian",alpha = 1,lambda = lasso_fit$lambda.min)#用λmin建立预测模型

lasso_best$a0#回归截距

k <- as.matrix(lasso_best$beta$rich_k)
k <- cbind(k,as.numeric(lasso_best$beta$PbBa))
k <- cbind(k,as.numeric(lasso_best$beta$weathering))
k <- as.data.frame(k)
names(k) <- c("高钾","铅钡","风化程度")#k为回归系数矩阵
print(k)


#####高钾玻璃的风化情况#####
richK_weathering_chemical <- read.csv("./richK_weathering_chemical.csv")
x <- as.matrix(richK_weathering_chemical[,c(-1,-16)])#定义回归自变量因变量
y <- as.matrix(richK_weathering_chemical[,16])#风化程度
fit <- glmnet(x,y,"gaussian",nlambda = 1000,alpha=0)#为了查看变量影响因素而进行的回归
plot(fit, xvar="lambda", label=TRUE)
lasso_fit <- cv.glmnet(x,y,family="gaussian",alpha=0,type.measure = "mse",nlambda=1000)#交叉验证,最佳lambda
plot(lasso_fit)
lasso_best <- glmnet(x,y,family="gaussian",alpha = 0,lambda = lasso_fit$lambda.min)#用λmin建立预测模型

lasso_best$a0#回归截距

k <- as.data.frame(as.matrix(lasso_best$beta))
names(k) <- "风化程度(0-2)"
print(k)

#####铅钡玻璃的风化情况#####
PbBa_weathering_chemical <- read.csv("./PbBa_weathering_chemical.csv")
x <- as.matrix(PbBa_weathering_chemical[,c(-1,-16)])#定义回归自变量因变量
y <- as.matrix(PbBa_weathering_chemical[,16])#风化程度
fit <- glmnet(x,y,"gaussian",nlambda = 1000,alpha=0)#为了查看变量影响因素而进行的回归
plot(fit, xvar="lambda", label=TRUE)
lasso_fit <- cv.glmnet(x,y,family="gaussian",alpha=0,type.measure = "mse",nlambda=1000)#交叉验证,最佳lambda
plot(lasso_fit)
lasso_best <- glmnet(x,y,family="gaussian",alpha = 0,lambda = lasso_fit$lambda.min)#用λmin建立预测模型

lasso_best$a0#回归截距

k <- as.data.frame(as.matrix(lasso_best$beta))
names(k) <- "风化程度(0-2)"
print(k)

第三小问

(3)根据风化点检测数据, 预测其风化前的化学成分含量

根据待预测对象的性质,并参考下面这张著名的sklearn cheat-sheet,笔者选择了岭回归模型对化学成分进行预测.

本题我们结合表单1和表单2之中的信息,以纹饰,风化程度,颜色,作为自变量,化学成分作为因变量

结合玻璃的类型,我们对大集合进行了划分.在预测风化前成分的时候,分为高钾和铅钡两类分别进行预测.我们导入数据,里面有风化和无风化的数据.提取所有的数据作为训练集;提取出所有的风化后的文物数据作为待预测的集合,并把预测集合的所有的风化程度设置为0.

铅钡玻璃风化前成分预测结果:

笔者在这里遇到了一些仍没有解决的问题:就是虽然预测之后的成分结果仍然是定和为1的,但是出现了负成分的数字.虽然可以直接将负成分去掉并把所有的成分重新定和,但我并不确定这样的做法靠不靠谱.别的预测结果可以通过运行代码直接执行出来,目前我还没有专门去保存.

第三小问R语言代码如下:

library(foreign)
library(glmnet)
##########高钾玻璃的风化前预测##########
#####数据准备#####
rich_K_weathering_chemical <- read.csv("./predict_b4_weathering.csv")
rich_K_train <- data.frame()#划分训练集与测试集.有风化的做测试集,并把风化位置为0后作为预测集
for (i in c(1:nrow(rich_K_weathering_chemical))) {#从大集合中取出高钾玻璃作为训练集
  if(rich_K_weathering_chemical[i,"rich_k"] == 1){
    rich_K_train <- rbind(rich_K_train,rich_K_weathering_chemical[i,])
  }
}
rich_K_predict <- data.frame()
for (i in c(1:nrow(rich_K_train))) {#从训练集中取出有风化的作为测试集,并把风化程度置为0
  if(rich_K_train[i,"weathering"] != 0){
    rich_K_predict <- rbind(rich_K_predict,rich_K_train[i,])
  }
}
rich_K_predict[,"weathering"] <- 0#并把风化程度置为0
train_y <- as.matrix(rich_K_train[,2:15])#定义回归训练自变量因变量
train_x <- as.matrix(rich_K_train[,16:30])#自变量
#####模型拟合#####
fit <- glmnet(train_x,train_y,"mgaussian",nlambda = 1000,alpha=0)#为了查看变量影响因素而进行的回归
plot(fit, xvar="lambda", label=TRUE)
lasso_fit <- cv.glmnet(train_x,train_y,family="mgaussian",alpha=0,type.measure = "mse",nlambda=1000)#交叉验证,最佳lambda
plot(lasso_fit)
lasso_best <- glmnet(train_x,train_y,family="mgaussian",alpha = 0,lambda = lasso_fit$lambda.min)#用λmin建立预测模型
#####模型预测#####
#定义回归预测自变量
predict_x <- as.matrix(rich_K_predict[,16:30])#自变量
predict_result <- as.data.frame(predict(lasso_best,predict_x))#进行预测
predict_result <- cbind(rich_K_predict[,1],predict_result)#优化预测表格格式
names(predict_result) <- names(rich_K_predict[,1:15])
#print(predict_result)



##########铅钡玻璃的风化前预测##########
#####数据准备#####
PbBa_weathering_chemical <- read.csv("./predict_b4_weathering.csv")
PbBa_train <- data.frame()#划分训练集与测试集.有风化的做测试集,并把风化位置为0后作为预测集
for (i in c(1:nrow(PbBa_weathering_chemical))) {#从大集合中取出铅钡玻璃作为训练集
  if(PbBa_weathering_chemical[i,"PbBa"] == 1){
    PbBa_train <- rbind(PbBa_train,PbBa_weathering_chemical[i,])
  }
}
PbBa_predict <- data.frame()
for (i in c(1:nrow(PbBa_train))) {#从训练集中取出有风化的作为测试集,并把风化程度置为0
  if(PbBa_train[i,"weathering"] != 0){
    PbBa_predict <- rbind(PbBa_predict,PbBa_train[i,])
  }
}
PbBa_predict[,"weathering"] <- 0#并把风化程度置为0
train_y <- as.matrix(PbBa_train[,2:15])#定义回归训练自变量因变量
train_x <- as.matrix(PbBa_train[,16:30])#自变量
#####模型拟合#####
fit <- glmnet(train_x,train_y,"mgaussian",nlambda = 1000,alpha=0)#为了查看变量影响因素而进行的回归
plot(fit, xvar="lambda", label=TRUE)
lasso_fit <- cv.glmnet(train_x,train_y,family="mgaussian",alpha=0,type.measure = "mse",nlambda=1000)#交叉验证,最佳lambda
plot(lasso_fit)
lasso_best <- glmnet(train_x,train_y,family="mgaussian",alpha = 0,lambda = lasso_fit$lambda.min)#用λmin建立预测模型
#####模型预测#####
#定义回归预测自变量
predict_x <- as.matrix(PbBa_predict[,16:30])#自变量
predict_result <- as.data.frame(predict(lasso_best,predict_x))#进行预测
predict_result <- cbind(PbBa_predict[,1],predict_result)#优化预测表格格式
names(predict_result) <- names(PbBa_predict[,1:15])
#print(predict_result)

第二题

笔者参考了他人第二大题的思路进行建模.第二大题被我分成如下三问:

(1)依据附件数据分析高钾玻璃、铅钡玻璃的分类规律;

(2)对于每个类别选择合适的化学成分对其进行亚类划分, 给出具体的划分方法及划分结果;

(3)并对分类结果的合理性和进行分析;

在第二题,由于化学成分数据维数很高,无法直观的理解聚类之后的关系,直接聚类也可能出现维数灾难.故我首先使用PCA进行降维,随后使用Kmeans进行聚类,为集合打上亚类标签.但由于降维后的主成分不具有良好的可解释性,所以另外利用化学成分,使用决策树对分类标签进行监督学习,从而得到可视化的易于理解的分类标准.

第一小问

(1)依据附件数据分析高钾玻璃、铅钡玻璃的分类规律

这是一个简单的二分类问题.在这里利用Python使用决策树算法进行分类.与之前的步骤基本一样.这里使用了一份同时含有化学成分,玻璃类别标签和风化程度的数据集来进行初步的二分类.

(文末链接文件中的class_weathering_chemical.csv)

以化学成分为自变量,玻璃类型标签为因变量构建决策树.随后进行决策树的可视化如图.

这是经过网格搜索之后,最优深度为1的决策树,交叉验证时的准确率能够高达100%,分类结果的准确率也达到了100%,各个节点中信息熵为0.也就是说,只通过氧化铅的含量即可判断该玻璃是高钾还是铅钡.如果氧化铅含量少于0.061,我们则认为它是高钾玻璃.如果高于0.061,我们认为它是铅钡玻璃.

第一小问(Python)代码如下:

import numpy as np
import pandas as pd
from sklearn.model_selection import GridSearchCV # 通过网格方式来搜索参数
from sklearn.tree import DecisionTreeClassifier as dtc
import matplotlib.pyplot as plt # 可视化
from matplotlib import rcParams # 图大小
from sklearn.tree import plot_tree # 树图

data = pd.read_csv("class_weathering_chemical.csv")#导入数据集
features = data.columns[1:15]
X = data[features]#设置待估Xy
y = data[data.columns[15:17]]

# 设置需要搜索的参数值，在这里寻找最优的决策树深度
parameters = {'max_depth':[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]}
model = dtc()  # 注意：在这里不用指定参数
# GridSearchCV
clf = GridSearchCV(model, parameters, cv=5)   
clf.fit(X, y)
# 输出最好的参数以及对应的准确率
print ("best score is: %.4f"%clf.best_score_, "  best param: ",clf.best_params_,)

model = dtc(max_depth=3,criterion = 'entropy')#一层无法可视化,于是构建两层
model = model.fit(X, y)

rcParams['figure.figsize'] = (5, 5)
target_names = ['rich_k','PbBa']
from sklearn.tree import plot_tree # 树图
plot_tree(
    model, 
    feature_names = features,  
    class_names= target_names,
    filled=True,
    rounded = True
)
plt.savefig('如何划分铅钡和高钾.png')

第二小问

(2)对于每个类别选择合适的化学成分对其进行亚类划分, 给出具体的划分方法及划分结果,

在这里,我们分别针对铅钡风化,铅钡无风化,高钾风化,高钾无风化进行亚类划分.

为什么不直接基于铅钡和高钾两类进行亚类划分,而是要先区分风化与不风化呢?原因在于由前文的分析可知,风化和无风化之间本来就存在化学成分差异.按照题目要求,选择合适的化学成分对其进行亚类划分的时候很可能会让模型再次”发现”一次风化和无风化玻璃之间的化学成分差异,导致亚类划分中做了一遍无用功.

我们把玻璃先划分为四个数据集,分别对应四类不同的玻璃:

然后对这四个数据集进行同样的操作.即可划分出四个类别里的亚类.笔者在此只对铅钡风化玻璃的亚类划分进行详细的解释,其它三类都是同理可得.

首先我们进行PCA降维,将十五个化学成分集成为两个主成分,把这两列加入到导入的data数据框中,并观察这两个主成分贡献了百分之多少的方差:

两列主成分贡献了接近85%的方差.我们认为这两列主成分可以作为原数据集的主成分.随后我们利用这两列主成分fac1,fac2作为聚类依据,对数据集进行KMeans聚类.考虑到数据集规模不大,亚类的种类不应过多,我们将再[2,3,4,5]中选取K的值.

for n, ax  in zip([2,3,4,5], sub.flatten()):
    kmean = KMeans(n_clusters=n)
    data.loc[:,str(n) + 'type'] = kmean.fit_predict(data.loc[:,['fac1','fac2']])

我们利用这两行代码分别建立K=2,3,4,5时的KMeans模型,随后分别用不同聚类数的模型对数据集进行聚类,并将聚类的标签加入原数据集data的新列中.执行完这些代码之后data数据集中除了化学成分列,主成分列fac1,fac2,还被加上了不同聚簇数量下某个样本的聚簇归属(最后四列).如下:

注:unnamed列因不明原因被引入,原数据集里并没有这一列.我有点懒没去删,

但这一列不会影响到程序运行与数据分析.

然后我们对这些聚类结果进行一个简单的可视化:

但是主成分分析之后的两列主成分和原成分之间不具有良好的可解释性,我们无法依据主成分给出良好的用于分亚类的合适化学成分.所以我们这边再次以化学成分作为自变量,亚类标签作为因变量,使用决策树去探讨亚类的具体分类标准.但是,要怎么知道在2,3,4,5个聚类中,几个聚类数最恰当呢?因为数据集量并不大,我们选择全部试一遍,从2分类到5分类,对不同的K利用决策树模型进行分类,用决策树模型交叉验证准确率寻找最优分类数量,使用准确率最高的那个聚类数作为亚类的数量.对于铅钡风化的玻璃,我们得到2分类的交叉验证准确率最高为1.随后便使用2分类的那行标签作为因变量,化学成分作为自变量去做决策树,并进行可视化得到如下结果:

分类器得出,铅钡风化玻璃以二氧化硫含量作为亚类划分标准.小于等于0.01的分为1类,高于0.01的分为0类.剩下的几类以此类推,可以得到亚类的划分结果.亚类划分结果总结:

分类标准总结
玻璃大类	亚类区分标志化学成分	亚类
铅钡风化	二氧化硫	(SO2<=0.01) class1
		(else) class0
铅钡无风化	二氧化硅	(SiO2>0.581) class0
		(0.398<=SiO2<=0.581) class2
		(SiO2<=0.398) class1
高钾风化	三氧化二铝	(Al2O3<=0.017) class0
		(else) class1
高钾无风化	二氧化硅	(SiO2<=0.744) class1
		(else) class0

(简单的发现一个规律:风化后都以微量元素作为分类标准.风化前都以二氧化硅作为划分标准.)

第三小问

(3)并对分类结果的合理性和敏感性进行分析.

合理性分析:主成分分析后再进行聚类,在避免维数灾难的同时保证了每种化学成分都具有恰当的贡献率.对不同的聚类数量利用分类模型进行交叉验证检验,保障了分类的准确性以及合理性.

敏感性分析:模型极为简单明了,只有一个关键化学成分和一到两条分类标准.只要超过阈值就进行分类.改变非关键变量对该条件下的分类不会产生影响.只有关键变量的波动会对分类会产生影响

由于第二小问四段代码都是同质的,在此便不展示全部的代码.只展示上面分析用到的代码.所有的代码可以在文末的附件中取用.

第二题第二小问的代码:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

data= pd.read_csv("./亚类分析/PbBa_weathering_chemical.csv")#铅钡风化的亚类分析

from sklearn.decomposition import PCA #使用PCA进行降维
X = data.iloc[:,1:15]
pca = PCA(n_components=2, copy=True,)
newX = pca.fit_transform(X)
#invX = pca.inverse_transform(newX)
#print(pca.components_)
print(pca.explained_variance_ratio_)
data['fac1'] = newX[:,0] #两个因子分别命名和储存到data中
data['fac2'] = newX[:,1]

from sklearn.cluster import KMeans #KMeans

fig, sub = plt.subplots(2,2,dpi = 300, figsize = (10,9))#将无监督聚类的散点进行可视化
plt.subplots_adjust(wspace=0.2, hspace=0.2)
for n, ax  in zip([2,3,4,5], sub.flatten()):
    kmean = KMeans(n_clusters=n)
    data.loc[:,str(n) + 'type'] = kmean.fit_predict(data.loc[:,['fac1','fac2']])
    sns.scatterplot(data = data, x = 'fac1', y = 'fac2', c =data.loc[:,str(n) + 'type'], ax =ax )
    ax.set_title(str(n)+'clusters distribution')
    plt.tight_layout()
    plt.savefig('铅钡风化聚类结果展示.png', dpi = 300, facecolor = 'w', bbox_inches = 'tight')

import numpy as np
from sklearn.model_selection import GridSearchCV # 通过网格方式来搜索参数
from sklearn.tree import DecisionTreeClassifier as dtc
import matplotlib.pyplot as plt # 可视化
from matplotlib import rcParams # 图大小
from sklearn.tree import plot_tree # 树图

features = data.columns[1:15]
X = data[features]#设置待估Xy
for i in range(2,6):#从二分类到5分类,对不同的y利用决策树进行网格搜索,寻找最优分类数量
    y = data.loc[:,str(i)+'type']
    # 设置需要搜索的参数值，在这里寻找最优的决策树深度
    parameters = {'max_depth':[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]}
    model = dtc()  
    clf = GridSearchCV(model, parameters, cv=5) #网格搜索  
    clf.fit(X, y)
    # 输出最好的参数以及对应的准确率
    #print ("在",i,"个聚类中","最优分数: %.4f"%clf.best_score_, "  最佳参数: ",clf.best_params_,)
    if(i == 2):#存储最好的各类参数
        current_best = clf.best_score_
        best_number = i
        depth = clf.best_params_['max_depth']
    else:
        if(current_best < clf.best_score_):
            best_number = i
            current_best = clf.best_score_
            depth = clf.best_params_['max_depth']
print("在",best_number,"个聚簇时决策树交叉验证准确率最高,为",current_best,"决策树最大深度",depth)
#深度为1的决策树无法画图.且希望多展示一层.
depth = depth + 1

classify_result = data[data.columns[0:15]]#输出最优分类结果
classify_result['sub_class'] = data[str(best_number)+'type']
classify_result.to_csv('铅钡风化亚类分类结果.csv')

#多次执行以上语句之后,选择最优分数对应的聚类数量(如相同分数取聚类数少的:奥卡姆剃刀)和最优决策树层数+1作为铅钡风化玻璃划分子类数量
y = data.loc[:,str(best_number)+'type']#利用上面的最优参数取对应的列作为y
class_names = []
for i in range(0,best_number):
    class_names.append(str(i))#为可视化创建类标签
model = dtc(max_depth=depth,criterion = 'entropy')
model = model.fit(X, y)
rcParams['figure.figsize'] = (5, 5)
from sklearn.tree import plot_tree # 树图
plot_tree(
    model, 
    feature_names = features,
    class_names = class_names,  
    filled=True,
)
plt.savefig('铅钡风化玻璃中的亚类划分.png')

第三题

(1)对附件表单3中未知类别玻璃文物的化学成分进行分析,鉴别其所属类型, 并对分类结果的敏感性进行分析.

该数据集给出了A1-A8八个玻璃文物,信息包括表面风化情况和各个需要我们判断其所属的类型.结合前面两大题的结果,我在划分的时候加入了对亚类的划分.首先对其进行大类的划分.根据下面这个决策树的条件,我们判断A1,A6,A7是高钾玻璃,其余都是铅钡玻璃.

之后,结合风化情况,我们对这几类玻璃进行亚类划分.划分的标准在亚类划分中已经展示过,在这里就不再展示了.最后的划分结果如下表所示:

文物编号	表面风化	大类	亚类
A1	无风化	高钾	高钾_无风化_class0(高二氧化硅)
A2	风化	铅钡	铅钡_风化_class1(低二氧化硫)
A3	无风化	铅钡	铅钡_无风化_class1(低二氧化硅)
A4	无风化	铅钡	铅钡_无风化_class1(低二氧化硅)
A5	风化	铅钡	铅钡_风化_class1(低二氧化硫)
A6	风化	高钾	高钾_风化_class0(低三氧化二铝)
A7	风化	高钾	高钾_风化_class1(高三氧化二铝)
A8	无风化	铅钡	铅钡_无风化_class2(中二氧化硅)

敏感性分析:由于这些分类都由决策树指定,且决策树层数不超过两层,所以这样的分类得以免受无关变量的干扰,只与关键的指标有关,只有关键化学成分含量的变化才会影响分类结果.

第四题

(1)针对不同类别的玻璃文物样品, 分析其化学成分之间的关联关系,

(2)并比较不同类别之间的化学成分关联关系的差异性.

第一小问

(1)针对不同类别的玻璃文物样品, 分析其化学成分之间的关联关系,

为了统计观察方便,我只对铅钡风化,铅钡无风化,高钾风化,高钾无风化四个大类进行了化学成分关联性的分析.分别针对这四类进行了相关性热力图的绘制,并只保留置信度95%以上的色块,并且分别对四幅图进行分析.

铅钡无风化

最显著的特点是铅钡和二氧化硅的负相关性.由资料可得铅钡玻璃中二氧化硅的成分占比较低.主要由氧化铅组成.可以合理推测添加氧化铅和氧化钡(使用助融剂)可导致二氧化硅水平的下降.

铅钡风化

相比于不风化的时候,多出了很多化学成分并出现了很多相关关系.最显著的特征是钙镁铝铁(除去铝和钙)这四大元素之间的正相关性相比于未风化时显著增强.由之前的分析,风化过程中铅钡玻璃的钡元素流失明显,可以推测镁,铁等元素是风化过程中的环境元素交换到玻璃中的结果.

高钾不风化

高钾玻璃最显著的特征是对角线附近的元素相关性.高钾玻璃中的钠钾钙之间的强正相关考虑是加入草木灰这样的助融剂而导致的.但助融剂会导致二氧化硅占比的下降,这样的规律在以铅矿石作为玻璃助融剂的时候同样存在.

高钾风化

高钾玻璃风化之后与环境交换的元素相对于铅钡玻璃来说较少,也能够体现出高钾玻璃的抗风化能力强于铅钡玻璃.由于高钾玻璃风化后的样本比较少,且检测不到的元素类型很多,所以热力图中出现了较多的空缺.氧化钙,三氧化二铝都与二氧化硅含量成负相关,而铝和钙之间成正相关.

该小问代码原理与之前展示过的相同,便不再展示.需要获取请点击文末链接.

第二小问

(2)并比较不同类别之间的化学成分关联关系的差异性.

高钾玻璃和铅钡玻璃的横向对比:主要区别在于,相比于高钾玻璃,铅钡玻璃含有的元素种类更多,元素之间的相关性也容易产生了.而且由图像可知,这两者之间作比较,横向的元素关联不大,相关系数的热力图甚至呈互补的趋势.而主要的联系在于,不论使用何种助融剂,都会降低二氧化硅的占比.

风化和未风化之间的横向对比:风化会导致原来的主要成分和特征元素含量下降(比如二氧化硅和铅钡玻璃中的钡元素),并且会稀释原有元素之间的关联性,使得风化前后的热力图呈互补趋势,这说明不同风化流失元素的流失速率并不相同.从外界交换进来的金属元素之间会产生新的关联.

高钾风化和铅钡风化之间的横向对比:在第二题中我们已经得出了这两者具有完全不同的风化进程的结论.在这里的分析进一步加强了这一个结论的可信度.铅钡玻璃更容易风化,更容易和外界交换元素.而高钾玻璃不容易风化,甚至很少引进新的元素,只是自身原来具有的元素之间关联度下降.

四:不足之处与待改进分析

1.过程中大量使用了相关系数热力图去衡量变量之间的关系.但是有时候相关的变量很多,往往没有分析出关键点.并且由于背景化学知识的缺乏,也无法从关系之中得出有效的结论.

2.分类模型大量采用决策树模型.虽然简单有效,但是仍需要质疑如此简单的分类标准是否合适.

3.预测风化前化学成分的时候表单中出现了负成分.是否有更好的方法和模型可以处理成分数据?直接去掉负值并重新定和是否可以接受还有待讨论.

4.PCA+KMeans+决策树的方法中,衡量模型性能的指标是交叉验证准确率.但随着聚类的增多,交叉验证确实准确率有可能下降,但这个时候并不一定代表聚类不合理.是否存在更合理的聚类方式还有待探讨.

5.文中的很多分析,包括描述统计和相关关系分析中还需要更多的数据支持.

五:代码&数据集资源

百度网盘资源链接:源代码&数据集

提取码: m5kb

GitHub连接:源代码&数据集

代码&数据集末次更新日期:2022年12月8日

文章出处登录后可见！

已经登录？立即刷新

数据科学导论大作业:数模国赛C题(古代玻璃文物分类)

一、问题重述

2022年高教社杯全国大学生数学建模竞赛C题

【1】项目背景

【2】问题重述: 古代玻璃制品的成分分析与鉴别

二:简述

三:过程记录

1.数据预处理

2.四个小题的完成笔记

第一题

第二题

第三题

第四题

四:不足之处与待改进分析

五:代码&数据集资源

相关推荐