目錄

一、寫在前面

二、數(shù)據(jù)

三、代碼-多元線性回歸

3.1 導入庫

3.2 導入數(shù)據(jù)

3.3 多元線性回歸模型

3.3.1 多元線性回歸-OLS

3.3.2 多元線性回歸模型預測值相對誤差

3.3.3 殘差圖

3.3.4 預測值與真實值分布圖

四、代碼-響應曲面

4.1 高斯消元

4.2 響應曲面繪制代碼

4.3 調用matching_3D繪制響應曲面

---

---

一、寫在前面

多項式回歸結合響應面分析方法，可以利用響應面圖直觀反應復雜的三維關系，從而清晰呈現(xiàn)兩個自變量和一個因變量之間關系的技術方法。

二、數(shù)據(jù)

除25組數(shù)據(jù)外，還加一組自然情況下的基準數(shù)據(jù)?：

No.?	?Temperature	pH?	Fe2+	Cu2+	Fe3+?	Y
0?	30	6.5	0?	0	0	7.55

三、代碼-多元線性回歸

3.1 導入庫

from mpl_toolkits.mplot3d.axes3d import Axes3D from matplotlib import cm from pylab import * from numpy import * import matplotlib.pyplot as plt import pandas as pd import math import numpy as np import copy plt.rcParams['axes.unicode_minus']=False #用于解決不能顯示負號的問題 mpl.rcParams['font.sans-serif'] = ['SimHei']

3.2 導入數(shù)據(jù)

#26 xArr = [[1,30,6.5,0,0,0],[1,30,1.5,0,0,0],[1,30,2,1,3,1],[1,30,2.25,3,5,3],[1,30,2.27,5,8,5],[1,30,2.41,8,10,8],[1,35,1.5,1,5,5],[1,35,2,3,8,8],[1,35,2.5,5,10,0],[1,35,3,8,0,1],[1,35,2.61,0,3,3],[1,40,1.5,3,10,1],[1,40,2,5,0,3],[1,40,2.5,8,3,5],[1,40,2.28,0,5,8],[1,40,3.23,1,8,0],[1,45,1.5,5,3,8],[1,45,2,8,5,0],[1,45,2.5,0,8,1],[1,45,2.44,1,10,3],[1,45,2.26,3,0,5],[1,50,1.5,8,8,3],[1,50,2,0,10,5],[1,50,2.17,1,0,8],[1,50,3,3,3,0],[1,50,2.82,5,5,1] ]#26 yArr = [7.55,7.14,7.2,7.05,6.82,6.51,6.73,6.69,6.46,6.75,6.70,6.55,6.3,6.21,6.18,5.97,5.95,5.9,5.5,5.72,5.6,5.62,5.29,5.57,5.30,5.21]# print(len(xArr),len(yArr))

3.3 多元線性回歸模型

3.3.1 多元線性回歸-OLS

這里插一點和題目不相關的東西，像這種有許多X與對應的Y的數(shù)據(jù)，可以考慮進行多元線性回歸

回歸代碼如下：（最小二乘法）

#最小二乘法 OLS def standRegres(xArr,yArr):xMat = mat(xArr)yMat = mat(yArr).TxTx = xMat.T*xMatif linalg.det(xTx) == 0.0:print("This matrix is singular, cannot do inverse")returnws = xTx.I * (xMat.T*yMat) # print(ws)return ws

3.3.2 多元線性回歸模型預測值相對誤差

采用多元線性回歸方程對驗證每組數(shù)據(jù)相對誤差：

mySum = 0 sse = 0 yPerList = []#采用全部數(shù)據(jù)進行訓練 ws = standRegres(xArr,yArr) #ws即為方程系數(shù) print(ws) for index,x in enumerate(xArr):yPer = float(x*ws) #yPer即為預測值yPerList.append(yPer) mySum += abs(yPer-yArr[index])*100sse = (yPer-yArr[index])**2 error = abs(yPer-yArr[index])/yArr[index]*100 #相對誤差 # print(yArr[index],round(yPer,2),str(round(error,2))+"%")plt.plot(index,error,"o")plt.title(" ",fontsize=13) #圖片上方留白 plt.rc('font',family='Arial') #設置字體 plt.rcParams['xtick.direction'] = 'in' #刻度線朝內 plt.rcParams['ytick.direction'] = 'in' plt.tick_params(labelsize=18) #刻度大小 plt.xlabel("Test No.",fontsize=18) plt.ylabel("Relative Error/%",fontsize=18) plt.savefig("線性回歸模型各擬合值相對誤差",dpi=500,bbox_inches = 'tight') #dpi-清晰度 plt.show()print("SSE=",sse,"平均相對誤差=",round(mySum/sum(yArr),2)) print(corrcoef(yPerList,yArr)[0][1])

?效果：

3.3.3 殘差圖

相關代碼：

#殘差圖 mySum = 0 sse = 0 yPerList = []#采用全部數(shù)據(jù)進行訓練 ws = standRegres(xArr,yArr) # print(ws) for index,x in enumerate(xArr):yPer = float(x*ws)residua = yArr[index] - yPerplt.plot(index,residua,"bo")x = np.linspace(0,25,100)plt.plot(x,np.zeros(len(x)),"r")plt.title(" ",fontsize=13) plt.rcParams['xtick.direction'] = 'in' plt.rcParams['ytick.direction'] = 'in' plt.rc('font',family='Arial') plt.tick_params(labelsize=18) plt.xlabel("Test No.",fontsize=18) plt.ylabel("Residua",fontsize=18) plt.ylim((-0.5, 0.5)) plt.savefig("殘差圖.jpg",dpi=500,bbox_inches = 'tight') plt.show()

效果：

3.3.4 預測值與真實值分布圖

相關代碼：

#分布圖 mySum = 0 sse = 0 yPerList = []#采用全部數(shù)據(jù)進行訓練 ws = standRegres(xArr,yArr) # print(ws) for index,x in enumerate(xArr):yPer = float(x*ws)yPerList.append(yPer) plt.plot(list(range(len(xArr))),yPerList,"bo",label="Fitted value") plt.plot(list(range(len(xArr))),yArr,"r*",label="Measurement value")plt.title(" ",fontsize=13) plt.rc('font',family='Arial') plt.rcParams['xtick.direction'] = 'in' plt.rcParams['ytick.direction'] = 'in' plt.tick_params(labelsize=18) plt.legend(framealpha=0,loc=(0.45, 0.55),fontsize=15) plt.xlabel("Test No.",fontsize=18) plt.ylabel("Oxygen solubility/mg·${{L^-}^1}$",fontsize=18) plt.savefig("分布圖.jpg",dpi=500,bbox_inches = 'tight') plt.show()

效果：

四、代碼-響應曲面

4.1 高斯消元

相關代碼：

#最小二乘法曲面擬合 def fun(x): round(x, 2)if x >= 0:return '+'+str(x)else:return str(x)def get_res(X, Y, Z, n):# 求方程系數(shù)sigma_x = 0for i in X: sigma_x += isigma_y = 0for i in Y: sigma_y += isigma_z = 0for i in Z: sigma_z += isigma_x2 = 0for i in X: sigma_x2 += i * isigma_y2 = 0for i in Y: sigma_y2 += i * isigma_x3 = 0for i in X: sigma_x3 += i * i * isigma_y3 = 0for i in Y: sigma_y3 += i * i * isigma_x4 = 0for i in X: sigma_x4 += i * i * i * isigma_y4 = 0for i in Y: sigma_y4 += i * i * i * isigma_x_y = 0for i in range(n):sigma_x_y += X[i] * Y[i]# print(sigma_xy)sigma_x_y2 = 0for i in range(n): sigma_x_y2 += X[i] * Y[i] * Y[i]sigma_x_y3 = 0for i in range(n): sigma_x_y3 += X[i] * Y[i] * Y[i] * Y[i]sigma_x2_y = 0for i in range(n): sigma_x2_y += X[i] * X[i] * Y[i]sigma_x2_y2 = 0for i in range(n): sigma_x2_y2 += X[i] * X[i] * Y[i] * Y[i]sigma_x3_y = 0for i in range(n): sigma_x3_y += X[i] * X[i] * X[i] * Y[i]sigma_z_x2 = 0for i in range(n): sigma_z_x2 += Z[i] * X[i] * X[i]sigma_z_y2 = 0for i in range(n): sigma_z_y2 += Z[i] * Y[i] * Y[i]sigma_z_x_y = 0for i in range(n): sigma_z_x_y += Z[i] * X[i] * Y[i]sigma_z_x = 0for i in range(n): sigma_z_x += Z[i] * X[i]sigma_z_y = 0for i in range(n): sigma_z_y += Z[i] * Y[i]# print("-----------------------")# 給出對應方程的矩陣形式a = np.array([[sigma_x4, sigma_x3_y, sigma_x2_y2, sigma_x3, sigma_x2_y, sigma_x2],[sigma_x3_y, sigma_x2_y2, sigma_x_y3, sigma_x2_y, sigma_x_y2, sigma_x_y],[sigma_x2_y2, sigma_x_y3, sigma_y4, sigma_x_y2, sigma_y3, sigma_y2],[sigma_x3, sigma_x2_y, sigma_x_y2, sigma_x2, sigma_x_y, sigma_x],[sigma_x2_y, sigma_x_y2, sigma_y3, sigma_x_y, sigma_y2, sigma_y],[sigma_x2, sigma_x_y, sigma_y2, sigma_x, sigma_y, n]])b = np.array([sigma_z_x2, sigma_z_x_y, sigma_z_y2, sigma_z_x, sigma_z_y, sigma_z])# 高斯消元解線性方程res = np.linalg.solve(a, b)return reslabelName = ["Oxygen solubility/mg·${{L^-}^1}$","T/$^\circ$C","pH","c(${{Fe^2}^+}$)/g·${{L^-}^1}$","c(${{Cu^2}^+}$)/g·${{L^-}^1}$","c(${{Fe^3}^+}$)/g·${{L^-}^1}$",] print(labelName)

4.2 響應曲面繪制代碼

繪圖的核心代碼在這里，感興趣的自己研究研究吧

變量res的系數(shù)就是用4.1節(jié)的高斯消元算法生成的系數(shù)

def matching_3D(X, Y, Z,xLabelIndex,yLabelIndex,name,arg1=37,arg2=-72):n = len(X)res = get_res(X, Y, Z, n)# 輸出方程形式print("z=%.6s*x^2%.6s*xy%.6s*y^2%.6s*x%.6s*y%.6s" % (fun(res[0]), fun(res[1]), fun(res[2]), fun(res[3]), fun(res[4]), fun(res[5])))# 畫曲面圖和離散點fig = plt.figure() # 建立一個空間ax = fig.add_subplot(111, projection='3d') # 3D坐標xgrid = np.linspace(min(X),max(X),100)ygrid = np.linspace(min(Y),max(Y),100)x,y = np.meshgrid(xgrid,ygrid)# 給出方程z = res[0] * x * x + res[1] * x * y + res[2] * y * y + res[3] * x + res[4] * y + res[5]# 畫出曲面sp = ax.plot_surface(x, y, z, rstride=3, cstride=3, cmap=cm.jet)ax.contourf(x,y,z,zdir='z',offset=5,cmap = plt.get_cmap('rainbow'))# 畫出點ax.scatter(X, Y, Z, c='r',label="實測點",alpha=0)plt.rc('font',family='Arial')plt.xlabel(labelName[xLabelIndex])plt.xticks(rotation=30,fontsize=9)plt.ylabel(labelName[yLabelIndex])ax.set_zlabel(labelName[0])# show_text(ax) # ax.legend()fig.colorbar(sp)ax.view_init(elev=arg1, azim=arg2)plt.rcParams['xtick.direction'] = 'in'plt.rcParams['ytick.direction'] = 'in'plt.savefig(name,dpi=500,bbox_inches = 'tight')fig.show()

4.3 調用matching_3D繪制響應曲面

代碼：

和4.2節(jié)matching_3D的參數(shù)對照看：

X與Y為自變量列表-分別為X軸與Y軸

yArr為因變量-Z軸坐標

1、4即為對應4.1節(jié)中l(wèi)abelName列表對應的名字 ，用于自動生成響應的坐標名稱

T-cu2+即為保存的圖片文件名

37，-72用于調整圖片視角

%matplotlib notebook X = [] Y = [] for x,y in zip(xArrMat[:,1],xArrMat[:,4]):X.append(float(x))Y.append(float(y)) matching_3D(X,Y,yArr,1,4,"T-cu2+曲面圖",37,-72)

效果：

?

總結

以上是生活随笔為你收集整理的论文绘图-教你如何绘制响应面的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。