机器学习(三)-多元线性回归(数学推导及代码实现)

前面讨论了 y = ax + b 考虑的只有一个 特征值(因素)的情况下,但在很多情况下 特征值不只有一个 打个比方 要预测房价 要考虑的不只是面积 还要有 地段 建造年代 户型 等等 ,此时就要用到多元线性回归了。

$({\theta }_0,{\theta }_1,{\theta }_2,{\theta }_3,.....,{\theta }_n)$ $\theta$代表一系列我们需要学习出来的参数
$(1,X_1^{(i)},X_2^{(i)}，X_3^{(i)},.....,X_n^{(i)}$ X代表了要训练的参数(特征) 比如 面积 朝向 户型 装修 交通 等等特征
决定一间房屋的价格可以由很多因素综合的出 这一组综合的权重就是要求解得出的。
房屋价格 = 面积*面积的权重 + 朝向 * 朝向的权重 + 户型 * 户型的权重 + 装修 * 装修的权重 + 交通 * 交通的权重
${\hat{y}}^{(i)}={\theta }_0{X}_0^{(i)}+{\theta }_1{X}_1^{(i)}+{\theta }_2{X}_2^{(i)}+{\theta }_3{X}_3^{(i)},.....,{\theta }_n{X}_n^{(i)}$

此处${\theta }_0{X}_0^{(i)}$ =b 为偏置(截距)

${\hat{y}}^{(i)}=b+{\theta }_1{X}_1^{(i)}+{\theta }_2{X}_2^{(i)}+{\theta }_3{X}_3^{(i)},.....,{\theta }_n{X}_n^{(i)}$

使用和一元线性同样的 损失函数:
$argmin{\sum }_{i=1}^m(y^{(i)}-{\hat{y}}^{(i)}{)}^2$

写成矩阵形式
$\theta =({\theta }_0,{\theta }_1,{\theta }_2,{\theta }_3,.....,{\theta }_n{)}^T$
$X^{(i)}=(1,X_1^{(i)},X_2^{(i)},X_3^{(i)},.....,X_n^{(i)}{)}^T$
方程于是就变成:

${\hat{y}}^{(i)}=\theta \cdot X^{(i)}$
相当于 $y=\theta X+b$ 的矩阵形式

输入的样本:
$X=\left[\begin{array}{ccccc}1& X_{(1)}^{(1)}& X_{(2)}^{(1)}& ....& X_{(n)}^{(1)}\\ 1& X_{(1)}^{(2)}& X_{(2)}^{(2)}& ....& X_{(n)}^{(2)}\\ .....& & & & \\ 1& X_{(1)}^{(m)}& X_{(2)}^{(m)}& ....& X_{(n)}^{(m)}\end{array}\right]$

待求得权重:$\theta ={\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ {\theta }_2\\ ....\\ {\theta }_n\end{array}\right]}^T$ 真实值:$y=\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ....\\ y_n\end{array}\right]$

$X为m\ast (n+1)的矩阵$
$\theta 为(n+1)\ast 1列的向量$
$y为(n+1)\ast 1列的向量$

下标(n) 代表特征(参数)数量
上表(m)代表 训练数据数量
X 第一列 为1 相当于 $\theta$为 偏置 b

$argmin{\sum }_{i=1}^m(y^{(i)}-{\hat{y}}^{(i)}{)}^2$ → $argmin(Y-X\cdot \theta {)}^2$ →$argmin(Y-X\cdot \theta {)}^T\cdot (Y-X\cdot \theta )$

所以 问题就转化成 求上述式子的最小值,求最小值是一个凸优化问题 遇到这种问题 一般 先需要证明 这个式子是有凸函数连续可导。

(以下推导可以不看 直接看结论)
一般推导思路

凸集定义

设集合$D\subset R^n$,如果对集合内任意元素$x,y\subset D$与任意的$a\subset [0,1]$,有$ax+(1-a)y\subset D$,则称集合D是凸集。

任意连2条直线 所形成的集合还在原空间中则称为凸集。

梯度定义

设n原函数$f(x)$ 对自变量$x=(x_1,x_2,...,x_n{)}^T$的个分量$x_i$的偏导数$\frac{\partial f(x)}{\partial x_i}i=(1,2,3...,n)$都存在,则函数 f(x)在x处一阶可导,并称向量
$$\nabla f(x)=\left[\begin{array}{c}\frac{\partial f(x)}{\partial x_1}\\ \frac{\partial f(x)}{\partial x_2}\\ ...\\ \frac{\partial f(x)}{\partial x_n}\end{array}\right]$$
为函数$f(x)$在x处的一节阶倒数或梯度,记为$\nabla f(x)$(列向量)

海塞矩阵

Hession(海塞)矩阵定义:设n元函数$f(x)$对自变量$x=(x_1,x_2,....,x_n{)}^T$的各分量$x_i$的二阶偏导数$\frac{{\partial }^{2}f(x)}{\partial x_i\partial x_j}(i,j=1,2,3.....,n)$都存在,则称函数$f(x)$在点$x$处二阶可导,并称矩阵
$${\nabla }^{2}f(x)=\nabla f(x)=\left[\begin{array}{cccc}\frac{{\partial }^{2}f(x)}{\partial x_1^2}& \frac{{\partial }^{2}f(x)}{\partial x_1\partial x_2}& ...& \frac{{\partial }^{2}f(x)}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_2^2}& ...& \frac{{\partial }^{2}f(x)}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f(x)}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f(x)}{\partial x_n\partial x_2}& ...& \frac{{\partial }^{2}f(x)}{\partial x_n\partial x_n}\end{array}\right]$$
为$f(x)$在$x$处的二阶倒数或 Hession 矩阵 ,记为${\nabla }^{2}f(x)$,若$f(x)$对$x$ 各变元的所有二阶偏导数都连续,则$\frac{{\partial }^{2}f(x)}{\partial x_i{x}_j}$ = $\frac{{\partial }^{2}f(x)}{\partial x_j{x}_i}$

多元实值函数凹凸型判定定理

设$D\subset R^n$是非空开凸集,$f:D\subset R^n,$且$f(x)$在$D$上二阶连续可微,如果$f(x)$的Hession矩阵${\nabla }^{2}f(x)$在$D$上是正定的,则$f(x)$是D上的严格凸函数。
凸充分性定理
若$f:R^n → R $ 是凸函数,且$f(x)$一阶连续可微,则 $x^{\ast }$是全局最优解(全局最小值)的充分必要条件是$\nabla f(x^{\ast })=0$,其中$f(x)$关于$x$的一阶导数(也称梯度)

[标量-向量]的矩阵微分公式为:
$$\frac{\partial y}{\partial x}=\left(\begin{array}{c}\frac{\partial y}{\partial x_1}\\ \frac{\partial y}{\partial x_2}\\ ⋮\\ \frac{\partial y}{\partial x_n}\end{array}\right)$$
$$(分母布局)$$
$$\frac{\partial y}{\partial x}=\left(\begin{array}{c}\frac{\partial y}{\partial x_1}\frac{\partial y}{\partial x_2}\dots \frac{\partial y}{\partial x_n}\end{array}\right)$$
$$(分子布局)$$
其中,$x=(x_1,x_2,.....,x_n{)}^T$为n为向量,$y$为 $x$的$n$元标量函数.
分子分母布局按照习惯选择。
由[标量 -向量]的矩阵微分公式可推得:
$\frac{\partial x^{T}a}{\partial x}=\frac{\partial a^{T}x}{\partial x}$ = $\frac{\partial y}{\partial x}=\left(\begin{array}{c}\frac{\partial (a_1{x}_1+a_2{x}_2+...+a_n{x}_n)}{\partial x_1}\\ \frac{\partial (a_1{x}_1+a_2{x}_2+...+a_n{x}_n)}{\partial x_2}\\ ⋮\\ \frac{\partial (a_1{x}_1+a_2{x}_2+...+a_n{x}_n)}{\partial x_n}\end{array}\right)$=$\left(\begin{array}{c}{a}_1\\ a_2\\ ⋮\\ a_n\end{array}\right)$=a
同理可推得:$\frac{\partial x^{T}Bx}{\partial x}=(B+B^T)x$

有了以上定理就可以开始证明了:

1. 证明损失函数$E\theta$ 是关于$\theta$的凸函数。

$\frac{\partial E\theta }{\partial \theta }=\frac{\partial ((Y-X\cdot \theta {)}^T(Y-X\cdot \theta ))}{\partial \theta }$

=$\frac{\partial }{\partial \theta }(Y-{\theta }^T{X}^T)(Y-X\theta )$

=$\frac{\partial }{\partial \theta }[-Y^{T}X\theta -{\theta }^T{X}^{T}Y+{\theta }^T{X}^{T}X\theta ]$

=$\frac{\partial Y^{T}X\theta }{\partial \theta }-\frac{\partial {\theta }^T{X}^{T}Y}{\partial \theta }+\frac{\partial {\theta }^T{X}^{T}X\theta }{\partial \theta }$

由$\frac{\partial x^{T}a}{\partial x}=\frac{\partial a^{T}x}{\partial x}$ = $\frac{\partial y}{\partial x}=a$ $\frac{\partial x^{T}Bx}{\partial x}=(B+B^T)x$可得:

$\frac{\partial E_{\theta }}{\partial \theta }=-X^{T}y-X^{T}y+(X^{T}X+X^{T}X\theta )$

=$2X^T(X\theta -Y)$ 为一阶偏导

再对一阶偏导数求二阶偏导数:
$\frac{{\partial }^{2}E\theta }{\partial \theta \partial {\theta }^T}$ = $\frac{\partial }{\partial \theta }(\frac{\partial E\theta }{\partial \theta })$ = $\frac{\partial }{\partial E\theta }(2X^{T}X\theta -2X^{T}Y)$=$2X^{T}X$ (即Hession矩阵)

假设$X^{T}X$为正定矩阵($E\theta$ 是关于$\theta$的凸函数)
令 $E\theta$ 关于$\theta$的一阶导数为0 前面已经证明过了:
$\frac{\partial E\theta }{\partial \theta }$=$2X^T(X\theta -Y)=0$
$2X^{T}X\theta -2X^{T}Y=0$
$2X^{T}X\theta =2X^{T}Y$
2边除以2:
$X^{T}X\theta =X^{T}Y$
2边乘以 $(X^{T}X{)}^{-1}$
$(X^{T}X{)}^{-1}{X}^{T}X\theta =(X^{T}X{)}^{-1}{X}^{T}Y$
得:$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$

以上证明毕：
$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$

问题:时间复杂度高:O($n^3$)(优化后O($n^{2.4}$)
优点:不需要对数据做归一化处理

import numpy as np
class LinearRegression:
    def __init__(self):
        """初始化linear Regression模型"""
        self.coef_ = None
        self.interception_ = None
        self._theta = None
    def fit_normal(self,X_train,y_train):
        assert X_train.shape[0] == y_train.shape[0], \
            "the size of x_train must be equal to the size of y_train"
        #传入的训练样本没有截距 在每一行加一个截距项
        X_b = np.hstack([np.ones((len(X_train),1)),X_train])
        #根据公式计算 theta 的值
        self._theta = np .linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train)
        self.interception_=self._theta[0]
        self.coef_ = self._theta[1:]
        return self
    def predict(self,X_predict):
        assert self.coef_ is not None and  self.interception_ is not None, \
            "you must run fit_normal before predict"
        assert X_predict.shape[1] == len(self.coef_), \
            "the size of X_predict must be equal to the size of X_train"
        X_predict_b = np.hstack([np.ones((len(X_predict),1)), X_predict])
        return X_predict_b.dot(self._theta)

    def score(self, X_test,y_test):
        y_predict = self.predict(X_test)
        return r_Squared(y_test,y_predict)

    def __repr__(self):
        return "LinearRegression()"
	def r_Squared(y_true,y_predict):
	    assert len(y_true)  == len(y_predict), \
	        "the size of y_true must be equal to the size of y_predict"
	    return 1- (mean_squared_error(y_true,y_predict) / np.var(y_true))

测试

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from ml_utils.data_split import train_test_split
from sklearn import datasets
boston = datasets.load_boston()
x = boston.data
y = boston.target
reg = LinearRegression()
reg.fit_normal(x_train,y_train)
print(reg.coef_)
print(reg.interception_)
print(reg.score(x_test,y_test))

Sklearn中的线性回归

from sklearn.linear_model import LinearRegression
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
boston = datasets.load_boston()
x =boston.data
y = boston.target
x_train,x_test,y_train,y_test = train_test_split(x,y,random_state=666);
lin_reg = LinearRegression()
lin_reg.fit(x_train,y_train)
#查看训练的权重参数
lin_reg.coef_
#截距
lin_reg.intercept_
#评分
lin_reg.score(x_test,y_test)

通过线性模型能找到 特征 比如 RM 是权重最高的 说明房间数量 和 房价正相关 而且 权重较高 CHAS 第二高的权重 表示 波士顿房子邻河不临河 临河和房价 也是正相关 而排在最后一个的参数 NOX 表示一氧化氮 的含量 和房价负相关 。这也说明线性回归法 对数据有可解释性 不管模型预测好坏 先用线性模型进行预测 能直观的观测处数据中存在的规律。