如何根據已有的數據擬合出多項式函數（多項式插值）？

02-12

實際項目中遇到了一個需求，就是根據一堆數據反過來求出一個函數，這個函數要表達的意思總結起來就是下面這樣：

x為10到無窮大；求一個函數f(x)，滿足當x=10, f(x) 約等於 0.98x=20, f(x) 約等於 0.96x=50, f(x) 約等於 0.95x=100, f(x) 約等於 0.90x=200, f(x) 約等於 0.85x=300, f(x) 約等於 0.80x=400, f(x) 約等於 0.75x=500, f(x) 約等於 0.70x>=1000, f(x) 約等於 0.60

有了這個函數後，我們就可以得到當x從10到無窮大的時候對應的f(x)值。那現在的問題就是如何根據這麼有限的幾個數據來擬合出所需的函數呢？

使用一次函數

我們使用numpy的polyfit函數來進行一次函數的擬合。

import numpy as npimport matplotlib.pyplot as pltCALIBRATE_X = [10, 20, 50, 100, 200, 300, 400, 500, 1000, 5000]CALIBRATE_Y = [0.98, 0.96, 0.95, 0.9, 0.85, 0.8, 0.75, 0.7, 0.6, 0.55]z1 = np.polyfit(CALIBRATE_X, CALIBRATE_Y, 1)p1 = np.poly1d(z1)#begin testx=range(6000)y=[]for i in x: y.append(p1(i))plt.plot(x,y)plt.show()

得出的圖像如下：

當然，一次函數是無法有效表達文章開頭提到的這些數據的特徵的，很明顯那些數據並不是線性的。

使用高次函數來擬合

我們用一個三次函數來進行擬合，那麼就有f(x) = $ax^{3} + bx^{2} + cx + d$ 。

於是我們得到下面的方程組：

$a10^{3} + b10^{2} + c10 + d = 0.98\ a20^{3} + b20^{2} + c20 + d = 0.96\ a50^{3} + b50^{2} + c50 + d = 0.95\ a100^{3} + b100^{2} + c100 + d = 0.90\ a200^{3} + b200^{2} + c200 + d = 0.85\ a300^{3} + b300^{2} + c300 + d = 0.80\ a400^{3} + b400^{2} + c400 + d = 0.75\ a500^{3} + b500^{2} + c500 + d = 0.70\ a1000^{3} + b1000^{2} + c1000 + d = 0.60$

然後我們就可以自然而然的使用矩陣運算來計算這個3次函數的係數：a、b、 c、 d。

import numpy as npimport matplotlib.pyplot as pltCALIBRATE_X = [10, 20, 50, 100, 200, 300, 400, 500, 1000]CALIBRATE_Y = [0.98, 0.96, 0.95, 0.9, 0.85, 0.8, 0.75, 0.7, 0.6]ORDER = 4x = np.array(CALIBRATE_X)y = np.array(CALIBRATE_Y)input_2d_x = np.ndarray((len(CALIBRATE_X), ORDER))#x^3 + x^2 + x^1 + x^0 if ORDER == 4for c in range(0, ORDER): input_2d_x[:, c] = x**cx_mat = np.mat(input_2d_x)y_mat = np.mat(y)#calculate the weigthsW = np.matmul(x_mat.I, y_mat.T)#begin the testdef fit(x): if x < 2: return 1 if x > 1000: return 0.6 input = [] for i in range(ORDER): input.append(x**i) x_mat = np.mat( input ) return float(np.matmul(x_mat, W))x = range(2000)y = []for i in x: y.append(fit(i))#show the image plt.plot(x, y)plt.show()