多元线性回归分析
什么是线性回归？
线性回归，如上图所示（这里用二维的例子比较好理解），我们知道许多的 ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x n , y n ) (x_1, y_1), (x_2, y_2), ..., (x_n, y_n) ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( x n , y n ) ，即图中红色的点，通过某种方法，得到图中蓝色的线（ y = w × x + b y = w\times x + b y = w × x + b ），即求 w ， b w， b w ， b 的值；然后可以使得未知数据 x n e w x_{new} x n e w 输入方程，使得结果 y n e w y_{new} y n e w 尽可能的接近真实值。
具体的关于线性回归的定义，这里是引用西瓜书上定义
定义：
给定由 d d d 个属性描述的示例 x = ( x 1 ; x 1 ; . . . ; x d ) x = (x_1; x_1; ... ; x_d) x = ( x 1 ; x 1 ; . . . ; x d ) ，其中 x i x_i x i 是 x \textbf{x} x 在第 i i i 个属性上的取值，线性模型（linear model）试图学得一个通过属性的组合来进行预测的函数，即
f ( x ) = w 1 x 1 + w 2 x 2 + . . . + w d x d + b f(x) = w_1x_1 + w_2x_2 + ... + w_dx_d + b f ( x ) = w 1 x 1 + w 2 x 2 + . . . + w d x d + b
一般用向量形式写成
f ( x ) = w T x + b f(x) = w^Tx + b f ( x ) = w T x + b
其中 w = ( w 1 ; w 2 ; . . . ; w 1 ; ) w = (w_1;w_2;...;w_1;) w = ( w 1 ; w 2 ; . . . ; w 1 ; ) 。 w w w 和 b b b 学习到以后，模型就得以确定
通常回归分为线性回归分析和非线性回归分析，判读二者最好的方法是看每一个变量的指数，如都是1次的就是线性，否则是非线性。
基于矩阵求解（解析解法或者公式法）
推导过程
假如我们的多元线性方程有n项，目标方程为
f ( x i ) = w 1 x 1 + w 2 x 2 + ⋯ + w n x n + b f\left( {{x_i}} \right) = {w_1}{x_1} + {w_2}{x_2} + \cdots + {w_n}{x_n} + b
f
(
x
i
)
=
w
1
x
1
+
w
2
x
2
+
⋯
+
w
n
x
n
+
b
有对应的关于数据D（
( x 1 , x 2 , . . . , x n , y ) (x_1,x_2,...,x_n,y)
(
x
1
,
x
2
,
.
.
.
,
x
n
,
y
)
）,需要根据输入的数据来求解（
w i , b w_i,b
w
i
,
b
）。这里对于最后的偏置项
b b
b
，可以认为是
w 0 x 0 w_0x_0
w
0
x
0
，其中
x 0 x_0
x
0
是一个恒为1的数，求解
w 0 w_0
w
0
相当于求解
b b
b
，即
f ( x i ) = w 0 x 0 + w 1 x 1 + w 2 x 2 + ⋯ + w n x n f\left( {{x_i}} \right) = {w_0}{x_0} + {w_1}{x_1} + {w_2}{x_2} + \cdots + {w_n}{x_n}
f
(
x
i
)
=
w
0
x
0
+
w
1
x
1
+
w
2
x
2
+
⋯
+
w
n
x
n
这样写成矩阵的形式为
f ( x ) = w T x f\left( x \right) = {w^T}{x}
f
(
x
)
=
w
T
x
其中
x x
x
是一个列向量，
w w
w
是一个列向量
在数据D中有多组数据，就可以写出矩阵的形式
w T [ x 1 , x 2 , … , x n ] = [ y 1 , y 2 , … , y n ] {w^T}{ \left[ \begin{array}{c} {{x^1}} ,{{x^2}},\dots ,{{x^n}} \end{array} \right ]} ={ \left[ \begin{array}{ccc} {{y^1}},{{y^2}}, \dots ,{{y^n}} \end{array} \right ]} w T [ x 1 , x 2 , … , x n ] = [ y 1 , y 2 , … , y n ]
如何确定
w w
w
，就是在于如何可以使得
y i y^i
y
i
(目标值)和
y ^ i {{\hat y}^i}
y
^
i
(预测值)直接最小，这里使用均方差来使得误差最小化，即定义下面的目标方程
w ∗ = arg min w ∑ i = 1 m ( y i − y ^ i ) 2 = arg min w ∑ i = 1 m ( y i − w T x i ) 2 \begin{array}{l} {w^*} &= \mathop {\arg \min }\limits_w \sum\limits_{i = 1}^m {({y^i} - } {{\hat y}^i}{)^2} \\ &= \mathop {\arg \min }\limits_w \sum\limits_{i = 1}^m {({y^i} - } {w^T}{x^i}{)^2} \end{array}
w
∗
=
w
ar
g
min
i
=
1
∑
m
(
y
i
−
y
^
i
)
2
=
w
ar
g
min
i
=
1
∑
m
(
y
i
−
w
T
x
i
)
2
其中
w ∗ {w^*}
w
∗
就是我们要找的最优的参数，
m m
m
表示有m组数据。根据上式可以知道，目标函数是一个关于
w w
w
的凸函数，对于凸函数而言，一阶导数等于零对应的
w w
w
就是目标函数的最优值；为了求导方便，我们再把式子变换一下，式子中
x i x^i
x
i
是一个列向量，如果有多组数据的话就可以把他写成一个矩阵
X = ( x 1 , x 2 , . . . , x m ) X=(x^1,x^2,...,x^m)
X
=
(
x
1
,
x
2
,
.
.
.
,
x
m
)
，相应的真实值
y = ( y 1 , y 2 , . . . , y n ) T y=({y}^1,{y}^2,...,{y}^n)^T
y
=
(
y
1
,
y
2
,
.
.
.
,
y
n
)
T
是一个列向量，然后目标函数就可以转换为
E = ( y − X T w ) T ( y − X T w ) E = {(y - {X^T}w)^T}(y - {X^T}w)
E
=
(
y
−
X
T
w
)
T
(
y
−
X
T
w
)
对
w w
w
求导
∂ E ∂ w = 2 X ( X T w − y ) \frac{{\partial E}}{{\partial w}} = 2X({X^T}w - y)
∂
w
∂
E
=
2
X
(
X
T
w
−
y
)
令
∂ E ∂ w = 0 \frac{{\partial E}}{{\partial w}} =0
∂
w
∂
E
=
0
，就可以求的
w = ( X X T ) − 1 X T y w = {(X{X^T})^{ - 1}}{X^T}y
w
=
(
X
X
T
)
−
1
X
T
y
代码实现
import
numpy
as
np
def
getParams
(
x
,
y
)
:
n
=
len
(
y
)
x_0
=
np
.
ones
(
n
)
x
=
np
.
matrix
(
np
.
vstack
(
[
x_0
,
x
]
)
.
T
)
y
=
np
.
matrix
(
y
)
params
=
(
(
x
.
T
*
x
)
.
I
*
x
.
T
*
(
y
.
T
)
)
b
=
params
[
0
]
w
=
params
[
1
:
]
return
w
,
b
if
__name__
==
"__main__"
:
n_samples
=
10000
x_1
,
x_2
=
np
.
random
.
random
(
n_samples
)
,
np
.
random
.
random
(
n_samples
)
theta_1
,
theta_2
,
b
=
map
(
float
,
input
(
)
.
split
(
)
)
x
,
y
=
np
.
vstack
(
[
x_1
,
x_2
]
)
,
theta_1
*
x_1
+
theta_2
*
x_2
+
b
params
=
getParams
(
x
,
y
)
梯度下降逼近法
线性回归梯度下降法，这里的梯度下降是指计算损失函数(loss function)的梯度；这个损失函数通常是均方差函数
E ( w ) = 1 2 m ∑ i = 1 m ( y ^ i − y i ) 2 = 1 2 m ∑ i = 1 m ( w T x i + b − y i ) 2 \begin{array}{l} E(w) &= \frac{1}{{2m}}\sum\limits_{i = 1}^m {({{\hat y}^i} - {y^i}} {)^2}\\ &= \frac{1}{{2m}}\sum\limits_{i = 1}^m {({w^T}{x^i} + b - {y^i}} {)^2} \end{array}
E
(
w
)
=
2
m
1
i
=
1
∑
m
(
y
^
i
−
y
i
)
2
=
2
m
1
i
=
1
∑
m
(
w
T
x
i
+
b
−
y
i
)
2
然后通过求解目标函数对于
w w
w
的梯度和对于
b b
b
的梯度，来不断的迭代更新他们的值，使得结果更加逼近真实值
∂ E ∂ w = 1 m ∑ i = 1 m ( w T x i + b − y i ) x i = 1 m ∑ i = 1 m ( y ^ i − y i ) x i ∂ E ∂ b = 1 m ∑ i = 1 m ( w T x i + b − y i ) x i = 1 m ∑ i = 1 m ( y ^ i − y i ) \begin{array}{l} \frac{{\partial E}}{{\partial w}} &= \frac{1}{m}\sum\limits_{i = 1}^m {({w^T}{x^i} + b - {y^i}} ){x^i}\\ &= \frac{1}{m}\sum\limits_{i = 1}^m {({{\hat y}^i} - {y^i}} ){x^i} \end{array} \\ \begin{array}{l} \frac{{\partial E}}{{\partial b}}& = \frac{1}{m}\sum\limits_{i = 1}^m {({w^T}{x^i} + b - {y^i}} ){x^i}\\ &= \frac{1}{m}\sum\limits_{i = 1}^m {({{\hat y}^i} - {y^i}} ) \end{array}
∂
w
∂
E
=
m
1
i
=
1
∑
m
(
w
T
x
i
+
b
−
y
i
)
x
i
=
m
1
i
=
1
∑
m
(
y
^
i
−
y
i
)
x
i
∂
b
∂
E
=
m
1
i
=
1
∑
m
(
w
T
x
i
+
b
−
y
i
)
x
i
=
m
1
i
=
1
∑
m
(
y
^
i
−
y
i
)
迭代更新:
w = w − α ∂ E ∂ w b = b − α ∂ E ∂ b \begin{array}{l} w{\rm{ = }}w{\rm{ - }}\alpha \frac{{\partial E}}{{\partial w}}\\ \\ b = b - \alpha \frac{{\partial E}}{{\partial b}} \end{array}
w
=
w
−
α
∂
w
∂
E
b
=
b
−
α
∂
b
∂
E
其中
α \alpha
α
表示学习率(learning rate)，用于控制每一步下降的步长(步长要选择适中，步长太小下降的速度比较慢，太大则可能造成不收敛；通常也可以选择动态步长的方法，即前期使用大步长用于快速下降，后期采用小步长使用收敛精度)，这里由于使用的是逼近的方法，所以我们可以设置通过循环的次数来结束求解过程，也可以使用设置阈值的方法结束求解。
代码实现
这里代码实现，分为使用numpy的方法和使用机器学习框架(pytorch)的方法
使用numpy
import
numpy
as
np
def
gradient_descent
(
X
,
y
,
lr
=
0.01
,
threshold
=
1e
-
3
)
:
y
=
y
[
:
,
np
.
newaxis
]
W
=
np
.
ones
(
shape
=
(
1
,
X
.
shape
[
1
]
)
)
b
=
np
.
array
(
[
[
1
]
]
)
n
=
y
.
shape
[
0
]
while
True
:
pred_y
=
np
.
dot
(
X
,
W
.
T
)
+
b
loss
=
np
.
dot
(
(
y
-
pred_y
)
.
T
,
(
y
-
pred_y
)
)
/
(
2
*
n
)
if
np
.
abs
(
loss
)
<
threshold
:
break
w_gradient
=
-
(
1.0
/
n
)
*
np
.
dot
(
(
y
-
pred_y
)
.
T
,
X
)
b_gradient
=
-
(
1.0
/
n
)
*
np
.
dot
(
(
y
-
pred_y
)
.
T
,
np
.
ones
(
shape
=
[
n
,
1
]
)
)
W
=
W
-
w_gradient
b
=
b
-
b_gradient
print
(
loss
)
return
W
[
0
]
[
0
]
,
W
[
0
]
[
1
]
,
b
[
0
]
[
0
]
if
__name__
==
"__main__"
:
n_samples
=
10000
x_0
=
np
.
ones
(
n_samples
)
x_1
,
x_2
=
np
.
random
.
random
(
n_samples
)
,
np
.
random
.
random
(
n_samples
)
theta_1
,
theta_2
,
b
=
map
(
float
,
input
(
)
.
split
(
)
)
X
,
y
=
np
.
vstack
(
[
x_1
,
x_2
]
)
.
T
,
theta_1
*
x_1
+
theta_2
*
x_2
+
b
params
=
gradient_descent
(
X
,
y
)
result
=
[
]
for
param
in
params
:
result
.
append
(
'{:.02f}'
.
format
(
param
)
)
print
(
' '
.
join
(
result
)
)
使用pytorch
import
torch
import
torch
.
nn
as
nn
import
numpy
as
np
# 线性回归
class
LinerRegression
(
nn
.
Module
)
:
def
__init__
(
self
)
:
super
(
LinerRegression
,
self
)
.
__init__
(
)
self
.
pred
=
torch
.
nn
.
Linear
(
2
,
1
)
def
forward
(
self
,
x
)
:
x
=
self
.
pred
(
x
)
return
x
# 产生10000个数据
n_samples
=
10000
x_1
,
x_2
=
np
.
random
.
random
(
n_samples
)
,
np
.
random
.
random
(
n_samples
)
x_train
=
np
.
vstack
(
[
x_1
,
x_2
]
)
.
T
# 设置想要的权重，产生y
theta_1
,
theta_2
,
b
=
8.69
,
-
20.66
,
76.12
y_train
=
theta_1
*
x_1
+
theta_2
*
x_2
+
b
# 把numpy array转为为tensor
x_train
=
torch
.
from_numpy
(
x_train
)
y_train
=
torch
.
from_numpy
(
y_train
)
# 为数据添加一个维度，和预测结果相匹配
y_train
=
y_train
.
unsqueeze
(
1
)
# 调用网络
net
=
LinerRegression
(
)
.
double
(
)
# 采用随机梯度下降法优化函数
optimizer
=
torch
.
optim
.
SGD
(
net
.
parameters
(
)
,
lr
=
0.01
)
# 使用均方差作为 loss function
criterion
=
nn
.
MSELoss
(
)
# 迭代循环，当loss小于一定的值时，结束
num
=
0
while
True
:
pred_y
=
net
(
x_train
)
loss
=
criterion
(
pred_y
,
y_train
)
optimizer
.
zero_grad
(
)
loss
.
backward
(
)
optimizer
.
step
(
)
if
(
num
+
1
)
%
20
==
0
:
print
(
'Loss: {:.6f}'
.
format
(
loss
.
data
)
)
if
loss
<
1e
-
3
:
break
num
+=
1
# 保存参数
torch
.
save
(
net
.
state_dict
(
)
,
'params.pth'
)
# 加载参数
net
.
load_state_dict
(
torch
.
load
(
'params.pth'
)
)
# 输出每一个参数的值
for
name
,
param
in
net
.
named_parameters
(
)
:
if
name
==
'pred.weight'
:
print
(
'W1:{:.02f}, W2:{:.02f}'
.
format
(
param
.
data
[
0
]
[
0
]
.
numpy
(
)
,
param
.
data
[
0
]
[
1
]
.
numpy
(
)
)
)
else
:
print
(
'b:{:.02f}'
.
format
(
param
.
data
[
0
]
.
numpy
(
)
)
)
注： 文中有写错的地方，欢迎大家不吝指正！！!
作者：黑暗主宰
邮箱：shengzhanhe@gmail.com
