线性回归

线性回归

一元线性回归

假设对于观测对象x和y我们收集到了一批数据$x = \begin{Bmatrix}x_1,x_2,x_3,\dots,x_n\end{Bmatrix}$ ，$y = \begin{Bmatrix}{y_1,y_2,y_3,\dots y_n}\end{Bmatrix}$

我们希望找到一个一元线性函数(一个因变量y和一个自变量x)

$$y_i = f(x_i) = wx_i + b \\ $$

使得函数（模型）预测出来的值和原来的值的误差平方和

$$S = \sum^{n}_{i=1}{(f(x_i)-y_i)^2}$$

最小，也即是它们的欧式距离最小，这样就会有$f(x_i) \approx y_i$

我们定义代价函数为（这里添加系数$\frac{1}{2}$是为了方便求导）

$$L(w,b)=\frac{1}{2}S =\frac{1}{2}\sum^{n}_{i=1}{(f(x_i)-y_i)^2}$$

所以问题变成了，寻找出$w$和$b$使得$L$最小，也即有$f(x_i)\approx y_i$，达到我们预测（拟合）的目的

即

$$\begin{align} (w^*,b^*) &= \underset{<w,b>}{\operatorname{arg min}} \frac{1}{2}\sum^{n}_{i=1}{(f(x_i)-y_i)^2}\\ &= \underset{<w,b>}{\operatorname{arg min}} \frac{1}{2}\sum^{n}_{i=1}{(wx_i+b-y_i)^2}\\ \end{align}$$

因为代价函数$L(w,b)$是一个凸函数，当它关于w和b的偏导都为0时，则可以取得$w$和$b$的最优解

令$L(w,b)$对$w$和$b$求偏导，得

$$\frac{\partial{L(w,b)}}{\partial w} = \sum^{n}_{i=1}{(wx_i+b-y_i)}x_i \\\frac{\partial{L(w,b)}}{\partial b} = \sum^{n}_{i=1}{(wx_i+b-y_i)} \\ $$

梯度下降法

梯度下降就是通过迭代，不断让函数的参数向着梯度下降的方向走一点点，不断的逼近最优解

设更新步长为$\alpha$，则有更新公式

$$w \leftarrow w - \alpha \frac{\partial{L}}{w} \\ b \leftarrow b - \alpha \frac{\partial{L}}{b} \\ $$

直接求解

我们也可以直接算出它的闭式解（解析解）。令上面两个偏导数等于0，就得到

$$w = \frac{\sum^{n}_{i=1}y_i(x_i-\overline{x})}{\sum{x^2_i-\frac{1}{n}(\sum^{n}_{i=1}wx_i)}} \\ b = \frac{1}{n}\sum_{i=1}^n(y_i-wx_i)$$

多元线性回归

多元线性回归就是具有多个自变量和一个因变量的回归模型，假设自变量x有m个特征，我们对x和y进行了n次观测，则有模型

$$f(x) = w_1x_1 + w_2x_2 + w_3x_3 +\dots + w_mx_m + b$$

把$yx_i$和$w$写成向量的形式(这里$x_i$代表对$x$的第$i$次观测得到的数据)

$$y = \begin{bmatrix} y_1\\ y_2\\ \dots\\ y_n\\ \end{bmatrix} , x_i = \begin{bmatrix} x_1\\ x_2\\ \dots\\ x_n\\ \end{bmatrix} , w = \begin{bmatrix} w_1\\ w_2\\ \dots\\ w_n\\ \end{bmatrix}$$

那我们可以把这个方程写成向量方程的形式

$$f(x_i) = w^Tx_i + b\\ $$

进一步的，对于所有的数据，有数据矩阵$X$

$$X = \begin{bmatrix} x_{11} \, x_{12} \, \dots \, x_{1m}\\ x_{21} \, x_{22} \, \dots \, x_{2m}\\ \dots\\ x_{n1} \, x_{n2} \, \dots \, x_{nm}\\ \end{bmatrix}$$

其中，每一行是一次观测，每一列是一个维度（特征）

然后，为了方便，再把常数项$b$纳入$w$中，并且在$X$中多加一列1

$$\hat{w} = \begin{bmatrix} w_1\\ w_2\\ \dots\\ w_n\\ b\\ \end{bmatrix} , X = \begin{bmatrix} x_{11} \,\, x_{12} \,\, \dots \,\, x_{1m} \,\, 1\\ x_{21} \,\, x_{22} \,\, \dots \,\, x_{2m} \,\, 1\\ \dots\\ x_{n1} \,\, x_{n2} \,\, \dots \,\, x_{nm} \,\, 1\\ \end{bmatrix}$$

则有矩阵方程

$$y = f(X) = X\hat{w}$$

则我们的优化目标就是

$$\hat{w}^* = \underset{<\hat{w}>}{\operatorname{arg min}} (y-X\hat{w})^T(y-X\hat{w})\\ $$

令

$$L(\hat{w}) = (y-X\hat{w})^T(y-X\hat{w}) = (y-X\hat{w})^2$$

对$\hat{w}$求偏导得

$$\frac{\partial L(\hat{w})}{\partial \hat{w}} = 2X^T(X\hat{w}-y)$$

我们的目标就是让$\frac{\partial L(\hat{w})}{\partial \hat{w}}=0$

梯度下降法

像一元线性回归那样，有

$$\hat{w} \leftarrow \hat{w} - \alpha \frac{\partial{L}}{\partial\hat{w}} \\ $$

正规方程法

令

$$2X^T(X\hat{w}-y) = 0$$

则当$X^TX$是正定或者满秩的时候，方程有唯一解（因为互不共线的向量只能找到唯一一个线性组合使其等0）

$$\hat{w} = (X^TX)^{-1}X^Ty$$

其实还有很多的，但是我很懒，不想写了

其实线性回归不单只可以用来拟合线性模型，还可以用来拟合多项式函数、对数函数、指数函数等等，只要通过一定的变换，把原来的问题转换成线性的问题就可以求解，本质上还是在优化一个凸函数，一个最小二乘的问题，其实也不一定是最小二乘，也可以用其他的，比如说误差绝对值，但这种东西是视情况而论的，就这样吧。

编程实现

理论理清楚了，编程就不会太难

# -*- coding: utf-8 -*-
"""
Created on Thu Jan 14 12:53:36 2021

@author: koneko
"""

'''
多元线性回归
'''

import numpy as np
import matplotlib.pyplot as plt

def lossf(X,w,y):
    return np.sum((y-np.dot(X,w))**2)

def init(X,y):
    if X.ndim == 1:
        X = X.reshape(X.size,1)
    if y.ndim == 1:
        y = y.reshape(y.size,1)
    #在x后面多加一列1
    X = np.c_[X,np.ones([X.shape[0],1])]
    n,m = X.shape
    w =  w = np.random.normal(1,0.1,m)
    w = w.reshape(w.size,1)
    return X,y,w

'''
使用正规方程来求
'''
def LRWithNormalEquation(x,y):
    X,y,w = init(x,y)
    inv = np.linalg.inv(np.dot(X.T,X))
    R = np.dot(X.T,y)
    w = np.dot(inv,R)
    return w

'''
通过迭代的方法来求
'''
def LRWithGradientDesc(x,y):
    #初始化
    X,y,w = init(x,y)
 
    delta = 0.001  #收敛系数
    alpha = 0.001  #学习速率
    max_step = 10000 #最大次数
    gradient = 1000
    err = 1000
    loss = []
    i = 1
    while err>delta and i < max_step:
        i += 1
        gradient = 2*np.dot(X.T,(np.dot(X,w)-y))
        w = w - alpha*gradient
        err = lossf(X,w,y)
        loss.append(err)
   
    plt.plot(loss)
    return w
   

def f(X,w):
    return np.dot(X,w)
    
x = np.array([0.50,0.75,1.00,1.25,1.50,1.75,1.75,2.00,
     2.25,2.50,2.75,3.00,3.25,3.50,4.00,4.25,4.50,4.75,5.00,5.50])

y = np.array([10,  26,  23,  43,  20,  22,  43,  50,  62, 50,  55,  75,  
     62,  78,  87,  76,  64,  85,  90,  98])
     
w1 = LRWithGradientDesc(x,y)
w2 = LRWithNormalEquation(x,y)

X,y,w = init(x,y)
y1 = f(X,w1)
y2 = f(X,w2)
plt.subplot(1,2,1)
plt.title('GradientDesc')
plt.scatter(x,y)
plt.plot(x,y1)
plt.subplot(1,2,2)
plt.scatter(x,y)
plt.plot(x,y2)
plt.title('NormalEquat')

参考资料