多項式迴歸

如果你有一張圖，其中的點代表著某些資料的散佈，例如，底下是個 50 x 50 的圖：

目前看來，這張圖有呈現某個線性的趨勢，可以找出該線性的方程式嗎？這可以透過迴歸（Regression）來求，簡單來說，你想要求得 f_Θ(x) = Θ₀ + Θ₁ * x 的 Θ₀ 與 Θ₁，然後畫出下圖：

這邊使用圖上的點來代表資料散佈，純綷就是我個人很懶得收集或建些假資料，用圖可以隨便在上面點一點，就拿來計算罷了。

對於這個需求，第一個要實現的就是從圖中取得點的座標並用 Matplotlib 畫為散佈圖，這並不難，唯一要注意的是，圖片座標與 Matplotlib 的座標，在 y 軸方向是相反的：

import numpy as np
import matplotlib.pyplot as plt
import cv2

img = cv2.imread('data.jpg', cv2.IMREAD_GRAYSCALE)
idx = np.where(img < 127)       # 黑點的索引

data_x = idx[1]
data_y = -idx[0] + img.shape[0] # 反轉 y 軸

plt.gca().set_aspect(1)
plt.scatter(data_x, data_y)

plt.show()

這可以畫出以下的圖形：

對於多項式方面的需求，NumPy 從 1.4 以後，建議使用 numpy.polynomial 模組，若要求多項式迴歸，可以透過子模組 numpy.polynomial.polynomial 的 polyfit 函式，名稱代表多項式擬合，polyfit 指定資料 x、y 以及多項式級數，函式會傳回多項式的係數，例如：

import numpy as np
import matplotlib.pyplot as plt
import cv2
import numpy.polynomial.polynomial as p

img = cv2.imread('data.jpg', cv2.IMREAD_GRAYSCALE)
idx = np.where(img < 127)               # 黑點的索引

data_x = idx[1]
data_y = -idx[0] + img.shape[0]         # 反轉 y 軸

t0, t1 = p.polyfit(data_x, data_y, 1)   # 求係數

x = np.array([0, 50])
y = t0  + t1 * x

plt.gca().set_aspect(1)
plt.scatter(data_x, data_y)
plt.plot(x, y)

plt.show()

這就可以畫出方才的線性迴歸圖，當然，資料不一定是線性趨勢，視不同的資料而定，需要指定不同的級數，例如這張圖：

看來有二次曲線的趨勢，這時可以使用：

import numpy as np
import matplotlib.pyplot as plt
import cv2
import numpy.polynomial.polynomial as p

img = cv2.imread('data.jpg', cv2.IMREAD_GRAYSCALE)
idx = np.where(img < 127)       # 黑點的索引

data_x = idx[1]
data_y = -idx[0] + img.shape[0] # 反轉 y 軸

plt.gca().set_aspect(1)
plt.scatter(data_x, data_y)

t0, t1, t2 = p.polyfit(data_x, data_y, 2) # 求二次多項式的係數

x = np.linspace(0, 50, 50)
y = t0  + t1 * x + t2 * (x ** 2)

plt.plot(x, y)

plt.show()

這可以畫出這張圖：

簡單來說，觀察資料分佈的趨勢，採用適當的級數來進行函式的擬合。

如果想控制更多選項，可以透過 numpy.polynomial.Polynomial 的類別方法 fit，方法名稱，表示它可用來進行多項式擬合，fit 可以指定資料 x、y、多項式級數、domain、window 等參數，傳回 Polynomial 實例，例如，若只是需要繪圖：

import numpy as np
import matplotlib.pyplot as plt
import cv2
from numpy.polynomial import Polynomial as P

img = cv2.imread('data.jpg', cv2.IMREAD_GRAYSCALE)
idx = np.where(img < 127)       # 黑點的索引

data_x = idx[1]
data_y = -idx[0] + img.shape[0] # 反轉 y 軸

plt.gca().set_aspect(1)
plt.scatter(data_x, data_y)

p = P.fit(data_x, data_y, 2)    # 傳回 `Polynomial` 實例
x, y = p.linspace()             # 傳回 domain 範圍內的 x 與對應的 y

plt.plot(x, y)
plt.show()

Polynomial 實例代表一個多項式，以上例來說，傳回的 Polynomial 實例，其 domain 會是 [0, 49]，也就是以給定的 data_x 為 domain 範圍，也就是資料的定義域，linspace 方法會傳回 domain 範圍內的 x 與對應的 y。

如果想要透過 Polynomial 實例取多項式係數，可以透過 coef 特性，不過要注意，係數是受到 domain、window 影響，domain 範圍會對應至 window，也就是便於你平移或縮放定義域，在不指定 window 的情況下，window 預設是 [-1, 1]。

就上例來說，就是將 domain 的 [0, 49] 預設對應至 window 的 [-1, 1]，這邊並不打算平移或縮放定義域，因此必須指定 window 與 domain 相同。例如：

import numpy as np
import matplotlib.pyplot as plt
import cv2
from numpy.polynomial import Polynomial as P

img = cv2.imread('data.jpg', cv2.IMREAD_GRAYSCALE)
idx = np.where(img < 127)       # 黑點的索引

data_x = idx[1]
data_y = -idx[0] + img.shape[0] # 反轉 y 軸

plt.gca().set_aspect(1)
plt.scatter(data_x, data_y)

p = P.fit(data_x, data_y, 2, window = [0, img.shape[1]]) # 指定 window
t0, t1, t2 = p.coef  # 取得係數

x = np.linspace(0, 50, 50)
y = t0  + t1 * x + t2 * (x ** 2)

plt.plot(x, y)
plt.show()

如此與來，取得的係數用來繪圖，就會與上圖相同。