多元线性回归【正规方程/sklearn】
- 1. 基本概念
- 1.1 线性回归
- 1.2 一元简单线性回归
- 1.3 最优解
- 1.4 多元线性回归
- 2. 正规方程求最优解
- 2.1 线性回归的损失函数（最小二乘法）
- 2.2 推导正规方程
- 2.3 正规方程练习
- 2.4 使用sklearn计算多元线性方程
- 2.5 凸函数
- 3. 线性回归算法推导
- 3.1 深入理解回归
- 3.2 误差分析
- 3.3 极大似然估计
- 3.4 高斯分布——概率密度函数
- 3.5 误差总似然
- 3.6 最小二乘法MSE
- 3.7 归纳总结
- 4. 线性回归实战
- 4.1 使用正规方程求解
- 4.1.1 简单线性回归
- 4.1.2 多元线性回归
- 4.2 sklearn求解
- 4.2.1 简单线性回归
- 4.2.2 多元线性回归
1. 基本概念
1.1 线性回归
线性回归是机器学习中有监督机器学习下的一种算法。 回归问题主要关注的是因变量（需要预测的值，可以是一个也可以是多个）和一个或多个数值型的自变量（预测变量）之间的关系。
-
需要预测的值：即目标变量，
target， y y y，连续值。 -
影响目标变量的因素： x 1 x_1 x1… x n x_n xn，可以是连续值也可以是离散值。
-
因变量和自变量之间的关系：即模型
model，是我们要求解的。
1.2 一元简单线性回归
1. 概念
- 前面提到过，算法说白了就是公式，简单线性回归属于一个算法，它所对应的公式：
y = w x + b y = wx + b y=wx+b
-
这个公式中， y y y 是目标变量即未来要预测的值， x x x 是影响 y y y 的因素， w , b w,b w,b 是公式上的参数即要求的模型。其实 b b b 就是咱们的截距， w w w 就是斜率，所以很明显如果模型求出来了，未来影响 y y y 值的就只有 x x x，也可以说影响 y y y 值的因素只有一个，所以称为一元，简单线性回归。
-
同时可以发现从 x x x 到 y y y 的计算， x x x 只是一次方，所以该算法叫线性回归。
2. 代码示例
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 10, 50)
y = 0.8 * x - 5
plt.plot(x, y, color='red')
1.3 最优解
1. 一些概念
-
Actual value：真实值，一般使用 y y y 表示。 -
Predicted value：预测值，是把已知的 x x x 带入到公式里面和学习出来的参数 w , b w,b w,b 计算得到的，一般使用 y ^ \hat{y} y^ 表示。 -
Error：误差，预测值和真实值的差距，一般使用 ε \varepsilon ε 表示。 -
Loss：整体的误差，Loss通过损失函数Loss function计算得到。 -
最优解：不断优化参数，尽可能的找到一个模型使得整体的误差
Loss最小。
2. 求解最优解的具体做法
- 通过一些具体的优化方法，找到
Loss function图像中的最低点或局部最低点。
1.4 多元线性回归
1. 一般表示
- 现实生活中，往往影响结果 y y y 的因素不止一个，这时 x x x 就从一个变成了 n 个， x 1 x_1 x1… x n x_n xn 同时简单线性回归的公式也就不在适用了。多元线性回归公式如下：
y ^ = w 1 x 1 + w 2 x 2 + … … + w n x n + b \hat{y} = w_1x_1 + w_2x_2 + …… + w_nx_n + b y^=w1x1+w2x2+……+wnxn+b
- b b b 是截距，也可以使用 w 0 w_0 w0来表示：
y ^ = w 1 x 1 + w 2 x 2 + … … + w n x n + w 0 \hat{y} = w_1x_1 + w_2x_2 + …… + w_nx_n + w_0 y^=w1x1+w2x2+……+wnxn+w0
y ^ = w 1 x 1 + w 2 x 2 + … … + w n x n + w 0 ∗ 1 \hat{y} = w_1x_1 + w_2x_2 + …… + w_nx_n + w_0 * 1 y^=w1x1+w2x2+……+wnxn+w0∗1
- 其中， b b b 可以看作是变量 x 0 x_0 x0恒为1的 w 0 ∗ x 0 w_0*x_0 w0∗x0。
2. 向量表示
- X X X 表示所有的变量，是一维向量（包含 x 0 x_0 x0，恒为1）； W W W 表示所有的系数（包含 w 0 w_0 w0），是一维向量，根据向量乘法规律，一般形式可改写为：
y ^ = W T X \hat{y} = W^TX y^=WTX
2. 正规方程求最优解
2.1 线性回归的损失函数（最小二乘法）
1. 线性回归的损失函数
J ( θ ) = 1 2 ∑ i = 0 n ( y i − y i ^ ) 2 J(\theta) = \frac{1}{2}\sum\limits_{i = 0}^n(y_i - \hat{y_i})^2 J(θ)=21i=0∑n(yi−yi^)2
- n n n：样本数量；
- y i y_i yi：第 i i i 个样本的真实值；
- y i ^ \hat{y_i} yi^：第 i i i 个样本的测量值。
此处的 1 2 \frac{1}{2} 21没有实际意义，只是为了后续求导方便计算，科学家们总结的。
2. 使用矩阵表示
J ( θ ) = 1 2 ( X θ − y ) T ( X θ − y ) J(\theta) = \frac{1}{2}(X\theta - y)^T(X\theta - y) J(θ)=21(Xθ−y)T(Xθ−y)
- X X X是总体样本；
- θ \theta θ是参数向量，也就是之前提到的 W W W；
- y y y是由真实值组成的向量 { y 1 , y 2 . . . , y n } \{y_1,y_2...,y_n\} {y1,y2...,yn}。
2.2 推导正规方程
1. 转置公式
-
( m A ) T = m A T (mA)^T = mA^T (mA)T=mAT，其中m是常数
-
( A + B ) T = A T + B T (A + B)^T = A^T + B^T (A+B)T=AT+BT
-
( A B ) T = B T A T (AB)^T = B^TA^T (AB)T=BTAT
-
( A T ) T = A (A^T)^T = A (AT)T=A
2. 求导公式如下
- ∂ X T ∂ X = I \frac{\partial X^T}{\partial X} = I ∂X∂XT=I，求解出来是单位矩阵
- ∂ X T A ∂ X = A \frac{\partial X^TA}{\partial X} = A ∂X∂XTA=A
- ∂ A X T ∂ X = A \frac{\partial AX^T}{\partial X} = A ∂X∂AXT=A
- ∂ A X ∂ X = A T \frac{\partial AX}{\partial X} = A^T ∂X∂AX=AT
- ∂ X A ∂ X = A T \frac{\partial XA}{\partial X} = A^T ∂X∂XA=AT
- ∂ X T A X ∂ X = ( A + A T ) X \frac{\partial X^TAX}{\partial X} = (A + A^T)X ∂X∂XTAX=(A+AT)X，A不是对称矩阵
- ∂ X T A X ∂ X = 2 A X \frac{\partial X^TAX}{\partial X} = 2AX ∂X∂XTAX=2AX，A是对称矩阵
3. 推导正规方程 θ \theta θ 的解
-
矩阵乘法公式展开：
-
J ( θ ) = 1 2 ( X θ − y ) T ( X θ − y ) J(\theta) = \frac{1}{2}(X\theta - y)^T(X\theta - y) J(θ)=21(Xθ−y)T(Xθ−y)
-
J ( θ ) = 1 2 ( θ T X T − y T ) ( X θ − y ) J(\theta) = \frac{1}{2}(\theta^TX^T - y^T)(X\theta - y) J(θ)=21(θTXT−yT)(Xθ−y)
-
J ( θ ) = 1 2 ( θ T X T X θ − θ T X T y − y T X θ + y T y ) J(\theta) = \frac{1}{2}(\theta^TX^TX\theta - \theta^TX^Ty -y^TX\theta + y^Ty) J(θ)=21(θTXTXθ−θTXTy−yTXθ+yTy)
-
-
进行求导（注意X、y是已知量， θ \theta θ 是未知数）：
- J ′ ( θ ) = 1 2 ( θ T X T X θ − θ T X T y − y T X θ + y T y ) ′ J'(\theta) = \frac{1}{2}(\theta^TX^TX\theta - \theta^TX^Ty -y^TX\theta + y^Ty)' J′(θ)=21(θTXTXθ−θTXTy−yTXθ+yTy)′
-
根据上面求导公式进行运算：
- J ′ ( θ ) = 1 2 ( X T X θ + ( θ T X T X ) T − X T y − ( y T X ) T ) J'(\theta) = \frac{1}{2}(X^TX\theta + (\theta^TX^TX)^T-X^Ty - (y^TX)^T) J′(θ)=21(XTXθ+(θTXTX)T−XTy−(yTX)T)
- J ′ ( θ ) = 1 2 ( X T X θ + X T X θ − X T y − X T y ) J'(\theta) = \frac{1}{2}(X^TX\theta + X^TX\theta -X^Ty - X^Ty) J′(θ)=21(XTXθ+XTXθ−XTy−XTy)
- J ′ ( θ ) = 1 2 ( 2 X T X θ − 2 X T y ) J'(\theta) = \frac{1}{2}(2X^TX\theta -2X^Ty) J′(θ)=21(2XTXθ−2XTy)
- J ′ ( θ ) = X T X θ − X T y J'(\theta) =X^TX\theta -X^Ty J′(θ)=XTXθ−XTy
- J ′ ( θ ) = X T ( X θ − y ) J'(\theta) =X^T(X\theta -y) J′(θ)=XT(Xθ−y)，矩阵运算分配律
-
令导数 J ′ ( θ ) = 0 J'(\theta) = 0 J′(θ)=0：
- 0 = X T X θ − X T y 0 =X^TX\theta -X^Ty 0=XTXθ−XTy
- X T X θ = X T y X^TX\theta = X^Ty XTXθ=XTy
-
矩阵没有除法，使用逆矩阵进行转化：
- ( X T X ) − 1 X T X θ = ( X T X ) − 1 X T y (X^TX)^{-1}X^TX\theta = (X^TX)^{-1}X^Ty (XTX)−1XTXθ=(XTX)−1XTy
- I θ = ( X T X ) − 1 X T y I\theta = (X^TX)^{-1}X^Ty Iθ=(XTX)−1XTy
- θ = ( X T X ) − 1 X T y \theta = (X^TX)^{-1}X^Ty θ=(XTX)−1XTy
4. 最重要的结论
- 实在推不出来也没有关系，只需要记住结论：
θ = ( X T X ) − 1 X T y \theta = (X^TX)^{-1}X^Ty θ=(XTX)−1XTy - 其中的 θ \theta θ 即是方程的解， θ \theta θ 取该值时，损失函数的值全局最小。
2.3 正规方程练习
1. 二元一次方程
{ x + y = 14 2 x − y = 10 \begin{cases} x + y=14\\ 2x - y = 10\\ \end{cases} {x+y=142x−y=10
- 上式也相当于，有两组数据，求函数的系数
w
1
w_1
w1，
w
2
w_2
w2：
{ w 1 + w 2 = 14 2 w 1 − w 2 = 10 \begin{cases} w_1 + w_2=14\\ 2w_1 - w_2 = 10\\ \end{cases} {w1+w2=142w1−w2=10
X = np.array([[1, 1], [2, -1]])
y = np.array([14, 10])
# linalg 线性代数，slove计算线性回归问题，内置的利用正规方程求解最优解的函数
theta = np.linalg.solve(X, y)
print(theta)
# 正规方程
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(theta)
"""
输出：
[8. 6.]
[8. 6.]
"""
- 得到函数： f ( x , y ) = 8 x + 6 y f(x,y)=8x+6y f(x,y)=8x+6y
2. 八元一次方程
{ 14 x 2 + 8 x 3 + 5 x 5 + − 2 x 6 + 9 x 7 + − 3 x 8 = 339 − 4 x 1 + 10 x 2 + 6 x 3 + 4 x 4 + − 14 x 5 + − 2 x 6 + − 14 x 7 + 8 x 8 = − 114 − 1 x 1 + − 6 x 2 + 5 x 3 + − 12 x 4 + 3 x 5 + − 3 x 6 + 2 x 7 + − 2 x 8 = 30 5 x 1 + − 2 x 2 + 3 x 3 + 10 x 4 + 5 x 5 + 11 x 6 + 4 x 7 + − 8 x 8 = 126 − 15 x 1 + − 15 x 2 + − 8 x 3 + − 15 x 4 + 7 x 5 + − 4 x 6 + − 12 x 7 + 2 x 8 = − 395 11 x 1 + − 10 x 2 + − 2 x 3 + 4 x 4 + 3 x 5 + − 9 x 6 + − 6 x 7 + 7 x 8 = − 87 − 14 x 1 + 4 x 3 + − 3 x 4 + 5 x 5 + 10 x 6 + 13 x 7 + 7 x 8 = 422 − 3 x 1 + − 7 x 2 + − 2 x 3 + − 8 x 4 + − 6 x 6 + − 5 x 7 + − 9 x 8 = − 309 \begin{cases}14x_2 + 8x_3 + 5x_5 + -2x_6 + 9x_7 + -3x_8 = 339\\-4x_1 + 10x_2 + 6x_3 + 4x_4 + -14x_5 + -2x_6 + -14x_7 + 8x_8 = -114\\-1x_1 + -6x_2 + 5x_3 + -12x_4 + 3x_5 + -3x_6 + 2x_7 + -2x_8 = 30\\5x_1 + -2x_2 + 3x_3 + 10x_4 + 5x_5 + 11x_6 + 4x_7 + -8x_8 = 126\\-15x_1 + -15x_2 + -8x_3 + -15x_4 + 7x_5 + -4x_6 + -12x_7 + 2x_8 = -395\\11x_1 + -10x_2 + -2x_3 + 4x_4 + 3x_5 + -9x_6 + -6x_7 + 7x_8 = -87\\-14x_1 + 4x_3 + -3x_4 + 5x_5 + 10x_6 + 13x_7 + 7x_8 = 422\\-3x_1 + -7x_2 + -2x_3 + -8x_4 + -6x_6 + -5x_7 + -9x_8 = -309\end{cases} ⎩ ⎨ ⎧14x2+8x3+5x5+−2x6+9x7+−3x8=339−4x1+10x2+6x3+4x4+−14x5+−2x6+−14x7+8x8=−114−1x1+−6x2+5x3+−12x4+3x5+−3x6+2x7+−2x8=305x1+−2x2+3x3+10x4+5x5+11x6+4x7+−8x8=126−15x1+−15x2+−8x3+−15x4+7x5+−4x6+−12x7+2x8=−39511x1+−10x2+−2x3+4x4+3x5+−9x6+−6x7+7x8=−87−14x1+4x3+−3x4+5x5+10x6+13x7+7x8=422−3x1+−7x2+−2x3+−8x4+−6x6+−5x7+−9x8=−309
# 上面八元一次方程对应的X数据
X = np.array([[0, 14, 8, 0, 5, -2, 9, -3],
[-4, 10, 6, 4, -14, -2, -14, 8],
[-1, -6, 5, -12, 3, -3, 2, -2],
[5, -2, 3, 10, 5, 11, 4, -8],
[-15, -15, -8, -15, 7, -4, -12, 2],
[11, -10, -2, 4, 3, -9, -6, 7],
[-14, 0, 4, -3, 5, 10, 13, 7],
[-3, -7, -2, -8, 0, -6, -5, -9]])
# 对应的y
y = np.array([339, -114, 30, 126, -395, -87, 422, -309])
theta = np.linalg.solve(X, y)
print(theta)
"""
输出：
[ 1. 5. 15. 3. 8. 4. 17. 12.]
"""
- 得到函数： f ( X ) = x 1 + 5 x 2 + 15 x 3 + 3 x 4 + 8 x 5 + 4 x 6 + 17 x 7 + 12 x 8 f(X)=x_1+5x_2+15x_3+3x_4+8x_5+4x_6+17x_7+12x_8 f(X)=x1+5x2+15x3+3x4+8x5+4x6+17x7+12x8
2.4 使用sklearn计算多元线性方程
1. 计算不带截距的线性方程
- 这里使用到的数据还是2.3中的八元一次方程：
from sklearn.linear_model import LinearRegression
# 上面八元一次方程对应的X数据
X = np.array([[0, 14, 8, 0, 5, -2, 9, -3],
[-4, 10, 6, 4, -14, -2, -14, 8],
[-1, -6, 5, -12, 3, -3, 2, -2],
[5, -2, 3, 10, 5, 11, 4, -8],
[-15, -15, -8, -15, 7, -4, -12, 2],
[11, -10, -2, 4, 3, -9, -6, 7],
[-14, 0, 4, -3, 5, 10, 13, 7],
[-3, -7, -2, -8, 0, -6, -5, -9]])
# 对应的y
y = np.array([339, -114, 30, 126, -395, -87, 422, -309])
# 实例化模型
model = LinearRegression(fit_intercept=False) # 设置fit_intercept=False表示不计算截距
# X: 数据, y: 目标值
model.fit(X, y) # 训练模型
theta = model.coef_ # 拿到方程的解，参数theta，不含斜率
print(theta)
b = model.intercept_ # 拿到截距
print(b)
"""
输出：
[ 1. 5. 15. 3. 8. 4. 17. 12.]
0.0
"""
2. 计算带截距的线性方程
- 改造数据，给八元一次方程增加截距，设截距为12：
{ 14 x 2 + 8 x 3 + 5 x 5 + − 2 x 6 + 9 x 7 + − 3 x 8 + 12 = 339 + 12 − 4 x 1 + 10 x 2 + 6 x 3 + 4 x 4 + − 14 x 5 + − 2 x 6 + − 14 x 7 + 8 x 8 + 12 = − 114 + 12 − 1 x 1 + − 6 x 2 + 5 x 3 + − 12 x 4 + 3 x 5 + − 3 x 6 + 2 x 7 + − 2 x 8 + 12 = 30 + 12 5 x 1 + − 2 x 2 + 3 x 3 + 10 x 4 + 5 x 5 + 11 x 6 + 4 x 7 + − 8 x 8 + 12 = 126 + 12 − 15 x 1 + − 15 x 2 + − 8 x 3 + − 15 x 4 + 7 x 5 + − 4 x 6 + − 12 x 7 + 2 x 8 + 12 = − 395 + 12 11 x 1 + − 10 x 2 + − 2 x 3 + 4 x 4 + 3 x 5 + − 9 x 6 + − 6 x 7 + 7 x 8 + 12 = − 87 + 12 − 14 x 1 + 4 x 3 + − 3 x 4 + 5 x 5 + 10 x 6 + 13 x 7 + 7 x 8 + 12 = 422 + 12 − 3 x 1 + − 7 x 2 + − 2 x 3 + − 8 x 4 + − 6 x 6 + − 5 x 7 + − 9 x 8 + 12 = − 309 + 12 \begin{cases}14x_2 + 8x_3 + 5x_5 + -2x_6 + 9x_7 + -3x_8+12= 339+12\\-4x_1 + 10x_2 + 6x_3 + 4x_4 + -14x_5 + -2x_6 + -14x_7 + 8x_8 +12= -114+12\\-1x_1 + -6x_2 + 5x_3 + -12x_4 + 3x_5 + -3x_6 + 2x_7 + -2x_8 +12= 30+12\\5x_1 + -2x_2 + 3x_3 + 10x_4 + 5x_5 + 11x_6 + 4x_7 + -8x_8+12 = 126+12\\-15x_1 + -15x_2 + -8x_3 + -15x_4 + 7x_5 + -4x_6 + -12x_7 + 2x_8+12= -395+12\\11x_1 + -10x_2 + -2x_3 + 4x_4 + 3x_5 + -9x_6 + -6x_7 + 7x_8 +12= -87+12\\-14x_1 + 4x_3 + -3x_4 + 5x_5 + 10x_6 + 13x_7 + 7x_8 +12= 422+12\\-3x_1 + -7x_2 + -2x_3 + -8x_4 + -6x_6 + -5x_7 + -9x_8 +12= -309+12\end{cases} ⎩ ⎨ ⎧14x2+8x3+5x5+−2x6+9x7+−3x8+12=339+12−4x1+10x2+6x3+4x4+−14x5+−2x6+−14x7+8x8+12=−114+12−1x1+−6x2+5x3+−12x4+3x5+−3x6+2x7+−2x8+12=30+125x1+−2x2+3x3+10x4+5x5+11x6+4x7+−8x8+12=126+12−15x1+−15x2+−8x3+−15x4+7x5+−4x6+−12x7+2x8+12=−395+1211x1+−10x2+−2x3+4x4+3x5+−9x6+−6x7+7x8+12=−87+12−14x1+4x3+−3x4+5x5+10x6+13x7+7x8+12=422+12−3x1+−7x2+−2x3+−8x4+−6x6+−5x7+−9x8+12=−309+12
# 上面八元一次方程对应的X数据
X = np.array([[0, 14, 8, 0, 5, -2, 9, -3],
[-4, 10, 6, 4, -14, -2, -14, 8],
[-1, -6, 5, -12, 3, -3, 2, -2],
[5, -2, 3, 10, 5, 11, 4, -8],
[-15, -15, -8, -15, 7, -4, -12, 2],
[11, -10, -2, 4, 3, -9, -6, 7],
[-14, 0, 4, -3, 5, 10, 13, 7],
[-3, -7, -2, -8, 0, -6, -5, -9]])
# 对应的y
y = np.array([339, -114, 30, 126, -395, -87, 422, -309])
# 设截距为12
y += 12 # 目标值向上移动，相当于b = 12
- 将截距 b b b 吸入矩阵 X X X：
# 在X矩阵的最后插入一列1
X = np.concatenate([X, np.full(shape=(8,1), fill_value=1)], axis=1)
print(X)
"""
输出：
[[ 0 14 8 0 5 -2 9 -3 1]
[ -4 10 6 4 -14 -2 -14 8 1]
[ -1 -6 5 -12 3 -3 2 -2 1]
[ 5 -2 3 10 5 11 4 -8 1]
[-15 -15 -8 -15 7 -4 -12 2 1]
[ 11 -10 -2 4 3 -9 -6 7 1]
[-14 0 4 -3 5 10 13 7 1]
[ -3 -7 -2 -8 0 -6 -5 -9 1]]
"""
- 使用正规方程和
sklearn分别对其进行计算：
# 正规方程
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print("正规方程：\n", theta)
# sklearn
model = LinearRegression(fit_intercept=False) # 还是不要计算截距，因为我们已经将截距吸入矩阵了，当成未知量系数来处理了
model.fit(X, y)
print("sklearn：\n", model.coef_, model.intercept_)
"""
输出：
正规方程：
[ 42.4453125 -6.3515625 83.359375 -1.9921875 10.9296875
-5.953125 -0.39648438 -0.78515625 -293.140625 ]
sklearn：
[-0.28902686 3.99875919 17.18947808 4.3391324 8.59756611 2.82031444
16.74961517 11.58693282 4.16860692] 0.0
"""
- 和我们期望的标准答案
f
(
X
)
=
x
1
+
5
x
2
+
15
x
3
+
3
x
4
+
8
x
5
+
4
x
6
+
17
x
7
+
12
x
8
+
12
f(X)=x_1+5x_2+15x_3+3x_4+8x_5+4x_6+17x_7+12x_8+12
f(X)=x1+5x2+15x3+3x4+8x5+4x6+17x7+12x8+12 相差甚远，这是为什么？
- 大家不妨思考线性代数中的一个常识，8个未知量，最少需要多少个等式才能有固定解？答案是8个。
- 当我们将 b b b 吸入矩阵 X X X 后，本质上是让8元一次方程，变成了9元一次方程，多了一个未知量。而我们只提供了8个等式，所以有无穷多个解。
- 为了解决这个问题，我们需要增加一个等式。
- 增加一个等式：
# 增加一个等式
w = np.array([1, 5, 15, 3, 8, 4, 17, 12]) # 我们期望的标准答案的系数
X9 = np.random.randint(-15, 15, size=8) # 随机生成一些自变量x
y9 = X9.dot(w) + 12 # 计算根据X9+标准答案计算出的y值
y = np.concatenate([y, [y9]]) # 将y9加入向量y
X9 = np.concatenate([X9, [1]]) # 给X9也加入一列1
X = np.concatenate([X, [X9]]) # 将X9放入X的最后一行
# 正规方程
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print("正规方程：\n", theta)
# sklearn
model = LinearRegression(fit_intercept=False)
model.fit(X, y)
print("sklearn：\n", model.coef_, model.intercept_)
"""
输出：
正规方程：
[ 1. 5. 15. 3. 8. 4. 17. 12. 12.]
sklearn：
[ 1. 5. 15. 3. 8. 4. 17. 12. 12.] 0.0
"""
sklearn提供了不需要吸入 b b b，带截距计算的方法（但是注意，它只是对吸入的操作进行了封装，本质上还是将 b b b 当做系数来处理了，有9个未知量）：
model = LinearRegression(fit_intercept=True)
model.fit(X[:, :-1], y) # 在上面数据的基础上，去掉了X最后的一列1
print("sklearn：\n", model.coef_, model.intercept_)
"""
输出：
sklearn：
[ 1. 5. 15. 3. 8. 4. 17. 12.] 12.000000000000028
"""
2.5 凸函数
1. 何为凸函数（主观）
- 主观上判断，只有一个凸起或凹陷，就是凸函数：
- 注意，凹凸是中文上的概念，而在英文中，只有“凸”的概念。英文中将左图的“凹”函数，称为“上凸”函数，右图则是“下凸”函数。
- 非凸函数就是不只有一个凹陷或凸起：
- 来一个更加立体的效果图（左边是凸函数，右边是非凸函数）：
2. 判定凸函数的好处
- 判定损失函数是凸函数的好处在于我们可能很肯定的知道我们求得的极值即最优解，一定是全局最优解。如果是非凸函数，那就不一定可以获取全局最优解。
3. 用客观的数学方式证明凸函数
-
判定凸函数的方式非常多，其中一个方法是看黑塞矩阵是否是半正定的。
-
黑塞矩阵(hessian matrix)是由目标函数在点 X 处的二阶偏导数组成的对称矩阵。
-
对于我们的式子来说就是在导函数的基础上再次对 θ \theta θ 来求偏导，结果就是 X T X X^TX XTX。所谓正定就是 X T X X^TX XTX 的特征值全为正数，半正定就是 X T X X^TX XTX 的特征值大于等于 0， 就是半正定。
-
这里我们对 J ( θ ) J(\theta) J(θ) 损失函数求二阶导数的黑塞矩阵是 X T X X^TX XTX ，得到的一定是半正定的，自己和自己做点乘嘛！
-
这里不用数学推导证明这一点。在机器学习中往往损失函数都是凸函数，到深度学习中损失函数往往是非凸函数，即找到的解未必是全局最优，只要模型堪用就好！
-
机器学习特点是：不强调模型 100% 正确，只要是有价值的，堪用的，就Okay！
3. 线性回归算法推导
3.1 深入理解回归
回归简单来说就是“回归平均值”(regression to the mean)。但是这里的 mean 并不是把历史数据直接当成未来的预测值，而是会把期望值当作预测值。追根溯源回归这个词是一个叫高尔顿的人发明的，他通过大量观察数据发现：
- 父亲比较高，儿子也比较高；父亲比较矮，那么儿子也比较矮！
- 父亲是 1.98，儿子肯定很高，但有可能不会达到1.98
- 父亲是 1.69，儿子肯定不高，但是有可能比 1.69 高
大自然让我们回归到一定的区间之内，这就是大自然神奇的力量。
人类社会很多事情都被大自然这种神奇的力量只配置：身高、体重、智商、相貌……
这种神秘的力量就叫正态分布。大数学家高斯，深入研究了正态分布，最终推导出了线性回归的原理：最小二乘法！所以正态分布又叫高斯分布。
3.2 误差分析
误差 ε i \varepsilon_i εi 等于第 i i i 个样本实际的值 y i y_i yi 减去预测的值 y ^ \hat{y} y^ ，公式可以表达为如下：
ε
i
=
∣
y
(
i
)
−
y
^
∣
\varepsilon_i = |y^{(i)} - \hat{y}|
εi=∣y(i)−y^∣
ε
i
=
∣
y
(
i
)
−
W
T
x
(
i
)
∣
\varepsilon_i = |y^{(i)} - W^Tx^{(i)}|
εi=∣y(i)−WTx(i)∣
假定所有的样本的误差都是独立的，有上下的震荡。震荡认为是随机变量，足够多的随机变量叠加之后形成的分布，它服从的就是正态分布，因为它是正常状态下的分布，也就是高斯分布！
均值是某一个值，方差是某一个值。 方差我们先不管，均值我们总有办法让它去等于零 0 ，可以通过截距项
b
b
b 来平移整体分布的位置从而使得
μ
=
0
μ=0
μ=0， 所有误差我们就可以认为是独立分布的，1<=i<=n，服从均值为 0，方差为某定值的高斯分布。
机器学习中我们假设误差符合均值为0，方差为定值的正态分布！！！
3.3 极大似然估计
1. 概念
-
极大似然估计（Maximum Likelihood Estimation，简称MLE）是一种统计方法，用于估计统计模型的参数。它基于这样的思想：给定一组数据，我们想要找到模型参数的值，使得这些参数下观察到的数据出现的概率（似然性）最大。
-
具体来说，假设我们有一个概率模型，它依赖于一些未知参数 θ \theta θ 。极大似然估计的目标是找到参数 θ \theta θ 的值，使得给定这些参数时，观测到的数据出现的概率（似然函数）最大。
2. 举例
-
假如有一个罐子，里面有黑白两种颜色的球，数目多少不知，两种颜色的比例也不知。我们想知道罐中白球和黑球的比例，但我们不能把罐中的球全部拿出来数。现在我们可以每次任意从已经摇匀的罐中拿一个球出来，记录球的颜色，然后把拿出来的球再放回罐中。这个过程可以重复，我们可以用记录的球的颜色来估计罐中黑白球的比例。假如在前面的一百次重复记录中，有七十次是白球，请问罐中白球所占的比例最有可能是多少？大家凭直觉，也能猜出来是 70%。
-
下面是详细推导过程：
-
设取出白球的概率是 p p p，取出黑球是 1 − p 1-p 1−p（罐子中非黑即白）
-
罐子中取一个请问是白球的概率是多少？
p p p -
罐子中取两个球，两个球都是白色，概率是多少？
p 2 p^2 p2 -
罐子中取5个球都是白色，概率是多少？
p 5 p^5 p5 -
罐子中取10个球，9个是白色，一个是黑色，概率是多少呢？
C 10 9 p 9 ( 1 − p ) 或 C 10 1 p 9 ( 1 − p ) C_{10}^9p^9(1-p) 或 C_{10}^1p^9(1-p) C109p9(1−p)或C101p9(1−p) -
罐子取100个球，70次是白球，30次是黑球，概率是多少？
P = C 100 30 p 70 ( 1 − p ) 30 P = C_{100}^{30}p^{70}(1-p)^{30} P=C10030p70(1−p)30 -
根据极大似然估计，我们要求的一个 p p p 值，使得 P P P 最大。
-
对上式求导：
P ′ = 70 p 69 ∗ ( 1 − p ) 30 + p 70 ∗ 30 ∗ ( 1 − p ) 29 ( − 1 ) P' = 70p^{69}*(1-p)^{30} + p^{70}*30*(1-p)^{29}(-1) P′=70p69∗(1−p)30+p70∗30∗(1−p)29(−1) -
令导数为0：
0 = 70 p 69 ∗ ( 1 − p ) 30 + p 70 ∗ 30 ∗ ( 1 − p ) 29 ( − 1 ) 0 = 70p^{69}*(1-p)^{30} +p^{70}*30*(1-p)^{29}(-1) 0=70p69∗(1−p)30+p70∗30∗(1−p)29(−1)
公式化简：
0 = 70 ∗ ( 1 − p ) − p ∗ 30 0 = 70*(1-p) - p*30 0=70∗(1−p)−p∗30
0 = 70 − 100 ∗ p 0 = 70 - 100*p 0=70−100∗p
p = 70 % p = 70\% p=70%
-
3.4 高斯分布——概率密度函数
最常见的连续概率分布是正态分布，也叫高斯分布，而这正是我们所需要的，其概率密度函数如下:
公式如下：
f ( x ∣ μ , σ 2 ) = 1 2 π σ e − ( x − μ ) 2 2 σ 2 f(x|\mu,\sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x - \mu)^2}{2\sigma^2}} f(x∣μ,σ2)=2πσ1e−2σ2(x−μ)2
随着参数 μ \mu μ 和 σ \sigma σ 变化，概率分布也产生变化。 下面重要的步骤来了，我们要把一组数据误差出现的总似然，也就是一组数据对应误差出现的整体可能性表达出来。数据的误差我们假设服从一个高斯分布，并且通过截距项来平移整体分布的位置从而使得 μ = 0 μ=0 μ=0，所以样本的误差我们可以表达其概率密度函数的值如下:
f ( ε ∣ μ = 0 , σ 2 ) = 1 2 π σ e − ( ε − 0 ) 2 2 σ 2 f(\varepsilon|\mu = 0,\sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(\varepsilon - 0)^2}{2\sigma^2}} f(ε∣μ=0,σ2)=2πσ1e−2σ2(ε−0)2
简化如下：
f ( ε ∣ 0 , σ 2 ) = 1 2 π σ e − ε 2 2 σ 2 f(\varepsilon| 0,\sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{\varepsilon^2}{2\sigma^2}} f(ε∣0,σ2)=2πσ1e−2σ2ε2
3.5 误差总似然
和前面黑球白球问题类似，也是一个累乘问题
P = ∏ i = 0 n f ( ε i ∣ 0 , σ 2 ) = ∏ i = 0 n 1 2 π σ e − ε i 2 2 σ 2 P = \prod\limits_{i = 0}^{n}f(\varepsilon_i|0,\sigma^2) = \prod\limits_{i = 0}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{\varepsilon_i^2}{2\sigma^2}} P=i=0∏nf(εi∣0,σ2)=i=0∏n2πσ1e−2σ2εi2
根据前面公式 ε i = ∣ y ( i ) − W T x ( i ) ∣ \varepsilon_i = |y^{(i)} - W^Tx^{(i)}| εi=∣y(i)−WTx(i)∣ 可以推导出来如下公式：
P = ∏ i = 0 n f ( ε i ∣ 0 , σ 2 ) = ∏ i = 0 n 1 2 π σ e − ( y ( i ) − W T x ( i ) ) 2 2 σ 2 P = \prod\limits_{i = 0}^{n}f(\varepsilon_i|0,\sigma^2) = \prod\limits_{i = 0}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)} - W^Tx^{(i)})^2}{2\sigma^2}} P=i=0∏nf(εi∣0,σ2)=i=0∏n2πσ1e−2σ2(y(i)−WTx(i))2
公式中的未知变量就是 W T W^T WT，即方程的系数，系数包含截距。如果，把上面当成一个函数，就是概率 P P P 关于 W W W 的函数！其余符号，都是常量！
P ( W ) = ∏ i = 0 n 1 2 π σ e − ( y ( i ) − W T x ( i ) ) 2 2 σ 2 P(W)= \prod\limits_{i = 0}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)}- W^Tx^{(i)})^2}{2\sigma^2}} P(W)=i=0∏n2πσ1e−2σ2(y(i)−WTx(i))2
现在问题，就变换成了，求最大似然问题。但是累乘的最大似然，求解是非常麻烦的！接下来，我们通过，求对数把累乘问题，转变为累加问题。
3.6 最小二乘法MSE
P ( W ) = ∏ i = 0 n 1 2 π σ e − ( y ( i ) − W T x ( i ) ) 2 2 σ 2 P(W) = \prod\limits_{i = 0}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)} - W^Tx^{(i)})^2}{2\sigma^2}} P(W)=i=0∏n2πσ1e−2σ2(y(i)−WTx(i))2
根据对数，单调性，对上面公式求自然底数e的对数，效果不变：
l o g e ( P ( W ) ) = l o g e ( ∏ i = 0 n 1 2 π σ e − ( y ( i ) − W T x ( i ) ) 2 2 σ 2 ) log_e(P(W)) = log_e(\prod\limits_{i = 0}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)} - W^Tx^{(i)})^2}{2\sigma^2}}) loge(P(W))=loge(i=0∏n2πσ1e−2σ2(y(i)−WTx(i))2)
接下来 l o g log log 函数继续为你带来惊喜，数学上连乘是个大麻烦，即使交给计算机去求解它也得哭出声来。惊喜是:
l
o
g
e
(
P
(
W
)
)
=
l
o
g
e
(
∏
i
=
0
n
1
2
π
σ
e
−
(
y
(
i
)
−
W
T
x
(
i
)
)
2
2
σ
2
)
log_e(P(W)) = log_e(\prod\limits_{i = 0}^{n}\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)} - W^Tx^{(i)})^2}{2\sigma^2}})
loge(P(W))=loge(i=0∏n2πσ1e−2σ2(y(i)−WTx(i))2)
=
∑
i
=
0
n
l
o
g
e
(
1
2
π
σ
e
−
(
y
(
i
)
−
W
T
x
(
i
)
)
2
2
σ
2
)
=\sum\limits_{i = 0}^{n}log_e(\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(y^{(i)} - W^Tx^{(i)})^2}{2\sigma^2}})
=i=0∑nloge(2πσ1e−2σ2(y(i)−WTx(i))2)累乘问题变成了累加问题。
乘风破浪，继续推导：
=
∑
i
=
0
n
(
l
o
g
e
1
2
π
σ
−
(
y
(
i
)
−
W
T
x
(
i
)
)
2
2
σ
2
)
=\sum\limits_{i = 0}^{n}(log_e\frac{1}{\sqrt{2\pi}\sigma} - \frac{(y^{(i)} - W^Tx^{(i)})^2}{2\sigma^2})
=i=0∑n(loge2πσ1−2σ2(y(i)−WTx(i))2)
=
∑
i
=
0
n
(
l
o
g
e
1
2
π
σ
−
1
σ
2
⋅
1
2
(
y
(
i
)
−
W
T
x
(
i
)
)
2
)
=\sum\limits_{i = 0}^{n}(log_e\frac{1}{\sqrt{2\pi}\sigma} - \frac{1}{\sigma^2}\cdot\frac{1}{2}(y^{(i)} - W^Tx^{(i)})^2)
=i=0∑n(loge2πσ1−σ21⋅21(y(i)−WTx(i))2)
上面公式是最大似然求对数后的变形，其中 π 、 σ \pi、\sigma π、σ都是常量，而 ( y i − W T x i ) 2 (y_i - WTx_i)2 (yi−WTxi)2肯定大于零！上面求最大值问题，即可转变为如下求最小值问题：
L
(
W
)
=
1
2
∑
i
=
0
n
(
y
(
i
)
−
W
T
x
(
i
)
)
2
L(W) = \frac{1}{2}\sum\limits_{i = 0}^n(y^{(i)} - W^Tx^{(i)})^2
L(W)=21i=0∑n(y(i)−WTx(i))2
L代表Loss，表示损失函数，损失函数越小，那么上面最大似然就越大。
有的书本上公式，也可以这样写，用 J ( θ ) J(\theta) J(θ)表示一个意思， θ \theta θ 的角色就是W：
J ( θ ) = 1 2 ∑ i = 1 n ( y ( i ) − θ T x ( i ) ) 2 = 1 2 ∑ i = 1 n ( θ T x ( i ) − y ( i ) ) 2 J(\theta) = \frac{1}{2}\sum\limits_{i = 1}^{n}(y^{(i)} - \theta^Tx^{(i)})^2 = \frac{1}{2}\sum\limits_{i = 1}^{n}(\theta^Tx^{(i)} - y^{(i)})^2 J(θ)=21i=1∑n(y(i)−θTx(i))2=21i=1∑n(θTx(i)−y(i))2
进一步提取：
J ( θ ) = 1 2 ∑ i = 1 n ( h θ ( x ( i ) ) − y ( i ) ) 2 J(\theta) = \frac{1}{2}\sum\limits_{i = 1}^n(h_{\theta}(x^{(i)}) - y^{(i)})^2 J(θ)=21i=1∑n(hθ(x(i))−y(i))2
写成矩阵形式就是：
J
(
θ
)
=
1
2
(
X
θ
−
y
)
T
(
X
θ
−
y
)
J(\theta) = \frac{1}{2}(X\theta - y)^T(X\theta - y)
J(θ)=21(Xθ−y)T(Xθ−y)
到这里，我们就已经推导出来了 MSE 损失函数
J
(
θ
)
J(\theta)
J(θ)，从公式我们也可以看出来 MSE 名字的来历，即mean squared error均方误差，上式也叫做最小二乘法！
3.7 归纳总结
这种最小二乘法估计，我们可以认为，假定了误差服从正太分布，认为样本误差的出现是随机的，独立的，使用最大似然估计思想，利用损失函数最小化 MSE 就能求出最优解！
所以反过来说，如果我们的数据误差不是互相独立的，或者不是随机出现的，那么就不适合去假设为正太分布，就不能去用正太分布的概率密度函数带入到总似然的函数中，故而就不能用 MSE 作为损失函数去求解最优解！所以，最小二乘法不是万能的。
还有譬如假设误差服从泊松分布，或其他分布那就得用其他分布的概率密度函数去推导出损失函数了。
所以有时我们也可以把线性回归看成是广义线性回归。比如，逻辑回归，泊松回归都属于广义线性回归的一种，这里我们线性回归可以说是最小二乘线性回归。
4. 线性回归实战
4.1 使用正规方程求解
4.1.1 简单线性回归
y = w x + b y = wx + b y=wx+b
一元一次方程，在机器学习中一元表示一个特征， b b b 表示截距， y y y 表示目标值。
import numpy as np
import matplotlib.pyplot as plt
# 转化成矩阵
X = np.linspace(0,10,num = 30).reshape(-1,1)
# 斜率和截距，随机生成
w = np.random.randint(1,5,size = 1)
b = np.random.randint(1,10,size = 1)
# 根据一元一次方程计算目标值y，并加上“噪声”，数据有上下波动
y = X * w + b + np.random.randn(30,1)
plt.scatter(X,y)
# 吸入b
X = np.concatenate([X,np.full(shape = (30,1),fill_value= 1)],axis = 1)
# 正规方程求解
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print('一元一次方程真实的斜率和截距是：',w, b)
print('通过正规方程求解的斜率和截距是：',theta)
# 根据求解的斜率和截距绘制线性回归线型图
plt.plot(X[:,0],X.dot(theta),color = 'green')
"""
输出：
一元一次方程真实的斜率和截距是： [4] [4]
通过正规方程求解的斜率和截距是： [[4.098669 ]
[3.70840987]]
"""
4.1.2 多元线性回归
y = w 1 x 1 + w 2 x 2 + b y = w_1x_1 + w_2x_2 + b y=w1x1+w2x2+b
二元一次方程， x 1 , x 2 x_1,x_2 x1,x2 相当于两个特征， b b b 是方程截距。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D # 绘制三维图像
# 转化成矩阵
x1 = np.random.randint(-150,150,size = (300,1))
x2 = np.random.randint(0,300,size = (300,1))
# 斜率和截距，随机生成
w = np.random.randint(1,5,size = 2)
b = np.random.randint(1,10,size = 1)
# 根据二元一次方程计算目标值y，并加上“噪声”，数据有上下波动
y = x1 * w[0] + x2 * w[1] + b + np.random.randn(300,1)
ax = plt.subplot(111, projection='3d')
ax.scatter(x1,x2,y) # 三维散点图
ax.view_init(elev=10, azim=-20) # 调整视角
# 吸入b
X = np.concatenate([x1,x2,np.full(shape = (300,1),fill_value=1)],axis = 1)
# 正规方程求解
theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print('二元一次方程真实的斜率和截距是：',w,b)
print('通过正规方程求解的斜率和截距是：',theta.reshape(-1, 1))
# 根据求解的斜率和截距绘制线性回归线型图
x = np.linspace(-150,150,100)
y = np.linspace(0,300,100)
z = x * theta[0] + y * theta[1] + theta[2]
ax.plot(x,y,z ,color = 'red')
"""
输出：
二元一次方程真实的斜率和截距是： [4 2] [5]
通过正规方程求解的斜率和截距是： [[4.00101176]
[1.99914996]
[5.06319223]]
"""
4.2 sklearn求解
4.2.1 简单线性回归
from sklearn.linear_model import LinearRegression
import numpy as np
import matplotlib.pyplot as plt
# 转化成矩阵
X = np.linspace(0,10,num = 30).reshape(-1,1)
# 斜率和截距，随机生成
w = np.random.randint(1,5,size = 1)
b = np.random.randint(1,10,size = 1)
# 根据一元一次方程计算目标值y，并加上“噪声”，数据有上下波动
y = X * w + b + np.random.randn(30,1)
plt.scatter(X,y)
# 使用scikit-learn中的线性回归求解
model = LinearRegression() # 默认fit_intercept=True
model.fit(X,y)
w_ = model.coef_
b_ = model.intercept_
print('一元一次方程真实的斜率和截距是：',w, b)
print('通过scikit-learn求解的斜率和截距是：',w_, b_)
plt.plot(X, X.dot(w_) + b_, color = 'green')
"""
输出：
一元一次方程真实的斜率和截距是： [1] [4]
通过scikit-learn求解的斜率和截距是： [[0.98841254]] [4.28075802]
"""
4.2.2 多元线性回归
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
# 转化成矩阵
x1 = np.random.randint(-150,150,size = (300,1))
x2 = np.random.randint(0,300,size = (300,1))
# 斜率和截距，随机生成
w = np.random.randint(1,5,size = 2)
b = np.random.randint(1,10,size = 1)
# 根据二元一次方程计算目标值y，并加上“噪声”，数据有上下波动
y = x1 * w[0] + x2 * w[1] + b + np.random.randn(300,1)
ax = plt.subplot(111, projection='3d')
ax.scatter(x1,x2,y) # 三维散点图
ax.view_init(elev=10, azim=-20) # 调整视角
# 构造X矩阵
X = np.concatenate([x1,x2],axis = 1)
# 使用scikit-learn中的线性回归求解
model = LinearRegression()
model.fit(X,y)
w_ = model.coef_.reshape(-1)
b_ = model.intercept_
print('二元一次方程真实的斜率和截距是：',w,b)
print('通过scikit-learn求解的斜率和截距是：',w_,b_)
# # 根据求解的斜率和截距绘制线性回归线型图
x = np.linspace(-150,150,100)
y = np.linspace(0,300,100)
z = x * w_[0] + y * w_[1] + b_
ax.plot(x,y,z ,color = 'green')
"""
输出：
二元一次方程真实的斜率和截距是： [4 4] [4]
通过scikit-learn求解的斜率和截距是： [4.00016457 3.99992066] [4.04071407]
"""
