【python机器学习】普通最小二乘法多元线性回归

普通最小二乘法线性回归

若数据集$D$由$n$个属性描述，则线性回归的假设函数为：
$$h_{\mathbf{w},b}(\mathbf{x})=\sum _{i=1}^n{w}_i{x}_i+b={\mathbf{w}}^{\mathrm{T}}\mathbf{x}+b$$
其中，$\mathbf{w}\in {\mathbb{R}}^n$与$b\in \mathbb{R}$为模型参数。

为了方便，我们通常将$b$纳入权向量$\mathbf{w}$，作为$w_0$，同时为输入向量$\mathbf{x}$添加一个常数1，作为$x_0$.
$$\begin{array}{c}\mathbf{w}={\left(b,w_1,w_2,\dots w_n\right)}^{\mathrm{T}}\\ \mathbf{x}={\left(1,x_1,x_2,\dots x_n\right)}^{\mathrm{T}}\end{array}$$

此时，假设函数为：
$$h_{\mathbf{w}}(\mathbf{x})=\sum _{i=0}^n{w}_i{x}_i={\mathbf{w}}^{\mathrm{T}}\mathbf{x}$$

其中，$\mathbf{w}\in {\mathbb{R}}^{n+1}$，通过训练确定模型参数$\mathbf{w}$后，便可使用模型对新的输入实例进行预测。

使用均方误差（MSE）作为损失函数，假设训练集$D$有$m$个样本，均方误差损失函数定义为
$$\begin{array}{rl}J(\mathbf{w})& =\frac{1}{2m}\sum _{i=1}^m{\left(h_{\mathbf{w}}\left({\mathbf{x}}_i\right)-y_i\right)}^2\\ & =\frac{1}{2m}\sum _{i=1}^m{\left({\mathbf{w}}^{\mathrm{T}}\mathbf{x}-y_i\right)}^2\end{array}$$

损失函数$J(w)$最小值点是其极值点，可先求$J(w)$对$w$的梯度并令其为0，再通过解方程求得。

计算$J(\mathbf{w})$的梯度：
$$\begin{array}{rl}\nabla J(\mathbf{w})& =\frac{1}{2m}\sum _{i=1}^m\frac{\partial }{\partial \mathbf{w}}{\left({\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i-y_i\right)}^2\\ & =\frac{1}{2m}\sum _{i=1}^{m}2\left({\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i-y_i\right)\frac{\partial }{\partial \mathbf{w}}\left({\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i-y_i\right)\\ & =\frac{1}{m}\sum _{i=1}^m\left({\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i-y_i\right){\mathbf{x}}_i\end{array}$$

以上公式使用矩阵运算描述形式更为简洁，设：
$$\mathbf{X}=\left[\begin{array}{ccccc}1,& x_{11},& x_{12}& \dots & x_{1n}\\ 1,& x_{21}& x_{22}& \dots & x_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1,& x_{m1}& x_{m2}& \dots & x_{mn}\end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}_1^{\mathrm{T}}\\ {\mathbf{x}}_2^{\mathrm{T}}\\ ⋮\\ {\mathbf{x}}_m^{\mathrm{T}}\end{array}\right]$$
$$\mathbf{y}=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$$

$$\mathbf{w}=\left[\begin{array}{c}b\\ w_1\\ w_2\\ ⋮\\ w_n\end{array}\right]$$

那么，梯度计算公式可写为：
$$\nabla J(\mathbf{w})=\frac{1}{m}\sum _{i=1}^m\left({\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_i-y_i\right){\mathbf{x}}_i$$
$$=\left[\begin{array}{c}{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_m\end{array}\right]\left[\begin{array}{c}{\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_1-y_1\\ {\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_2-y_2\\ ⋮\\ {\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_m-y_m\end{array}\right]$$
$$=\left[\begin{array}{c}{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_m\end{array}\right]\left(\left[\begin{array}{c}{\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_1\\ {\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_2\\ ⋮\\ {\mathbf{w}}^{\mathrm{T}}{\mathbf{x}}_m\end{array}\right]-\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]\right)$$
$$=\frac{1}{m}{\mathbf{X}}^{\mathrm{T}}(\mathbf{X}\mathbf{w}-\mathbf{y})$$
令梯度为0，解得：
$$\hat{\mathbf{w}}={\left({\mathbf{X}}^{\mathrm{T}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathrm{T}}\mathbf{y}$$

$\hat{\mathbf{w}}$即为使得损失函数（均方误差）最小的$\mathbf{w}$。以上求解最优$\mathbf{w}$的方法被称为普通最小二乘法（Ordinary Least Squares，OLS）。

import numpy as np


class OLSLinearRession:
    def _ols(self, X, y):
        '''普通最小二乘法估算w'''
        tmp = np.linalg.inv(np.matmul(X.T, X))
        tmp = np.matmul(tmp, X.T)
        w = np.matmul(tmp, y)
        return w


    def _preprocess_data(self, X):
        '''数据预处理:添加x0=1'''
        m, n = X.shape
        X_ = np.ones((m, n + 1))
        X_[:, 1:] = X
        return X_

    def train(self, X, y):
        '''训练模型'''
        X = self._preprocess_data(X)
        self.w = self._ols(X, y)

    def predict(self, X):
        '''预测'''
        X = self._preprocess_data(X)
        y = np.matmul(X, self.w)
        return y