Python中多元5度多項式迴歸的曲面圖

Question

我正在用Python實現一篇文章，最初是在MATLAB中實現的。該論文指出，使用來自一組採樣數據點的曲線擬合找到五次多項式。我不想使用它們的多項式，所以我開始使用樣本數據點（在紙中給出），並嘗試使用sklearn多項式特徵和linear_model來查找5度多項式。因爲它是一個多元方程f（x，y），其中x和y是某個池塘的長度和寬度，f是污染物的初始濃度。

所以我的問題是，sklearn多項式特徵將測試和訓練數據點轉換爲n多項式點（我認爲據我瞭解）。但是當我需要clf.predict函數（其中clf是訓練模型）只取x和y值，因爲當我從Matplotlib繪製表面圖時，它需要meshgrid，所以當我meshgrid my sklean轉換測試點，它的形狀變得像NxN，而預測函數需要Nxn（其中n是它轉換數據的多項式的次數），N是行數。

是否有任何可能的方法來繪製這個多項式的網格點？

紙業鏈接：http://www.ajer.org/papers/v5(11)/A05110105.pdf 文章：基於二維 對流擴散模型

如果可能的話在兼廢水穩定塘生物需氧量的數學模型，請看圖5和6的紙（上面的鏈接）。這是我想要實現的。

謝謝 enter code here

from math import exp 
import numpy as np 
from operator import itemgetter 
from sklearn.preprocessing import PolynomialFeatures 
from sklearn import linear_model 
from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt 
from matplotlib import cm 
from matplotlib.ticker import LinearLocator, FormatStrFormatter 

fig = plt.figure() 
ax = fig.gca(projection='3d') 


def model_BOD (cn): 
    cnp1 = [] 
    n = len(cn) 
    # variables: 
    dmx = 1e-5 
    dmy = 1e-5 
    u = 2.10e-4 
    v = 2.10e-4 
    obs_time = 100 
    dt = 0.1 

    for t in np.arange(0.1,obs_time,dt): 
     for i in range(N): 
      for j in range(N): 
       d = j + (i-1)*N 
       dxp1 = d + N 
       dyp1 = d + 1 
       dxm1 = d - N 
       dym1 = d - 1 

       cnp1.append(t*(((-2*dmx/dx**2)+(-2*dmy/dy**2)+(1/t))*cn[dxp1] + (dmx/dx**2)*cn[dyp1] \ 
           + (dmy/dy**2)*cn[dym1] - (u/(2*dx))*cn[dxp1] + (u/(2*dx))*cn[dxm1] \ 
           - (v/(2*dy))*cn[dyp1] + (v/(2*dy))*cn[dym1])) 
     cn = cnp1 
     cnp1 = [] 
    return cn 

N = 20 
Length = 70 
Width = 77 
dx = Length/N 
dy = Width/N 

deg_of_poly = 5 

datapoints = np.array([ 
    [12.5,70,81.32],[25,70,88.54],[37.5,70,67.58],[50,70,55.32],[62.5,70,56.84],[77,70,49.52], 
    [0,11.5,71.32],[77,57.5,67.20], 
    [0,23,58.54],[25,46,51.32],[37.5,46,49.52], 
    [0,34.5,63.22],[25,34.5,48.32],[37.5,34.5,82.30],[50,34.5,56.42],[77,34.5,48.32],[37.5,23,67.32], 
    [0,46,64.20],[77,11.5,41.89],[77,46,55.54],[77,23,52.22], 
    [0,57.5,93.72], 
    [0,70,98.20],[77,0,42.32] 
    ]) 

X = datapoints[:,0:2] 
Y = datapoints[:,-1] 

predict_x = [] 
predict_y = [] 
for i in np.linspace(0,Width,N): 
    for j in np.linspace(0,Length,N): 
     predict_x.append([i,j]) 

predict = np.array(predict_x) 

poly = PolynomialFeatures(degree=deg_of_poly) 
X_ = poly.fit_transform(X) 

predict_ = poly.fit_transform(predict) 
clf = linear_model.LinearRegression() 
clf.fit(X_, Y) 
prediction = [] 

for k,i in enumerate(predict_): 
    prediction.append(clf.predict(np.array([i]))[0]) 

prediction_ = model_BOD(prediction) 
print prediction_ 
XX = [] 
XX = predict[:,0] 
YY = [] 
YY = predict[:,-1] 
XX,YY = np.meshgrid(X,Y) 
Z = prediction 
##R = np.sqrt(XX**2+YY**2) 
##Z = np.tan(R) 

surf = ax.plot_surface(XX,YY,Z) 
plt.show() 

Answer 1

如果我理解正確的話，這裏的關鍵邏輯是從你的meshgrid生成多項式的特點，使預測與原meshgrid繪製預測。希望下面給你想要的東西：

import matplotlib.pyplot as plt 
from matplotlib import cm 
from mpl_toolkits.mplot3d import Axes3D 
import numpy as np 
from sklearn.preprocessing import PolynomialFeatures 
from sklearn import linear_model 

# The training set 
datapoints = np.array([ 
    [12.5,70,81.32], [25,70,88.54], [37.5,70,67.58], [50,70,55.32], 
    [62.5,70,56.84], [77,70,49.52], [0,11.5,71.32], [77,57.5,67.20], 
    [0,23,58.54], [25,46,51.32], [37.5,46,49.52], [0,34.5,63.22], 
    [25,34.5,48.32], [37.5,34.5,82.30], [50,34.5,56.42], [77,34.5,48.32], 
    [37.5,23,67.32], [0,46,64.20], [77,11.5,41.89], [77,46,55.54], 
    [77,23,52.22], [0,57.5,93.72], [0,70,98.20], [77,0,42.32] 
    ]) 
X = datapoints[:,0:2] 
Y = datapoints[:,-1] 
# 5 degree polynomial features 
deg_of_poly = 5 
poly = PolynomialFeatures(degree=deg_of_poly) 
X_ = poly.fit_transform(X) 
# Fit linear model 
clf = linear_model.LinearRegression() 
clf.fit(X_, Y) 

# The test set, or plotting set 
N = 20 
Length = 70 
predict_x0, predict_x1 = np.meshgrid(np.linspace(0, Length, N), 
            np.linspace(0, Length, N)) 
predict_x = np.concatenate((predict_x0.reshape(-1, 1), 
          predict_x1.reshape(-1, 1)), 
          axis=1) 
predict_x_ = poly.fit_transform(predict_x) 
predict_y = clf.predict(predict_x_) 

# Plot 
fig = plt.figure(figsize=(16, 6)) 
ax1 = fig.add_subplot(121, projection='3d') 
surf = ax1.plot_surface(predict_x0, predict_x1, predict_y.reshape(predict_x0.shape), 
         rstride=1, cstride=1, cmap=cm.jet, alpha=0.5) 
ax1.scatter(datapoints[:, 0], datapoints[:, 1], datapoints[:, 2], c='b', marker='o') 

ax1.set_xlim((70, 0)) 
ax1.set_ylim((0, 70)) 
fig.colorbar(surf, ax=ax1) 
ax2 = fig.add_subplot(122) 
cs = ax2.contourf(predict_x0, predict_x1, predict_y.reshape(predict_x0.shape)) 
ax2.contour(cs, colors='k') 
fig.colorbar(cs, ax=ax2) 
plt.show()