线性回归的原理

• 多元回归的一般式

$$\hat{y}(\theta ,x)={\theta }_0+{\theta }_1{x}_1+\dots +{\theta }_p{x}_p$$

其中：
$\hat{y}$:预测值

${\theta }_0,{\theta }_1,\dots {\theta }_p$:参数

$x_0,x_1,\dots x_p$:自变量

• 方程的拟合思路

最小二乘法

• 通过求导求极值解法

下面通过一元线性回归模型做推导：

样本回归模型：

$${\hat{y}}_i={\hat{\beta }}_0+{\hat{\beta }}_1{x}_i$$

残差平方和(损失函数):

$$Q=\sum _{i=1}^n{e}_i^2=\sum _{i=1}^n{\left(y_i-{\hat{y}}_i\right)}^2=\sum _{i=1}^n{\left(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i\right)}^2$$

通过使得Q最小时的参数确定直线，即把它们看作是Q的函数，就变成了一个求极值的问题，可以通过求导数得到。求Q对两个待估参数的偏导数

$$\{\begin{array}{rl}\frac{\partial Q}{\partial {\hat{\beta }}_0}& =2\sum _{i=1}^n\left(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i\right)(-1)=0\\ \frac{\partial Q}{\partial {\hat{\beta }}_1}& =2\sum _{i=1}^n\left(y_i-{\hat{\beta }}_0-{\hat{\beta }}_1{x}_i\right)\left(-x_i\right)=0\end{array}$$

解得：

$${\hat{\beta }}_1=\frac{n\sum x_i{y}_i-\sum x_i\sum y_i}{n\sum x_i^2-{\left(\sum x_i\right)}^2}$$

$${\hat{\beta }}_0=\frac{\sum x_i^2\sum y_i-\sum x_i\sum x_i{y}_i}{n\sum x_i^2-{\left(\sum x_i\right)}^2}$$

其中：

$${\hat{\beta }}_0=\bar{y}-{\hat{\beta }}_1\bar{x}$$

$${\hat{\beta }}_1=\frac{{\sum }_{i=1}^n{x}_i{y}_i-n\overline{xy}}{{\sum }_{i=1}^n{x}_i^2-n{\bar{x}}^2}$$

• 通过矩阵求解

定义：
$$h_{\theta }\left(x_1,x_2,\dots x_{n-1}\right)={\theta }_0+{\theta }_1{x}_1+\dots +{\theta }_n{x}_{n-1}$$
现在有m个样本，每个样本有n−1维特征将所有样本点代入模型中得：

$$\begin{array}{l}{h}_1={\theta }_0+{\theta }_1{x}_{1,1}+{\theta }_2{x}_{1,2}+\dots +{\theta }_{n-1}{x}_{1,n-1}\\ h_2={\theta }_0+{\theta }_1{x}_{2,1}+{\theta }_2{x}_{2,2}+\dots +{\theta }_{n-1}{x}_{2,n-1}\\ ⋮\\ h_m={\theta }_0+{\theta }_1{x}_{m,1}+{\theta }_2{x}_{m,2}+\dots +{\theta }_{n-1}{x}_{m,n-1}\end{array}$$

为方便用矩阵表示令$x_0=1$于是上述方程可以用矩阵表示为：

$$\mathbf{h}=\mathbf{X}\theta$$
其中，h为m x 1的向量, 代表模型的理论值，θ 为n x 1的向量，X为 m x n维的矩阵，m代表样本的个数,n代表样本的特征数，于是目标损失函数用矩阵表示为：

$$J(\theta )=\Vert \mathbf{h}-\mathbf{Y}{\Vert }^2=\Vert \mathbf{X}\theta -\mathbf{Y}{\Vert }^2=(\mathbf{X}\theta -\mathbf{Y}{)}^T(\mathbf{X}\theta -\mathbf{Y})$$

关于矩阵的求导补充以下知识点：

$$\begin{array}{rl}\frac{\partial x^{T}a}{\partial x}& =\frac{\partial a^{T}x}{\partial x}=a\\ \frac{\partial x^{T}Ax}{\partial x}& =Ax+A^{T}x\end{array}$$
又由于如果矩阵A是对称的：

$$Ax+A^{T}x=2Ax$$
对目标函数化简：

$$J(\theta )={\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}Y-Y^{T}X\theta +Y^{T}Y$$
则θ最小值：

$$\begin{array}{rl}\hat{\theta }=\underset{\theta }{\mathrm{argmin}}L(\theta )& ⟶\frac{\partial }{\partial \theta }L(\theta )=0\\ & ⟶2X^{T}X\hat{\theta }-2X^{T}Y=0\\ & ⟶\hat{\theta }={\left(X^{T}X\right)}^{-1}{X}^{T}Y=X^{+}Y\end{array}$$
这个式子中${\left(X^{T}X\right)}^{-1}{X}^T$又被称为伪逆写作$X^{+}$,对于行满秩或者列满秩的X ，可以直接求解。

• 几何意义法求解

求解 b 在A的列向量空间中的投影。

在几何上，最小二乘法相当于模型（这里就是直线）和试验值的距离的平方求和，假设我们的试验样本张成一个 维空间（满秩的情况）模型写成$f(w)=X\beta$而最小二乘法就是说希望 Y 和这个模型距离越小越好，于是它们的差应该与这个张成的空间垂直：

$$X^T\cdot (Y-X\beta )=0⟶\beta ={\left(X^{T}X\right)}^{-1}{X}^{T}Y$$

python实现最小二乘法（通过矩阵求解）

class LinearRegression:
    '''使用Python实现的线性回归.（最小二乘法）'''
    def fit(self,X,y):
        '''根据提提供的训练数据X，对模型进行训练。
        
        Parameters
        ------
        X ：类数组类型。形状：[样本数量，特征数量]
        特征矩阵，用来对模型进行训练
        
        y = 类数组类型，形状：[样本数量]
        
        '''
        # 如果X是数组对象的一部分，而不是完整的对象数据（eg:切片）则无法完成矩阵的转换
        X = np.asmatrix(X.copy())  #这里创建X的拷贝对象，避免转换成矩阵时候失败
        # y 为一维结构（行向量或列向量）没有限制，可以不用拷贝
        # 注意：现在要进行矩阵的运算，需要是二维结构，通过reshape方法进行转换
        y = np.asmatrix(y).reshape(-1,1)    #    (-1,1)意思：把数组转换为一列，依据数据自定义行的矩阵
        # 通过最小二乘公式，求解出最佳的权重值
        self.w_ = (X.T * X ).I * X.T * y
    def predict(self,X):
        '''根据参数传递的样本X，delattr样本数据进行预测。
    
         Parameters
         ------
         X：类数组类型。形状[样本数量，特征数量]
         待预测的样本特征（属性）
       
         Returns
         ------
         result:数组类型
                预测的结果
        '''
        # 将X转换为矩阵 注意：需要对X进行拷贝
        X = np.asmatrix(X.copy())
        result = X * self.w_
        # 将矩阵转换为ndarray数组 ，进行扁平化处理 ，然后返回结果
        return np.array(result).ravel()

为什么会用纯Python来实现呢？这种陶（hao）冶（wu）情（yong）操（chu）的做法肯定也不是我自己闲着蛋疼搞着玩的。。。  明明有更好的科学计算库不用，非要纯手写正面硬刚，除了老师用挂科作为威胁还能有谁？

下面是一大纯手打的堆理推导... 写的逻辑有些混乱，后续有时间再慢慢整理 觉得麻烦的小伙伴就不用看了，反正我代码里也写的有注释，代码稍后传上github后更新连接

---------------------------------------------------------------------------------------------------------------
 

多元线性回归 最小二乘法实现

线性回归包括一元线性回归和多元线性回归，一元的是只有一个x和一个y。多元的是指有多个x和一个y。

对于平面中的这n个点，可以使用无数条曲线来拟合。要求样本回归函数尽可能好地拟合这组值。综合起来看，这条直线处于样本数据的中心位置最合理。 选择最佳拟合曲线的标准可以确定为：使总的拟合误差（即总残差）达到最小。有以下三个标准可以选择：

        （1）用“残差和最小”确定直线位置是一个途径。但很快发现计算“残差和”存在相互抵消的问题。
        （2）用“残差绝对值和最小”确定直线位置也是一个途径。但绝对值的计算比较麻烦。
        （3）最小二乘法的原则是以“残差平方和最小”确定直线位置。用最小二乘法除了计算比较方便外，得到的估计量还具有优良特性。这种方法对异常值非常敏感。

样本x为一个多维向量，x=[x1,x2,…xd]

那么权重W也是与X同维度的一个向量，W=[w1,w2,…wd]

所以Y对x的预测可以写为Y=W ∙ x+b(这里为内积的形式)，但通常将b移到W向量里面去，写为W=[b,w1,w2,…wd]，改变后的W用θ表示，θ=[θ0,θ1,θ2,…θd]

因为W的维度变了所以我们也要给x对应的多加一个维度，此时的x为[x0,x1,…xd]，使x0等于1，让其在矩阵相乘时匹配刚才在W中加入的b

现在Y=W ∙ x+b 就等价为 Y=θ∙x

在机器学习的矩阵运算当中，有个约定俗称的规定，当说某某向量时都指的是一个列向量，所以θ和x都是一个d+1维的列向量[θ0⋮θd]、[x0⋮xd]，同时也都是一个(d+1)x1维的矩阵

这里的维度就对应方程里面的元数，有几维就是几元

所以这里用大写的X表示有m个样本，X=[⋯x1T⋯⋯⋮⋯⋯xmT⋯]=[x01⋯xd1⋮⋯⋮x0m⋯xdm]，（这里的X是一个矩阵，包含了所有的样本点x1,x2,…xm而每一个样本点都是有相同的多维的， [⋯x1T⋯⋯⋮⋯⋯xmT⋯]* [θ0⋮θd]= y0⋮ym=Y，得到一个m*1维的矩阵Y，这个矩阵中的每一个参数代表给定的θ参数下对每一个样本xi的预测值yi，Y是对应样本的标签值，也是一个列向量。

得到Y（预测值）过后，就可以算损失函数了，L=i=1m（yi-yi）2，y是对应样本的真实值

L=（y1-y1）2+（y2-y2）2+…+（ym-ym）2

 =[ y1-y1,⋯,ym-ym] ∙ [y1-y1⋮ym-ym]

 =(y1,⋯,ym-[y1,⋯,ym]) ∙ [y1-y1⋮ym-ym]

 =(x∙θ-Y)T ∙  (x∙θ-Y)

因为x、θ和Y都是一个m+1维的列向量所以第一项(x∙θ-Y)的结果也是一个列项量，所以加了一个转置

得到损失函数L=(x∙θ-Y)T (x∙θ-Y)，通过改变θ的值来最小化L

对θ求偏导，∂L∂θ=∂(θTxT-YT)(xθ-Y)∂θ=θTxTxθ-θTxTY-YTxθ+YTY∂θ，每一项都对θ求偏导然后加起来

 ∂L∂θ=2xTxθ-xTY-xTY，要使损失函数最小，就是偏导为0的点

所以推出xTxθ=xTY，θ=( xTx)-1xTY，这就要求(xTx)是可逆的

最小二乘法与梯度下降法的区别：

实现方法和结果不同：最小二乘法是直接对∆求导找出全局最小，是非迭代法。而梯度下降法是一种迭代法，先给定一个β ，然后向∆下降最快的方向调整β ，在若干次迭代之后找到局部最小。梯度下降法的缺点是到最小点的时候收敛速度变慢，并且对初始点的选择极为敏感，其改进大多是在这两方面下功夫。