题目链接
http://www.mathorcup.org/detail/2260
本文仅演示数据预处理环节。
理论基础:https://www.cnblogs.com/fighterkaka22/p/14052346.html
数据预处理
本文取每个水池中,A、B两个采样点各理化因子的实测值的均值作为各理化因子的计算值。总磷、总淡、氨氮15周的数据可以参考附件一。而附件二中COD、溶氧、PH值间隔两周采集一次,与附件一数据不对称,不足以建立合理的模型,因此考虑利用现有数据插值以补充数据。
插值方法选用三次B样条插值,该方法可以很好的保持数据光滑性和连续性,减少信息量的损失。最终得到的数据如下。
源代码:(以轮虫-时间为例)
#暂时手动导入数据,下周开始学习自动从EXCEL表读取数据
#轮虫-时间
import numpy as np
import matplotlib.pyplot as plt
#进行样条插值
import scipy.interpolate as spi
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
#数据准备
X=[1,3,5,7,9,11,13,15]
Y1=[1913,1945,1920,2205,2260,2302,2385,2420]
Y2=[591,670,760,837,860,880,885,912]
Y3=[4.46,3.94,2.61,2.66,2.12,1.02,1.35,1.12]
Y4=[8.49,8.45,7.92,7.13,3.49,3.12,2.38,1.56]
#定义整周数点
xpoint=np.arange(1,15.1,1)
#定义曲线x点
xline=np.arange(1,15.1,0.1)
#进行三次样条拟合点
xpoint_r=spi.splrep(X,Y1,k=3) #样本点导入,生成参数
ypoint1=spi.splev(xpoint,xpoint_r) #根据观测点和样条参数,生成插值
xpoint_r=spi.splrep(X,Y2,k=3) #样本点导入,生成参数
ypoint2=spi.splev(xpoint,xpoint_r) #根据观测点和样条参数,生成插值
xpoint_r=spi.splrep(X,Y3,k=3) #样本点导入,生成参数
ypoint3=spi.splev(xpoint,xpoint_r) #根据观测点和样条参数,生成插值
xpoint_r=spi.splrep(X,Y4,k=3) #样本点导入,生成参数
ypoint4=spi.splev(xpoint,xpoint_r) #根据观测点和样条参数,生成插值
#进行三次样条拟合曲线
xline_r=spi.splrep(X,Y1,k=3) #样本点导入,生成参数
yline1=spi.splev(xline,xline_r) #根据观测点和样条参数,生成插值
xline_r=spi.splrep(X,Y2,k=3) #样本点导入,生成参数
yline2=spi.splev(xline,xline_r) #根据观测点和样条参数,生成插值
xline_r=spi.splrep(X,Y3,k=3) #样本点导入,生成参数
yline3=spi.splev(xline,xline_r) #根据观测点和样条参数,生成插值
xline_r=spi.splrep(X,Y4,k=3) #样本点导入,生成参数
yline4=spi.splev(xline,xline_r) #根据观测点和样条参数,生成插值
#输出数据表格
plt.figure(figsize=(15,10))
cell_text=[np.round(ypoint1,2),np.round(ypoint2,2),np.round(ypoint3,2),np.round(ypoint4,2)]
the_table=plt.table(cellText=cell_text,
rowLabels=["1号池","2号池","3号池","4号池"],
colLabels=xpoint,
)
#作图
plt.scatter(xpoint,ypoint1)
plt.plot(xline,yline1,label='1号池')
plt.scatter(xpoint,ypoint2)
plt.plot(xline,yline2,label='2号池')
plt.scatter(xpoint,ypoint3)
plt.plot(xline,yline3,label='3号池')
plt.scatter(xpoint,ypoint4)
plt.plot(xline,yline4,label='4号池')
plt.xticks([])
#plt.xlabel('时间(周)')
plt.ylabel('轮虫(10^6/L)')
plt.title('轮虫-时间(周)曲线图')
plt.legend()#显示标签
plt.show()