手写Python实现

笔者将利用numpy封装一个线性回归的算法类，其中会实现公式求解、矩阵求解和梯度下降求解三种方法，并分别对一元线性回归和多元线性回归的数据进行测试。

先看下算法类的主体，在train函数中可以选择"formula", "matrix", "gradient"三种方法进行模型的训练，分别代表"公式求解法", "矩阵求解法", "梯度下降法"，其中公式求解法仅限用于一元线性回归的求解。

import numpy as np

import matplotlib.pyplot as plt

class LinearRegression(object):

    """simple linear regression & multivariate linear regression"""

    def __init__(self):

        self.w = 0                  #斜率

        self.b = 0                  #截距

        self.sqrLoss = 0            #最小均方误差

        self.trainSet = 0           #训练集特征

        self.label = 0              #训练集标签

        self.learning_rate = None   #学习率

        self.n_iters = None         #实际迭代次数

        self.lossList = []          #梯度下降每轮迭代的误差列表

    def train(self, X, y, method, learning_rate=0.1, n_iters=1000):

        if X.ndim < 2:

            raise ValueError("X must be 2D array-like!")

        self.trainSet = X

        self.label = y

        if method.lower() == "formula":

            self.__train_formula()

        elif method.lower() == "matrix":

            self.__train_matrix()

        elif method.lower() == "gradient":

            self.__train_gradient(learning_rate, n_iters)

        else:

            raise ValueError("method value not found!")

        return

公式求解法的函数如下：

    #公式求解法(仅适用于一元线性回归)

    def __train_formula(self):

        n_samples, n_features = self.trainSet.shape

        X = self.trainSet.flatten()

        y = self.label

        Xmean = np.mean(X)

        ymean = np.mean(y)

        #求w

        self.w = (np.dot(X, y) - n_samples*Xmean*ymean)/(np.power(X,2).sum() - n_samples*Xmean**2)

        #求b

        self.b = ymean - self.w*Xmean

        #求误差

        self.sqrLoss = np.power((y-np.dot(X,self.w) - self.b), 2).sum()

        return

大林上位机机器视觉，_苏州电工培训_苏州PLC培训_苏州机器视觉培训_苏州上位机培训_苏州工业机器人培训，最适合电工及plc编程人员学习的上位机机器视觉课程 大林老师：15861139266（微信同号）

矩阵求解法的函数如下：

    #矩阵求解法

    def __train_matrix(self):

        n_samples, n_features = self.trainSet.shape

        X = self.trainSet

        y = self.label

        #合并w和b，在X尾部添加一列全是1的特征

        X2 = np.hstack((X, np.ones((n_samples, 1))))

        #求w和b

        EX = np.linalg.inv(np.dot(X2.T,X2))

        what = np.dot(np.dot(EX,X2.T),y)

        self.w = what[:-1]

        self.b = what[-1]

        self.sqrLoss = np.power((y-np.dot(X2,what).flatten()), 2).sum()

        return

梯度下降法相对前两者，参数上多了学习率，迭代次数，最小误差参数，迭代停止条件设置为“迭代达到指定次数”，或“迭代后的误差小于最小误差”，“参数更新前后的误差之差小于最小误差”。函数如下：

    #梯度下降法

    def __train_gradient(self, learning_rate, n_iters, minloss=1.0e-6):

        n_samples, n_features = self.trainSet.shape

        X = self.trainSet

        y = self.label

        #初始化迭代次数为0，初始化w0，b0为1，初始化误差平方和以及迭代误差之差

        n = 0

        w = np.ones(n_features)

        b = 1

        sqrLoss0 = np.power((y-np.dot(X,w).flatten()-b), 2).sum()

        self.lossList.append(sqrLoss0)

        deltaLoss = np.inf

        while (nminloss) and (abs(deltaLoss)>minloss):

            #求w和b的梯度

            ypredict = np.dot(X, w) + b

            gradient_w = -1.*np.dot((y - ypredict), X)/n_samples

            gradient_b = -1.*sum(y - ypredict)/n_samples

            #更新w和b

            w = w - learning_rate * gradient_w

            b = b - learning_rate * gradient_b

            #求更新后的误差和更新前后的误差之差

            sqrLoss1 = np.power((y-np.dot(X,w).flatten()-b), 2).sum()

            deltaLoss = sqrLoss0 - sqrLoss1

            sqrLoss0 = sqrLoss1

            self.lossList.append(sqrLoss0)

            n += 1

        print("第{}次迭代，损失平方和为{}，损失前后差为{}".format(n, sqrLoss0, deltaLoss))

        self.w = w

        self.b = b

        self.sqrLoss = sqrLoss0

        self.learning_rate = learning_rate

        self.n_iters = n+1

        return

对编写好的代码分别进行一元线性回归和多元线性回归的测试。

if __name__ == "__main__":

    #1、先用公式法和矩阵法测试下一元线性回归

    simpleLR(1.34, 2.08)

    #2、再用矩阵法和梯度下降法测试下多元线性回归

    multivariateLR()

一元线性回归测试：自己指定w和b，造了一组数据，并加入了误差，然后分别公式法和矩阵法进行求解。

def simpleLR(w, b, size=100):

    X = np.expand_dims(np.linspace(-10, 10, size), axis=1)

    y = X.flatten()*w + b + (np.random.random(size)-1)*3

    #公式法求解

    lr1 = LinearRegression()

    lr1.train(X, y, method='formula')

    print("【formula方法】\nw:{}, b:{}, square loss:{}".format(lr1.w, lr1.b, lr1.sqrLoss))

    #矩阵法求解

    lr2 = LinearRegression()

    lr2.train(X, y, method='Matrix')

    print("【matrix方法】\nw:{}, b:{}, square loss:{}".format(lr2.w, lr2.b, lr2.sqrLoss))

    #画图

    fig = plt.figure()

    ax = fig.add_subplot(1,1,1)

    ax.scatter(X, y)

    ax.plot(X, X*lr2.w+lr2.b, color='r', linewidth=3)

    plt.show()

    return