多项式回归梯度下降法

前言：先看博主的一元线性回归有助于理解本篇思路

点击此处→一元线性回归梯度下降法(通俗易懂，初学专属)

一、数据准备

生成N个训练样本（X，Y），其中自变量X取值服从均值为0，方差为1 的正太分布.

多项式方程的式子为：

 以下我们以Y =5+2X+3X^2+4*X^3， N=20为例.

Y =5+2X+3X^2+4*X^3 + e_r, 其中 e_r为误差项，其取值服从均值为0，方差为1 的正太分布

创建数据后查看数据散点图，代码及数据散点图如下：

#创建数据集
a = 0
b=1
num=20
#创建样本，均值a=0，方差b=1，数量num为100
X = np.random.normal(a,b,num)
e_r = np.random.normal(a,b,num)
Y = 5 + 2*X + 3*X*X + 4*X*X*X + e_r
plt.scatter(X,Y)
plt.show()

 二、定义一个MultinomialModal 类

这次我们用矩阵运算，先看下图：

         一眼便能看出我们的计算思路。w和x用矩阵表示，w_matrix是(1，M)的矩阵，x_matrix是(M，n)的矩阵，这两矩阵相乘的结果正好就是(1, n)的矩阵，正好代表的就是y的值。

所以，我们接下来的步骤就是：

1.构建x_matrix矩阵

        构建该矩阵，我们使用拼凑的方法，numpy中有一个函数vstack()可以将两个数组拼凑一起，如下展示：

 

 我们一开始创建的X数据集是一个[X1, X2, X3,......,Xn]的一维数组，分别计算这个数组的0次方，1次方，2次方，然后拼凑成一个二维数组。写成一个函数，代码如下：

#数据准备，创建x的矩阵
def prepare(x, M):
    M = M +1        # M代表x的最高次方，但创建出来的矩阵应当是M+1行
    x_matrix = x**0
    if M ==1:
        x_matrix = np.atleast_2d(x_matrix)   # 升到二维
    for i in range(1, M):
        _ = x**i
        x_matrix = np.vstack((x_matrix,_))
    return x_matrix

M代表x的最高次方，但创建出来的矩阵应当是M+1行，比如x的最高次方为二，那么第一行应该是x的零次方，第二行是x的一次方，第三行是x的二次方。所以此处M要先加1。（此处还有一个问题，如果x的最高次数是0，即传进来的M是0，加一后M为1，会跳过循环直接以一维数组的形式返回，一维数据后面不能和二维数据进行乘法运算，所以要使用升维函数将其升到二维。）

2.根据梯度下降法的步骤，先求误差值loss

y = np.dot(w_matrix, x_matrix)

 3.求loss对每个w的偏导

        对w2求偏导的推导过程如图所示

 

        用同样的方式求w3，w4的偏导会发现规律 

        所以对w求偏导的公式为

 

        我们需要把x_matrix转置一下，让x的相同次方位于同一列。

        如下图仅列出第一列的求和结果，结果值就是w0的偏导

 

4.更新w_matrix

        公式与一元回归一样，对于每个w

        wj = wj - LearningRate*kwj

        kwj代表wj的偏导

5.正式定义MultinomialModal类

传入参数：数据x和真实值y，x的最高次方M，学习率LearningRate，构建好的x矩阵x_matrix

class MultinomialModal():
    def __init__(self, x, y, M, LearningRate, x_matrix):
        self.x = x
        self.y = y
        self.M = M
        self.LearningRate = LearningRate
        self.n = len(x)
        self.w_matrix = np.ones((1, M))      #初始化w矩阵
        self.x_matrix = x_matrix
    
        
    def train(self):
        for i in range(100000):
            #求预测值
            y_predict = np.dot(self.w_matrix,self.x_matrix)
            
            # 求误差值loss
            loss = np.sum((y_predict-self.y)**2)/self.n
            
            # 当误差值处于可接受范围则退出训练
            if loss <= 0.000001:
                break
                
            # 将x_matrix矩阵转置
            x_matrix_T = self.x_matrix.T
            # 求偏导
            dw_matrix = np.dot((y_predict-self.y), x_matrix_T)
            # 更新w_matrix
            self.w_matrix = self.w_matrix - self.LearningRate*dw_matrix
            
            

三、预测

#预测函数
def predict(w_matrix, x_matrix):
    y_predict = np.dot(w_matrix, x_matrix)
    return y_predict.flatten()

M = 3
# 准备x的矩阵
x_matrix = prepare(X, M)
# 实例化类
MM = MultinomialModal(X,Y, M+1, 0.0001, x_matrix)
# 开始训练
MM.train()
# 获取训练结果的w_matrix
w_matrix = MM.w_matrix
print(w_matrix)

# 绘图
plt.scatter(X,Y)
# 将x排序
X_sort = np.sort(X)
# 创建排序后的x的矩阵
X_sort_matrix = prepare(X_sort, M)
# 预测
y_predict = predict(w_matrix, X_sort_matrix)

plt.plot(X_sort, y_predict)
plt.show()

 可见w分别趋近5，2，3，4.证明拟合效果极佳

自此多项式回归梯度下降法完成

---

如有错误还请指教