題目連結

http://www.mathorcup.org/detail/2260
本文僅演示資料預處理環節

資料預處理

本文取每個水池中，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()