线性回归算法

秦枟

已于 2022-03-16 15:36:04 修改

阅读量348

点赞数

文章标签：机器学习

于 2019-04-22 16:32:00 首次发布

本文链接：https://blog.csdn.net/a766543984/article/details/89341350

版权

线性回归

理论
实现

理论

思想：用一条直线拟合样本点，让残差平方和最小。

模型： $h(x)=X\theta$

参数： $\theta=(\theta_0,...,\theta_n)^{'}$

目标函数： $J(\theta)=\frac{1}{2m}\sum^{m}_{i=1}(y^{(i)}-h(x^{(i)}))^2=\frac{1}{2m}||y-X\theta||_2^2$

梯度： $\bigtriangledown J(\theta)=\frac{1}{m}X^T(h(x)-y)$

最优参数： $\hat{\theta}=(X^TX)^{-1}X^Ty$

优点：

建模速度快，不需要很复杂的计算，在数据量大的情况下依然运行速度很快。
可以根据系数给出每个变量的解释。

缺点：

不能很好地拟合非线性数据。所以需要先判断变量之间是否是线性关系。

1.最小二乘法(代数视角)

一种思路是，对于样本，有观测值 $y_i$ 和理论值 $\hat{y_i}$ ，两者的差就是残差 $e_i$ 。我们让总残差最小，比如使均方误差(MSE)最小，这就是最小二乘法的主要思想。

小样本OLS基本假定

线性假定：每个解释变量对被解释变量的边际效应均为常数。
严格外生性： $E(\epsilon | X) = E(\epsilon_i | x_1, ..., x_n) = 0$ ，即扰动项与所有解释变量均不相关。
不存在严格多重共线性：即X满列秩，则 $X^{'}X$ 为正定矩阵。
球形扰动项：扰动项满足‘同方差’、‘无自相关’等性质。
$var(\epsilon | X) = E(\epsilon ^{'}\epsilon | X) = \sigma ^2I = \left\{ \begin{matrix} \sigma^2 & \cdots & 0\ \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \sigma^2 \end{matrix} \right\}$

大样本OLS的基本假定

线性假定（不变）
同期外生性：不再要求扰动项对所有解释变量正交，只是要求扰动项对同期的解释变量正交。
无严格多重共线性（不变）
渐近独立的平稳过程：保证样本均值是总体均值的一致估计。

公式推导

$J(\beta) = \frac{1}{m}\sum^{m}_{i=1}(y^{(i)}-X^{(i)}\hat{\beta})^2\\ =\frac{1}{m}(e^{'})^2\\ =\frac{1}{m}\boldsymbol{e^{'}e}\\ =\frac{1}{m}\boldsymbol{(Y - X\hat{\beta})^{'}(Y - X\hat{\beta})}\\ =\frac{1}{m}(\boldsymbol{Y^{'}Y-Y^{'}X\hat{\beta}-\hat{\beta}^{'}X^{'}Y+\hat{\beta}^{'}X^{'}X\hat{\beta}})$

$\frac{\partial J(\beta)}{\partial \hat{\beta}}=\frac{\partial (\boldsymbol{Y^{'}Y-Y^{'}X\hat{\beta}-\hat{\beta}^{'}X^{'}Y+\hat{\beta}^{'}X^{'}X\hat{\beta}}}{\partial \hat{\beta}}\\ =\frac{\partial tr(\boldsymbol{Y^{'}Y-Y^{'}X\hat{\beta}-\hat{\beta}^{'}X^{'}Y+\hat{\beta}^{'}X^{'}X\hat{\beta}})}{\partial \hat{\beta}}\\ =-(Y^{'}X)^{'}-X^{'}Y+X^{'}X\hat{\beta}+(\hat{\beta}^{'}X^{'}X)^{'}\\ =-2X^{'}Y+2X^{'}X\hat{\beta}$

令 $\frac{\partial J(\beta)}{\partial \hat{\beta}}=0$ ，则 $\hat{\beta}=(X^{'}X)^{-1}X^{'}Y$

矩阵的迹求导

2.梯度下降法

梯度下降法是一种基于搜索的最优化方法，作用是最小化一个损失函数(最大化效用函数用梯度上升法)。

我们首先来明晰两个概念
方向导数：函数 $z = f (x, y)$ 在点 $x_0,y_0)$ 沿方向 $\vec{l}$ 的方向导数

$\frac{\partial f}{\partial l}=\lim_{\rho \rightarrow 0^+}\frac{f(x_0+\rho \cos{\alpha}, y_0+\rho \cos{\beta})-f(x_0,y_0)}{\rho}$

$\frac{\partial f}{\partial l}$ 是函数 $z$ 对点 $x_0,y_0)$ 沿方向 $\vec{l}$ 对 $\rho$ 的变化率，也是曲面 $z$ 在点 $x_0,y_0)$ 沿方向 $\vec{l}$ 的倾斜程度。

梯度：向量 $(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})$ 是使 $f (x, y)$ 在一点增加最快的方向，称向量 $(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})$ 为可微函数 $z = f (x, y)$ 在点 $(x, y)$ 处的梯度向量，简称梯度。

记作：
$f=(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})=\frac{\partial f}{\partial x}i+\frac{\partial f}{\partial y}j=\bigtriangledown f$

梯度 $g r a d f$ 是一个向量，是可微函数 $z = f (x, y)$ 在点 $(x, y)$ 处取得最大方向导数的方向(即函数增加最快的方向)。

原理

对于损失函数 $J(\beta) = \frac{1}{2m}\sum^{m}_{i=1}(y^{(i)}-X^{(i)}\beta)^2=\frac{1}{2m}\sum^{m}_{i=1}(y^{(i)}-\sum^{n}_{j=0}\beta_jx^{(i)}_j)^2$ ，

梯度 $\bigtriangledown J(\beta)=(\frac{\partial J}{\partial \beta_0},...,\frac{\partial J}{\partial \beta_n})^{'}=-\frac{1}{m}(...,\sum^{m}_{i=1}(y^{(i)}-X^{(i)}\beta)x^{(i)}_j,...)^{'}j\in{(0,n)}$

用矩阵形式表示

$\bigtriangledown J(\beta)=\frac{1}{m}X^{'}(X\beta - Y)$

梯度的方向就是损失函数 $J$ 上升最快的方向，要使 $J$ 最小，则 $\beta := \beta - \eta \bigtriangledown J$ ，其中 $\eta$ 是步长，也称学习率。即 $\beta$ 沿损失函数 $J$ 下降最快的方向移动 $\eta \bigtriangledown J$ ，直至梯度为0。

3.最大似然法(概率视角)

最大似然估计的思想，在于已知模型，最有可能产生观测的样本的参数是多少。

$X\beta + \epsilon$ ，假定 $var(\epsilon | X) = \sigma ^2I$ ，则 $\sim N(X\beta, \sigma^2I)$

$f(y^{(i)}|X^{(i)})=\frac{1}{\sqrt{2\pi} \sigma}\exp\{-\frac{(y^{(i)}-X^{(i)'}\beta)^2}{2\sigma^2}\}$

效用函数：

$L(\beta)=\prod^{m}_{i=1}f(y^{(i)}|X^{(i)})=(2\pi \sigma^2)^{-\frac{m}{2}}\exp\{-\frac{(Y-X\beta)^{'}(Y-X\beta)}{2\sigma^2}\}$

$\ln{L(\beta}）=-\frac{m}{2}\ln{2\pi \sigma^2}-\frac{(Y-X\beta)^{'}(Y-X\beta)}{2\sigma^2}$

对 $\beta$ 求导

$\frac{\partial \ln{L(\beta}）}{\partial \beta}=-\frac{1}{2\sigma^2}\frac{\partial tr((Y-X\beta)^{'}(Y-X\beta))}{\partial \beta}=0$

则 $\hat{\beta}=(X^{'}X)^{-1}X^{'}Y$

4.从几何视角看线性回归

对于 $X$ 构成的n维特征空间(每个特征向量构成一个轴)，向量 $Y$ 在空间上的投影为 $X\beta$ ，要使残差 $Y-X\beta$ 最小，则 $Y-X\beta$ 与超平面 $X$ 正交。

即 $X^{'}(Y-X\beta)=0(两个向量垂直，即a^{'}b=0)$

得 $\hat{\beta}=(X^{'}X)^{-1}X^{'}Y$

实现

class myLinReg():
    def __init__(self, eta=1e-1, thread=1e-3, max_iter=1000):
        self._X = None
        self._y = None
        self._theta = None
        self._eta = eta
        self._thread = thread
        self._max_iter = max_iter
        
    def fit(self, X_train, y_train):
        self._X = np.c_[np.ones(X_train.shape[0]), X_train]
        self._y = np.array(y_train).reshape(-1, 1)
        self._theta = np.linalg.inv(self._X.T @ self._X) @ self._X.T @ self._y
    
    # 目标函数
    def _J(self):
        h = self._X @ self._theta
        m = len(self._y)
        return np.sum((self._y - h)**2) / 2*m
    
    # 梯度
    def _dJ(self):
        h = self._X @ self._theta
        m = len(self._y)
        return (self._X.T @ (h - self._y)) / m
    
    def fit_gd(self, X_train, y_train):
        self._X = np.c_[np.ones(X_train.shape[0]), X_train]
        self._y = np.array(y_train).reshape(-1, 1)
        self._theta = np.zeros(self._X.shape[1]).reshape(-1,1)
        count = 1
        while count < self._max_iter:
            old_J = self._J()
            dJ = self._dJ()
            self._theta -= self._eta * dJ
            new_J = self._J()
            if np.abs(new_J - old_J) < self._thread:
                break
            count +=1
        return self
    
    def predict(self, X_test):
        X_b = np.c_[np.ones(X_test.shape[0]), X_test]
        return X_b @ self._theta
    
    def score(self, X_test, y_test):
        y_pred = self.predict(X_test).reshape(-1,1)
        y_test = np.array(y_test).reshape(-1,1)
        sse = np.sum((y_test - y_pred)**2)
        sst = np.sum((y_test - y_test.mean())**2)
        return 1 - sse/sst