线性回归

线性回归

一元线性回归

假设对于观测对象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')

参考资料

posted @ 2021-01-03 21:40  裏表異体  阅读(299)  评论(0编辑  收藏  举报