线性回归

一、形式化定义

在统计学中，线性回归是一种线性方法，用于建模标量响应与一个或多个解释变量之间的关系。

当只有一个解释变量时，我们称之为简单线性回归：

$\begin{aligned} &数据假设：y_i\in R\\ &模型假设：y_i=x_i^Tw+\epsilon_i\\ &\epsilon_i\sim N(0,\sigma^2),y_i|x_i\sim N(x_i^Tw,\sigma^2) \end{aligned}$

根据对分布的模型假设，我们最终得到的条件概率为：

$P(y_i|x_i,w)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x_i^Tw-y_i)^2}{2\sigma^2}}$

我们假设模型是一条过原点的直线（类似于感知机通过升维即可使得函数过原点）而对于每个特征$x_i$，它的标签$y_i$则从一个平均值为$w^Tx_i$，方差为$σ^2$的高斯分布中得出，我们在线性回归建模中的任务便是根据数据集估计斜率$w$。

单标量预测变量$x$和单标量响应变量$y$的最简单情况称为简单线性回归，而多个预测变量$x$参与的预测为多元线性回归，几乎所有现实世界的回归模型都涉及多个预测因子，线性回归的基本描述通常用多元回归模型来表述。

二、参数估计

线性回归模型中的重要参数就是“斜率”$w$，我们同样可以通过$MLE$和$MAP$两种方式对它进行估计。

进行下列推导之前，我们先了解几个常用的矩阵求导公式：

$\frac{\partial Ax}{\partial x}=A^T,\frac{\partial x^TA}{\partial x}=A,\frac{\partial A^TxB}{\partial x}=AB^T\\\frac{\partial A^Tx^TB}{\partial x}=BA^T,\frac{\partial A^Txx^TB}{\partial x}=(AB^T+BA^T)x,\frac{\partial x^TAx}{\partial x}=2Ax$

再了解一些矩阵的相关性质：

$(A^T)^T=A,(A+B)^T=A^T+B^T,(AB)^T=B^TA^T$

2.1 极大似然估计MLE

$\begin{aligned} &\hat{w}_{MLE}=\underset{w}{argmax}~P(Y,X|w)=\underset{w}{argmax}~\prod_{i=1}^nP(y_i,x_i|w)\\ &~~~~~~~~~~~=\underset{w}{argmax}~\prod_{i=1}^nP(y_i|x_i,w)P(x_i|w)=\underset{w}{argmax}~\prod_{i=1}^nP(y_i|x_i,w)\\ &~~~~~~~~~~~=\underset{w}{argmax}~\sum_{i=1}^n\log P(y_i|x_i,w)=\underset{w}{argmax}~\sum_{i=1}^n\log\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(w^Tx_i-y_i)^2}{2\sigma^2}}\\ &~~~~~~~~~~~=\underset{w}{argmin}~\sum_{i=1}^n(x_i^Tw-y_i)^2=\underset{w}{argmin}~\frac1n\sum_{i=1}^n(x_i^Tw-y_i)^2 \end{aligned}$

我们可以发现，线性回归问题的极大似然估计结果与最小二乘法（OLS）的形式完全相同。

此外我们还可以推导出$w$的闭合形式解，这样会使得计算更加方便：根据上述估计值，我们可以定义损失函数$l(w)$

$\begin{aligned} &l(w)=(Xw-Y)^2 \end{aligned}$

我们要求解$l(w)'=0$时所对应的$w$，$X$为$n\times d$维矩阵，$w$为$d\times 1$维矩阵，$Y$为$n\times1$维矩阵：

$\begin{aligned} &l(w)=(Xw-Y)^T(Xw-Y)=(w^TX^T-Y^T)(Xw-Y)\\ &~~~~~~~=w^TX^TXw-w^TX^TY-Y^TXw-Y^TY\\ &\frac{\partial l(w)}{\partial w}=2X^TXw-2X^TY=0\\ &则有：w=(X^TX)^{-1}X^TY \end{aligned}$

2.2 极大后验估计MAP

引入一个附加的模型假设：

$P(w)=\frac1{\sqrt{2\pi\tau^2}}e^{-\frac{w^Tw}{2\tau^2}}$

则极大后验估计的过程为：

$\begin{aligned} &\hat{w}_{MAP}=\underset{w}{argmax}~P(w|Y,X)=\underset{w}{argmax}~\frac{P(Y,X|w)P(w)}{P(Y,X)}\\ &~~~~~~~~~~~=\underset{w}{argmax}~P(Y,X|w)P(w)=\underset{w}{argmax}~\prod_{i=1}^nP(y_i,x_i|w)P(w)\\ &~~~~~~~~~~~=\underset{w}{argmax}~\sum_{i=1}^n\log P(y_i,x_i|w)+\log P(w)\\ &~~~~~~~~~~~=\underset{w}{argmax}~-\sum_{i=1}^n(x_i^Tw-y_i)^2-\frac{1}{2\tau^2}w^Tw\\ &~~~~~~~~~~~=\underset{w}{argmin}~\sum_{i=1}^n(x_i^Tw-y_i)^2+\lambda||w||_2^2 \end{aligned}$

我们将这一回归称为岭回归，它同样存在闭合形式：我们定义损失函数为：

$l(w)=(Xw-Y)^2+\lambda w^Tw$

我们要求解$l(w)'=0$时所对应的$w$：

$\begin{aligned} &\frac{\partial l(w)}{\partial w}=2X^TXw-2X^TY+2\lambda Iw=0\\ &则有：w=(X^TX+\lambda I)^{-1}X^TY \end{aligned}$

三、模型实现

对于简单的线性回归模型，我们通过闭合形式可以直接求解，将分别通过手动实现和调库的方式构造模型。

3.1 数据集

我们使用波士顿房价数据集，查看数据集的内容并对其进行存储依靠以下部分代码：

'''屏蔽warning'''
import warnings
warnings.filterwarnings("ignore")
'''导入重要的库'''
import numpy as np
import pandas as pd
from sklearn.datasets import load_boston

boston=load_boston()
print(boston.data)#输出数据集内容
print(boston.target)#输出每个数据集的标签
print(boston.feature_names)#输出数据集每一列对应的特征
print(len(boston.data))#输出样本数量
print(len(boston.feature_names))#输出特征数量
'''将数据集存储为csv文件'''
df=pd.DataFrame(data=boston.data,columns=boston.feature_names)#构建二维表格
df['label']=boston.target
df.to_csv("data/boston_data.csv",index=False)

波士顿房价数据集有506个样本，每个样本有13个特征，接下来我们利用线性回归模型对数据集进行拟合。

3.2 手动实现模型

我们的评估指标选择$r2$系数，其定义如下：

$R^2=1-\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{\sum_{i=1}^n(y_i-\overline{y})^2}=1-\frac{MSE}{Var}$

显然，预测越准$r2$系数越趋近于1。

'''屏蔽warning'''
import warnings
warnings.filterwarnings("ignore")
'''导入重要的库'''
import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score

class LinearRegression(object):
    '''初始化'''
    def __init__(self,b,lam):
        self.b=b        #偏差项
        self.lam=lam    #岭回归所要用到的岭系数
    '''求解闭合形式'''
    def fit_OLS(self,X,Y):  	#最小二乘法的闭合形式
        self.w=np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)),X.T),Y)
    def fit_Ridge(self,X,Y):    #岭回归的闭合形式
        self.w=np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)+self.lam*np.eye(len(X[0]))),X.T),Y)
    '''预测结果'''
    def predict(self,X_test):
        return np.dot(X_test,self.w)

'''载入数据'''
boston=load_boston()
X=boston.data
Y=boston.target
'''数据集划分'''
X_train,X_test,y_train,y_test=train_test_split(X,Y,test_size=0.2)#选取20%的数据作为测试集
'''模型拟合'''
linear_regression=LinearRegression(0,1)
# linear_regression.fit_OLS(X_train,y_train)
linear_regression.fit_Ridge(X_train,y_train)
'''模型预测'''
y_pred=linear_regression.predict(X_test)
print(r2_score(y_test,y_pred)) #0.7797762734166735
for i in range(len(y_test)):
    print("true:",y_test[i],"pred:",y_pred[i])

3.3 调库实现模型

调库主要用到的是sklearn中的LinearRegression模型，其具体参数请查阅文档

'''屏蔽warning'''
import warnings
warnings.filterwarnings("ignore")
'''导入重要的库'''
import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.linear_model import LinearRegression


'''载入数据'''
boston=load_boston()
X=boston.data
Y=boston.target
'''数据集划分'''
X_train,X_test,y_train,y_test=train_test_split(X,Y,test_size=0.2)#选取20%的数据作为测试集
'''模型拟合'''
linear_regression=LinearRegression()
linear_regression.fit(X_train,y_train)
'''模型预测'''
y_pred=linear_regression.predict(X_test)
print(r2_score(y_test,y_pred))	#0.748825851468819
for i in range(len(y_test)):
    print("true:",y_test[i],"pred:",y_pred[i])