线性回归的概念
1、线性回归的原理
2、线性回归损失函数、代价函数、目标函数
3、优化方法(梯度下降法、牛顿法、拟牛顿法等)
4、线性回归的评估指标
5、sklearn参数详解
1、线性回归的原理
进入一家房产网,可以看到房价、面积、厅室呈现以下数据:
面积($x_1$) | 厅室数量($x_2)$ | 价格(万元)(y) |
---|---|---|
64 | 3 | 225 |
59 | 3 | 185 |
65 | 3 | 208 |
116 | 4 | 508 |
…… | …… | …… |
我们可以将价格和面积、厅室数量的关系习得为 f ( x ) = θ 0 + θ 1 x 1 + θ 2 x 2 f(x)=\theta_0+\theta_1x_1+\theta_2x_2 f(x)=θ0+θ1x1+θ2x2,使得 f ( x ) ≈ y f(x)\approx y f(x)≈y,这就是一个直观的线性回归的样式。
小练习:这是国内一个房产网站上任意搜的数据,有兴趣可以找个网站观察一下,还可以获得哪些可能影响到房价的因素?可能会如何影响到实际房价呢?线性回归的一般形式:
有数据集
{
(
x
1
,
y
1
)
,
(
x
2
,
y
2
)
,
.
.
.
,
(
x
n
,
y
n
)
}
\{(x_1,y_1),(x_2,y_2),...,(x_n,y_n)\}
{(x1,y1),(x2,y2),...,(xn,yn)},其中,
x
i
=
(
x
i
1
;
x
i
2
;
x
i
3
;
.
.
.
;
x
i
d
)
,
y
i
∈
R
x_i = (x_{i1};x_{i2};x_{i3};...;x_{id}),y_i\in R
xi=(xi1;xi2;xi3;...;xid),yi∈R
其中n表示变量的数量,d表示每个变量的维度。
可以用以下函数来描述y和x之间的关系:
如何来确定 θ \theta θ的值,使得 f ( x ) f(x) f(x)尽可能接近y的值呢?均方误差是回归中常用的性能度量,即:
J
(
θ
)
=
1
2
∑
j
=
1
n
(
h
θ
(
x
(
i
)
)
−
y
(
i
)
)
2
J(\theta)=\frac{1}{2}\sum_{j=1}^{n}(h_{\theta}(x^{(i)})-y^{(i)})^2
J(θ)=21j=1∑n(hθ(x(i))−y(i))2
我们可以选择 θ \theta θ,试图让均方误差最小化:
极大似然估计(概率角度的诠释)
下面我们用极大似然估计,来解释为什么要用均方误差作为性能度量
我们可以把目标值和变量写成如下等式:
y ( i ) = θ T x ( i ) + ϵ ( i ) y^{(i)} = \theta^T x^{(i)}+\epsilon^{(i)} y(i)=θTx(i)+ϵ(i)
ϵ \epsilon ϵ表示我们未观测到的变量的印象,即随机噪音。我们假定 ϵ \epsilon ϵ是独立同分布,服从高斯分布。(根据中心极限定理)
p ( ϵ ( i ) ) = 1 2 π σ e x p ( − ( ϵ ( i ) ) 2 2 σ 2 ) p(\epsilon^{(i)}) = \frac{1}{\sqrt{2\pi}\sigma}exp\left(-\frac{(\epsilon^{(i)})^2}{2\sigma^2}\right) p(ϵ(i))=2πσ1exp(−2σ2(ϵ(i))2)
因此,
p ( y ( i ) ∣ x ( i ) ; θ ) = 1 2 π σ e x p ( − ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) p(y^{(i)}|x^{(i)};\theta) = \frac{1}{\sqrt{2\pi}\sigma}exp\left(-\frac{(y^{(i)}-\theta^T x^{(i)})^2}{2\sigma^2}\right) p(y(i)∣x(i);θ)=2πσ1exp(−2σ2(y(i)−θTx(i))2)
我们建立极大似然函数,即描述数据遵从当前样本分布的概率分布函数。由于样本的数据集独立同分布,因此可以写成
L ( θ ) = p ( y ⃗ ∣ X ; θ ) = ∏ i = 1 n 1 2 π σ e x p ( − ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) L(\theta) = p(\vec y | X;\theta) = \prod^n_{i=1}\frac{1}{\sqrt{2\pi}\sigma}exp\left(-\frac{(y^{(i)}-\theta^T x^{(i)})^2}{2\sigma^2}\right) L(θ)=p(y∣X;θ)=i=1∏n2πσ1exp(−2σ2(y(i)−θTx(i))2)
选择 θ \theta θ,使得似然函数最大化,这就是极大似然估计的思想。
为了方便计算,我们计算时通常对对数似然函数求最大值:
显然,最大化 l ( θ ) l(\theta) l(θ)即最小化 1 2 ∑ i = 1 n ( ( y ( i ) − θ T x ( i ) ) 2 \frac{1}{2}\sum^n_{i=1}((y^{(i)}-\theta^T x^{(i)})^2 21∑i=1n((y(i)−θTx(i))2。
这一结果即均方误差,因此用这个值作为代价函数来优化模型在统计学的角度是合理的。
2、线性回归损失函数、代价函数、目标函数
- 损失函数(Loss Function):度量单样本预测的错误程度,损失函数值越小,模型就越好。
- 代价函数(Cost Function):度量全部样本集的平均误差。
- 目标函数(Object Function):代价函数和正则化函数,最终要优化的函数。
常用的损失函数包括:0-1损失函数、平方损失函数、绝对损失函数、对数损失函数等;常用的代价函数包括均方误差、均方根误差、平均绝对误差等。
当模型复杂度增加时,有可能对训练集可以模拟的很好,但是预测测试集的效果不好,出现过拟合现象,这就出现了所谓的“结构化风险”。结构风险最小化即为了防止过拟合而提出来的策略,定义模型复杂度为 J ( F ) J(F) J(F),目标函数可表示为:
m i n f ∈ F 1 n ∑ i = 1 n L ( y i , f ( x i ) ) + λ J ( F ) \underset{f\in F}{min}\, \frac{1}{n}\sum^{n}_{i=1}L(y_i,f(x_i))+\lambda J(F) f∈Fminn1i=1∑nL(yi,f(xi))+λJ(F)
有以上6个房价和面积关系的数据点,可以看到,当设定
f
(
x
)
=
∑
j
=
0
5
θ
j
x
j
f(x)=\sum_{j=0}^{5}\theta_jx_j
f(x)=∑j=05θjxj时,可以完美拟合训练集数据,但是,真实情况下房价和面积不可能是这样的关系,出现了过拟合现象。当训练集本身存在噪声时,拟合曲线对未知影响因素的拟合往往不是最好的。
通常,随着模型复杂度的增加,训练误差会减少;但测试误差会先增加后减小。我们的最终目的时试测试误差达到最小,这就是我们为什么需要选取适合的目标函数的原因。
3、线性回归的优化方法
1、梯度下降法
设定初始参数
θ
\theta
θ,不断迭代,使得
J
(
θ
)
J(\theta)
J(θ)最小化:
θ
j
:
=
θ
j
−
α
∂
J
(
θ
)
∂
θ
\theta_j:=\theta_j-\alpha\frac{\partial{J(\theta)}}{\partial\theta}
θj:=θj−α∂θ∂J(θ)
即:
θ j = θ j + α ∑ i = 1 n ( y ( i ) − f θ ( x ) ( i ) ) x j ( i ) \theta_j = \theta_j + \alpha\sum_{i=1}^{n}(y^{(i)}-f_\theta(x)^{(i)})x_j^{(i)} θj=θj+αi=1∑n(y(i)−fθ(x)(i))xj(i)
注:下标j表示第j个参数,上标i表示第i个数据点。
将所有的参数以向量形式表示,可得:
θ = θ + α ∑ i = 1 n ( y ( i ) − f θ ( x ) ( i ) ) x ( i ) \theta = \theta + \alpha\sum_{i=1}^{n}(y^{(i)}-f_\theta(x)^{(i)})x^{(i)} θ=θ+αi=1∑n(y(i)−fθ(x)(i))x(i)
由于这个方法中,参数在每一个数据点上同时进行了移动,因此称为批梯度下降法,对应的,我们可以每一次让参数只针对一个数据点进行移动,即:
θ = θ + α ( y ( i ) − f θ ( x ) ( i ) ) x ( i ) \theta = \theta + \alpha(y^{(i)}-f_\theta(x)^{(i)})x^{(i)} θ=θ+α(y(i)−fθ(x)(i))x(i)
这个算法成为随机梯度下降法,随机梯度下降法的好处是,当数据点很多时,运行效率更高;缺点是,因为每次只针对一个样本更新参数,未必找到最快路径达到最优值,甚至有时候会出现参数在最小值附近徘徊而不是立即收敛。但当数据量很大的时候,随机梯度下降法经常优于批梯度下降法。
当J为凸函数时,梯度下降法相当于让参数 θ \theta θ不断向J的最小值位置移动
梯度下降法的缺陷:如果函数为非凸函数,有可能找到的并非全局最优值,而是局部最优值。
2、最小二乘法矩阵求解
令
X
=
[
(
x
(
1
)
)
T
(
x
(
2
)
)
T
…
(
x
(
n
)
)
T
]
X = \left[ \begin{array} {cccc} (x^{(1)})^T\\ (x^{(2)})^T\\ \ldots \\ (x^{(n)})^T \end{array} \right]
X=⎣⎢⎢⎡(x(1))T(x(2))T…(x(n))T⎦⎥⎥⎤
其中,
x ( i ) = [ x 1 ( i ) x 2 ( i ) … x d ( i ) ] x^{(i)} = \left[ \begin{array} {cccc} x_1^{(i)}\\ x_2^{(i)}\\ \ldots \\ x_d^{(i)} \end{array} \right] x(i)=⎣⎢⎢⎢⎡x1(i)x2(i)…xd(i)⎦⎥⎥⎥⎤
由于
Y = [ y ( 1 ) y ( 2 ) … y ( n ) ] Y = \left[ \begin{array} {cccc} y^{(1)}\\ y^{(2)}\\ \ldots \\ y^{(n)} \end{array} \right] Y=⎣⎢⎢⎡y(1)y(2)…y(n)⎦⎥⎥⎤
h θ ( x ) h_\theta(x) hθ(x)可以写作
h θ ( x ) = X θ h_\theta(x)=X\theta hθ(x)=Xθ
对于向量来说,有
z T z = ∑ i z i 2 z^Tz = \sum_i z_i^2 zTz=i∑zi2
因此可以把损失函数写作
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 ( θ ) J(\theta) J(θ),对 θ \theta θ求导可得:
令偏导数等于零,由于最后一项和 θ \theta θ无关,偏导数为0。
因此,
∂ J ( θ ) ∂ θ = 1 2 ∂ ∂ θ θ T X T X θ − ∂ ∂ θ θ T X T Y \frac{\partial{J(\theta)}}{\partial\theta} = \frac{1}{2}\frac{\partial}{\partial\theta} \theta^TX^TX\theta - \frac{\partial}{\partial\theta} \theta^T X^TY ∂θ∂J(θ)=21∂θ∂θTXTXθ−∂θ∂θTXTY
利用矩阵求导性质,
∂ x ⃗ T α ∂ x ⃗ = α \frac{\partial \vec x^T\alpha}{\partial \vec x} =\alpha ∂x∂xTα=α
和 和 和
∂ A T B ∂ x ⃗ = ∂ A T ∂ x ⃗ B + ∂ B T ∂ x ⃗ A \frac{\partial A^TB}{\partial \vec x} = \frac{\partial A^T}{\partial \vec x}B + \frac{\partial B^T}{\partial \vec x}A ∂x∂ATB=∂x∂ATB+∂x∂BTA
∂ J ( θ ) ∂ θ = X T X θ − X T Y \frac{\partial{J(\theta)}}{\partial\theta} = X^TX\theta - X^TY ∂θ∂J(θ)=XTXθ−XTY
令导数等于零,
X T X θ = X T Y X^TX\theta = X^TY XTXθ=XTY
θ = ( X T X ) ( − 1 ) X T Y \theta = (X^TX)^{(-1)}X^TY θ=(XTX)(−1)XTY
3、牛顿法
f ( θ ) ′ = f ( θ ) Δ , Δ = θ 0 − θ 1 f(\theta)' = \frac{f(\theta)}{\Delta},\Delta = \theta_0 - \theta_1 f(θ)′=Δf(θ),Δ=θ0−θ1
可 求 得 , θ 1 = θ 0 − f ( θ 0 ) f ( θ 0 ) ′ 可求得,\theta_1 = \theta_0 - \frac {f(\theta_0)}{f(\theta_0)'} 可求得,θ1=θ0−f(θ0)′f(θ0)
重复迭代,可以让逼近取到 f ( θ ) f(\theta) f(θ)的最小值
当我们对损失函数 l ( θ ) l(\theta) l(θ)进行优化的时候,实际上是想要取到 l ′ ( θ ) l'(\theta) l′(θ)的最小值,因此迭代公式为:
θ : = θ − l ′ ( θ ) l ′ ′ ( θ ) \theta :=\theta-\frac{l'(\theta)}{l''(\theta)} θ:=θ−l′′(θ)l′(θ)
当 θ 是 向 量 值 的 时 候 , θ : = θ − H − 1 Δ θ l ( θ ) 当\theta是向量值的时候,\theta :=\theta - H^{-1}\Delta_{\theta}l(\theta) 当θ是向量值的时候,θ:=θ−H−1Δθl(θ)
其中,
Δ
θ
l
(
θ
)
\Delta_{\theta}l(\theta)
Δθl(θ)是
l
(
θ
)
l(\theta)
l(θ)对
θ
i
\theta_i
θi的偏导数,
H
H
H是
J
(
θ
)
J(\theta)
J(θ)的海森矩阵,
H
i
j
=
∂
2
l
(
θ
)
∂
θ
i
∂
θ
j
H_{ij} = \frac{\partial ^2l(\theta)}{\partial\theta_i\partial\theta_j}
Hij=∂θi∂θj∂2l(θ)
问题:请用泰勒展开法推导牛顿法公式。
Answer:将 f ( x ) f(x) f(x)用泰勒公式展开到第二阶,
f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + 1 2 f ′ ′ ( x 0 ) ( x − x 0 ) 2 f(x) = f(x_0) + f'(x_0)(x - x_0)+\frac{1}{2}f''(x_0)(x - x_0)^2 f(x)=f(x0)+f′(x0)(x−x0)+21f′′(x0)(x−x0)2
对上式求导,并令导数等于0,求得x值
f ′ ( x ) = f ′ ( x 0 ) + f ′ ′ ( x 0 ) x − f ′ ′ ( x 0 ) x 0 = 0 f'(x) = f'(x_0) + f''(x_0)x -f''(x_0)x_0 = 0 f′(x)=f′(x0)+f′′(x0)x−f′′(x0)x0=0
可以求得,
x = x 0 − f ′ ( x 0 ) f ′ ′ ( x 0 ) x = x_0 - \frac{f'(x_0)}{f''(x_0)} x=x0−f′′(x0)f′(x0)
牛顿法的收敛速度非常快,但海森矩阵的计算较为复杂,尤其当参数的维度很多时,会耗费大量计算成本。我们可以用其他矩阵替代海森矩阵,用拟牛顿法进行估计,
4、拟牛顿法
拟牛顿法的思路是用一个矩阵替代计算复杂的海森矩阵H,因此要找到符合H性质的矩阵。
要求得海森矩阵符合的条件,同样对泰勒公式求导 f ′ ( x ) = f ′ ( x 0 ) + f ′ ′ ( x 0 ) x − f ′ ′ ( x 0 ) x 0 f'(x) = f'(x_0) + f''(x_0)x -f''(x_0)x_0 f′(x)=f′(x0)+f′′(x0)x−f′′(x0)x0
令 x = x 1 x = x_1 x=x1,即迭代后的值,代入可得:
f ′ ( x 1 ) = f ′ ( x 0 ) + f ′ ′ ( x 0 ) x 1 − f ′ ′ ( x 0 ) x 0 f'(x_1) = f'(x_0) + f''(x_0)x_1 - f''(x_0)x_0 f′(x1)=f′(x0)+f′′(x0)x1−f′′(x0)x0
更一般的,
f ′ ( x k + 1 ) = f ′ ( x k ) + f ′ ′ ( x k ) x k + 1 − f ′ ′ ( x k ) x k f'(x_{k+1}) = f'(x_k) + f''(x_k)x_{k+1} - f''(x_k)x_k f′(xk+1)=f′(xk)+f′′(xk)xk+1−f′′(xk)xk
f ′ ( x k + 1 ) − f ′ ( x k ) = f ′ ′ ( x k ) ( x k + 1 − x k ) = H ( x k + 1 − x k ) f'(x_{k+1}) - f'(x_k) = f''(x_k)(x_{k+1}- x_k)= H(x_{k+1}- x_k) f′(xk+1)−f′(xk)=f′′(xk)(xk+1−xk)=H(xk+1−xk)
x k x_k xk为第k个迭代值
即找到矩阵G,使得它符合上式。
常用的拟牛顿法的算法包括DFP,BFGS等.
4、线性回归的评价指标
均方误差(MSE): 1 m ∑ i = 1 m ( y ( i ) − y ^ ( i ) ) 2 \frac{1}{m}\sum^{m}_{i=1}(y^{(i)} - \hat y^{(i)})^2 m1∑i=1m(y(i)−y^(i))2
均方根误差(RMSE): M S E = 1 m ∑ i = 1 m ( y ( i ) − y ^ ( i ) ) 2 \sqrt{MSE} = \sqrt{\frac{1}{m}\sum^{m}_{i=1}(y^{(i)} - \hat y^{(i)})^2} MSE=m1∑i=1m(y(i)−y^(i))2
平均绝对误差(MAE):$\frac{1}{m}\sum^{m}_{i=1} | (y^{(i)} - \hat y^{(i)} | $
但以上评价指标都无法消除量纲不一致而导致的误差值差别大的问题,最常用的指标是 R 2 R^2 R2,可以避免量纲不一致问题
R 2 : = 1 − ∑ i = 1 m ( y ( i ) − y ^ ( i ) ) 2 ∑ i = 1 m ( y ˉ − y ^ ( i ) ) 2 = 1 − 1 m ∑ i = 1 m ( y ( i ) − y ^ ( i ) ) 2 1 m ∑ i = 1 m ( y ˉ − y ^ ( i ) ) 2 = 1 − M S E V A R R^2: = 1-\frac{\sum^{m}_{i=1}(y^{(i)} - \hat y^{(i)})^2}{\sum^{m}_{i=1}(\bar y - \hat y^{(i)})^2} =1-\frac{\frac{1}{m}\sum^{m}_{i=1}(y^{(i)} - \hat y^{(i)})^2}{\frac{1}{m}\sum^{m}_{i=1}(\bar y - \hat y^{(i)})^2} = 1-\frac{MSE}{VAR} R2:=1−∑i=1m(yˉ−y^(i))2∑i=1m(y(i)−y^(i))2=1−m1∑i=1m(yˉ−y^(i))2m1∑i=1m(y(i)−y^(i))2=1−VARMSE
可以把 R 2 R^2 R2理解为,回归模型可以成功解释的数据方差部分在数据固有方差中所占的比例, R 2 R^2 R2越接近1,表示可解释力度越大,模型拟合的效果越好。
5、sklearn.linear_model参数详解:
fit_intercept : 默认为True,是否计算该模型的截距。如果使用中心化的数据,可以考虑设置为False,不考虑截距。注意这里是考虑,一般还是要考虑截距
normalize: 默认为false. 当fit_intercept设置为false的时候,这个参数会被自动忽略。如果为True,回归器会标准化输入参数:减去平均值,并且除以相应的二范数。当然啦,在这里还是建议将标准化的工作放在训练模型之前。通过设置sklearn.preprocessing.StandardScaler来实现,而在此处设置为false
copy_X : 默认为True, 否则X会被改写
n_jobs: int 默认为1. 当-1时默认使用全部CPUs
可用属性:
coef_:训练后的输入端模型系数,如果label有两个,即y值有两列。那么是一个2D的array
intercept_: 截距
可用的methods:
fit(X,y,sample_weight=None):
X: array, 稀疏矩阵 [n_samples,n_features]
y: array [n_samples, n_targets]
sample_weight: 权重 array [n_samples]
在版本0.17后添加了sample_weight
get_params(deep=True): 返回对regressor 的设置值
predict(X): 预测 基于 R^2值
score: 评估
参考https://blog.csdn.net/weixin_39175124/article/details/79465558
生成数据
#生成数据
import numpy as np
#生成随机数
np.random.seed(1234)
x = np.random.rand(500,3)
#构建映射关系,模拟真实的数据待预测值,映射关系为y = 4.2 + 5.7*x1 + 10.8*x2,可自行设置值进行尝试
y = x.dot(np.array([4.2,5.7,10.8]))
1、先尝试调用sklearn的线性回归模型训练数据
import numpy as np
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
%matplotlib inline
# 调用模型
lr = LinearRegression(fit_intercept=True)
# 训练模型
lr.fit(x,y)
print("估计的参数值为:%s" %(lr.coef_))
# 计算R平方
print('R2:%s' %(lr.score(x,y)))
# 任意设定变量,预测目标值
x_test = np.array([2,4,5]).reshape(1,-1)
y_hat = lr.predict(x_test)
print("预测值为: %s" %(y_hat))
估计的参数值为:[ 4.2 5.7 10.8]
R2:1.0
预测值为: [85.2]
2、最小二乘法的矩阵求解
class LR_LS():
def __init__(self):
self.w = None
def fit(self, X, y):
# 最小二乘法矩阵求解
#============================= show me your code =======================
self.w = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
#============================= show me your code =======================
def predict(self, X):
# 用已经拟合的参数值预测新自变量
#============================= show me your code =======================
y_pred = X.dot(self.w)
#============================= show me your code =======================
return y_pred
if __name__ == "__main__":
lr_ls = LR_LS()
lr_ls.fit(x,y)
print("估计的参数值:%s" %(lr_ls.w))
x_test = np.array([2,4,5]).reshape(1,-1)
print("预测值为: %s" %(lr_ls.predict(x_test)))
估计的参数值:[ 4.2 5.7 10.8]
预测值为: [85.2]
3、梯度下降法
class LR_GD():
def __init__(self):
self.w = None
def fit(self,X,y,alpha=0.02,loss = 1e-10): # 设定步长为0.002,判断是否收敛的条件为1e-10
y = y.reshape(-1,1) #重塑y值的维度以便矩阵运算
[m,d] = np.shape(X) #自变量的维度
self.w = np.zeros((d)) #将参数的初始值定为0
tol = 1e5
#============================= show me your code =======================
while tol > loss:
h_f = X.dot(self.w).reshape(-1,1)
theta = self.w + alpha*np.mean(X*(y - h_f),axis=0) #计算迭代的参数值
tol = np.sum(np.abs(theta - self.w))
self.w = theta
#============================= show me your code =======================
def predict(self, X):
# 用已经拟合的参数值预测新自变量
y_pred = X.dot(self.w)
return y_pred
if __name__ == "__main__":
lr_gd = LR_GD()
lr_gd.fit(x,y)
print("估计的参数值为:%s" %(lr_gd.w))
x_test = np.array([2,4,5]).reshape(1,-1)
print("预测值为:%s" %(lr_gd.predict(x_test)))
估计的参数值为:[ 4.20000001 5.70000003 10.79999997]
预测值为:[85.19999995]
参考
吴恩达 CS229课程
周志华 《机器学习》
李航 《统计学习方法》
https://hangzhou.anjuke.com/
https://www.jianshu.com/p/e0eb4f4ccf3e
https://blog.csdn.net/qq_28448117/article/details/79199835
https://blog.csdn.net/weixin_39175124/article/details/79465558n