紫梦沁香
0%

多项式回归

发表于 更新于 分类于 阅读次数： Valine：
本文字数： 3.7k 阅读时长 ≈ 7 分钟

算法思想

如果你的数据比直线更复杂怎么办？例如抛物线。显然 线性回归 无法胜任这个工作，但我们只需要稍加修改，将每个特征的幂次方作为一个新特征即可，大致与线性回归类似，不过是多元一次和多元多次的区别。

线性回归假设目标拟合曲线函数为$y = w_1x_1 + w_2x_2 + w_3x_3 + … + b$，其中 $w_n$ 为系数 $x_n$ 为特征。

而多项式回归假设目标拟合曲线为$y = w_1x_1^n + w_2x_1^{n-1}x_2 + w_3x_1^{n-2}x_2^2 + … + w_nx_2^n$，这里出于简便只写了两个未知数的情况。


实例

这里我用的是三次函数的图像外加了点噪声，用线性回归肯定是拟合不了的。

1
2
3
4
5
6
7
8
9
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.arange(-10, 10, 0.1)
y = 2.6 * x ** 3 - 3.14 * x ** 2 - 8.8 * x + np.random.normal(0, 700, (200, ))

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

1
2
3
4
5
X = x.reshape([-1, 1])
X = np.hstack([X ** 3, X ** 2, X])
X_b = np.hstack([np.ones([len(y), 1]), X])
# X [200, 3]
print(X.shape)

如算法思想所述，梯度下降的代码是没有改变的，只是添加了特征的幂指数作为额外的特征值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def J(X_b, y, theta):
    return np.mean((X_b.dot(theta) - y) ** 2);

def DJ(X_b, y, theta):
    return X_b.T.dot((X_b.dot(theta) - y)) * 2 / len(y);

def gradient_descent(X_b, y, theta, eta = 0.01, n_iters = 1e4, epsilon = 1e-6):
    i_iter = 0
    
    while i_iter < n_iters:
        gradient = DJ(X_b, y, theta)
        last_theta = theta
        theta = theta - eta * gradient
        if np.abs(J(X_b, y, theta) - J(X_b, y, last_theta)) < epsilon:
            break
        i_iter += 1
        
    return theta
1
2
initial_theta = np.zeros(X_b.shape[1])
theta = gradient_descent(X_b, y, initial_theta, eta = 1e-6, n_iters = 1e5)

得出 theta = array([-12.65120346, 2.89719563, -3.12234806, -15.72650818])

也基本对应了我预设的函数系数[2.6, -3.14, -8.8]

注：theta 每一项分别对应了截距，$x^3$ 的系数，$x^2$ 的系数，$x$ 的系数。

1
2
3
4
pred_y = theta[1] * x ** 3 + theta[2] * x ** 2 + theta[3] * x + theta[0]
plt.scatter(x, y)
plt.plot(x, pred_y, color = "r")
plt.show()


代码实现

需要注意的是多项式化是成阶乘式的增长，比指数爆炸更可怕！

$\frac{(n + d)!}{d!n!}$，其中 n 为原本特征数，d 为最高次数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
class PolynomialRegession:
    
    # 生成多项式特征
    def generatePolynomialFeatures(self, X, degree):
        self.degree = degree
        if degree == 0:
            return np.ones_like(X)
        
        n = X.shape[1]
        ans = np.array(X)
        if n > 1:
            for i in range(n):
                ans = np.hstack([ans, (X[:, i] ** degree).reshape(-1,1)])
                for j in range(i+1, n):
                    for k in range(1, degree):
                        ans = np.hstack([ans, (X[:, i] ** k * X[:, j] ** (degree - k)).reshape(-1,1)])
        else:
            for k in range(2, degree + 1):
                ans = np.hstack([ans, X ** k])

        return np.hstack([ans, np.ones(ans.shape[0]).reshape([-1, 1])])
    
    # 均值方差归一化 -> avg = 0, std = 1
    def standardization(self, X):
        tmp = X[:, :-1]
        return np.hstack([(tmp - np.mean(tmp, axis=0)) / np.std(tmp, axis = 0), np.ones([len(X), 1])])
    
    # 使用MSE作为损失函数
    def f(self, X_b, y, theta):
        return np.mean((X_b @ theta - y) ** 2)

    # 损失函数的导数
    def df(self, X_b, y, theta):
        return X_b.T @(X_b @ theta - y) * 2 / len(y)

    # 梯度下降
    def gd(self, X_b, y, theta, eta = 0.01, n_iters = 1e4, epsilon = 1e-6, plot_process = True):
        i_iter = 0

        while i_iter < n_iters:
            if ((plot_process and (i_iter + 1) % 1000 == 0) and (i_iter + 1) <= 10000) or (i_iter + 1) == n_iters:
                display.clear_output(wait=True)
                plt.scatter(x, y1)
                plt.plot(x, (X_b @ theta).reshape(x.shape), c="green", label = 'predict')
                plt.plot(x, y2, c="red", label = 'truth')
                plt.legend()
                plt.title("degree = " + str(self.degree))
                plt.show()
                plt.pause(0.5)

            gradient = self.df(X_b, y, theta)
            last_theta = theta
            theta = theta - eta * gradient
            if abs(self.f(X_b, y, last_theta) - self.f(X_b, y, theta)) < epsilon:
                break
            i_iter += 1

        return theta
    
    # 拟合
    def fit(self, X, y, degree, eta = 1e-4, n_iters = 1e4, epsilon = 1e-6, plot_process = True):
        self.X_b = self.generatePolynomialFeatures(X, degree)
        self.X_b = self.standardization(self.X_b)
        initial_theta = np.zeros_like(self.X_b[1])
        self.theta = self.gd(self.X_b, y, initial_theta, eta=eta, n_iters=n_iters, epsilon=epsilon, plot_process=plot_process)
        self.coefficient = self.theta[:-1]
        self.intercept = self.theta[-1]
    
    # 预测
    def predict(self, X):
        return self.standardization(self.generatePolynomialFeatures(X, self.degree) @ self.theta)
    
    # 预测分析
    def score(self, X, y):
        return np.mean((predict(X) - y) ** 2)

测试

1
2
3
4
5
6
7
8
poly_reg = PolynomialRegression(3)
poly_reg.fit(x.reshape([-1, 1]), y, eta = 1e-6, n_iters = 1e5)

theta = poly_reg.theta
pred_y = theta[3] * x ** 3 + theta[2] * x ** 2 + theta[1] * x + theta[0]
plt.scatter(x, y)
plt.plot(x, pred_y, color = "r")
plt.show()

结果如出一辙。