读源码理解线性回归

  提到线性回归，似乎感觉很简单，不就是用一条线作为一个模型来预测，实际上根本没有我们想的那么简单。

 为何这么说呢?

主要就是我们怎么选择一条合适的线做为模型呢?什么时候迭代结束呢？

loss函数会给我们答案，但是loss又不是仅仅是最小二乘法就足够的，因为我们要考虑过拟合，这时就会用到正则项。

也就是我们真正在训练的过程中的loss=loss（损失）+L（正则项）

正则项分为两种一个是L1，一个是L2

应用L1称为LASSO回归

应用L2称为岭回归

和逻辑回归很像，我们需要用权重将很多个x进行加权后求和，然后再和真实值进行对比，整个训练的过程就是调整权重的过程，因此调整权重也是使用的梯度方法，唯一不同的就是有了正则项，正则项看成损失就可以很好的类比了

L1的的损失函数就是将权重的绝对值进行求和，然后乘alpha

L2的的损失函数就是将权重进行平方后求和，然后乘alpha再除以2

梯度也很好理解，x的绝对值求导就是sign(x)，x的平方求导2x

 

不仅仅是loss

  我们在进行线性回归的时候，也是可以选择不进行迭代的，此时不需要loss，我们仅仅通过矩阵的变换计算得到，这种方法叫做正规方程。我们最终的目的是求权重w，我们可以推到一下，利用矩阵的方法：

 我们可以从公式中看到我们需要的权重其实是可以通过矩阵变换得到的，当然这仅仅限于可逆对于这个矩阵

  如果矩阵不是特别大，其实用这个效果不仅仅好，而且计算也不会很大，如果矩阵很大这个方法就不适用了，毕竟需要计算很庞大的矩阵运算。

附上人家的代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import make_regression

from utils import train_test_split
from utils import mean_squared_error, Plot

# L1正则化
class l1_regularization():
    def __init__(self, alpha):
        self.alpha = alpha

    # L1正则化的方差
    def __call__(self, w):
        loss = np.sum(np.fabs(w))
        return self.alpha * loss

    # L1正则化的梯度
    def grad(self, w):
        return self.alpha * np.sign(w)


# L2正则化
class l2_regularization():
    def __init__(self, alpha):
        self.alpha = alpha

    # L2正则化的方差
    def __call__(self, w):
        loss = w.T.dot(w)
        return self.alpha * 0.5 * float(loss)

    # L2正则化的梯度
    def grad(self, w):
        return self.alpha * w


class LinearRegression():
    """
    Parameters:
    -----------
    n_iterations: int
        梯度下降的轮数
    learning_rate: float
        梯度下降学习率
    regularization: l1_regularization or l2_regularization or None
        正则化
    gradient: Bool
        是否采用梯度下降法或正规方程法。
        若使用了正则化，暂只支持梯度下降
    """

    def __init__(self, n_iterations=3000, learning_rate=0.00005, regularization=None, gradient=True):
        self.n_iterations = n_iterations
        self.learning_rate = learning_rate
        self.gradient = gradient
        if regularization == None:
            self.regularization = lambda x: 0
            self.regularization.grad = lambda x: 0
        else:
            self.regularization = regularization

    def initialize_weights(self, n_features):
        # 初始化参数
        limit = np.sqrt(1 / n_features)
        w = np.random.uniform(-limit, limit, (n_features, 1))
        b = 0
        self.w = np.insert(w, 0, b, axis=0)

    def fit(self, X, y):
        m_samples, n_features = X.shape
        self.initialize_weights(n_features)
        X = np.insert(X, 0, 1, axis=1)
        y = np.reshape(y, (m_samples, 1))
        self.training_errors = []
        if self.gradient == True:
            # 梯度下降
            for i in range(self.n_iterations):
                y_pred = X.dot(self.w)
                loss = np.mean(0.5 * (y_pred - y) ** 2) + self.regularization(self.w) #计算loss
                # print(loss)
                self.training_errors.append(loss)
                w_grad = X.T.dot(y_pred - y) + self.regularization.grad(self.w)  # (y_pred - y).T.dot(X)，计算梯度
                self.w = self.w - self.learning_rate * w_grad #更新权值w
        else:
            # 正规方程
            X = np.matrix(X)
            y = np.matrix(y)
            X_T_X = X.T.dot(X)
            X_T_X_I_X_T = X_T_X.I.dot(X.T)
            X_T_X_I_X_T_X_T_y = X_T_X_I_X_T.dot(y)
            self.w = X_T_X_I_X_T_X_T_y

    def predict(self, X):
        X = np.insert(X, 0, 1, axis=1)
        y_pred = X.dot(self.w)
        return y_pred





def main():
    X, y = make_regression(n_samples=100, n_features=1, noise=20)

    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
    n_samples, n_features = np.shape(X)

    # 可自行设置模型参数，如正则化，梯度下降轮数学习率等
    model = LinearRegression(n_iterations=3000, regularization=l2_regularization(alpha=0.5))

    model.fit(X_train, y_train)

    # Training error plot 画loss的图
    n = len(model.training_errors)
    training, = plt.plot(range(n), model.training_errors, label="Training Error")
    plt.legend(handles=[training])
    plt.title("Error Plot")
    plt.ylabel('Mean Squared Error')
    plt.xlabel('Iterations')
    plt.show()

    y_pred = model.predict(X_test)
    y_pred = np.reshape(y_pred, y_test.shape)

    mse = mean_squared_error(y_test, y_pred)
    print("Mean squared error: %s" % (mse))

    y_pred_line = model.predict(X)

    # Color map
    cmap = plt.get_cmap('viridis')

    # Plot the results，画拟合情况的图
    m1 = plt.scatter(366 * X_train, y_train, color=cmap(0.9), s=10)
    m2 = plt.scatter(366 * X_test, y_test, color=cmap(0.5), s=10)
    plt.plot(366 * X, y_pred_line, color='black', linewidth=2, label="Prediction")
    plt.suptitle("Linear Regression")
    plt.title("MSE: %.2f" % mse, fontsize=10)
    plt.xlabel('Day')
    plt.ylabel('Temperature in Celcius')
    plt.legend((m1, m2), ("Training data", "Test data"), loc='lower right')
    plt.show()


if __name__ == "__main__":
    main()