线性回归

线性回归

一元线性回归

假设对于观测对象x和y我们收集到了一批数据(x = egin{Bmatrix}x_1,x_2,x_3,dots,x_nend{Bmatrix}) ，(y = egin{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)，达到我们预测（拟合）的目的

即

[egin{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 = egin{bmatrix} y_1\ y_2\ dots\ y_n\ end{bmatrix} , x_i = egin{bmatrix} x_1\ x_2\ dots\ x_n\ end{bmatrix} , w = egin{bmatrix} w_1\ w_2\ dots\ w_n\ end{bmatrix} ]

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

[f(x_i) = w^Tx_i + b\ ]

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

[X = egin{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} = egin{bmatrix} w_1\ w_2\ dots\ w_n\ b\ end{bmatrix} , X = egin{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) = Xhat{w} ]

则我们的优化目标就是

[hat{w}^* = underset{<hat{w}>}{operatorname{arg min}} (y-Xhat{w})^T(y-Xhat{w})\ ]

令

[L(hat{w}) = (y-Xhat{w})^T(y-Xhat{w}) = (y-Xhat{w})^2 ]

对(hat{w})求偏导得

[frac{partial L(hat{w})}{partial hat{w}} = 2X^T(Xhat{w}-y) ]

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

梯度下降法

像一元线性回归那样，有

[hat{w} leftarrow hat{w} - alpha frac{partial{L}}{partialhat{w}} \ ]

正规方程法

令

[2X^T(Xhat{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')

参考资料