1、简单线性回归
所要求得的结果是一个具体的数值,而不是一个类别的话,则该问题是回归问题。只有一个特征的回归问题,成为简单线性回归。两个变量之间存在一次方函数关系,就称它们之间存在线性关系。拟合就是把平面上一系列的点,用一条光滑的曲线连接起来。
使用方程
y
=
a
x
+
b
y=a x+b
y=ax+b 最大程度拟合样本点。图中每一个点对应了一个样本特征
x
i
x^{i}
xi,表示第
i
i
i个数据点,样本对应的结果值为
y
i
y^{i}
yi。
如找到了拟合程度最大的方程,就是直到a,b的值。将
x
i
x^{i}
xi带入找到的这个方程,对应的结果值为:
y
^
i
\hat{y}^{i}
y^i。即:
y
^
(
i
)
=
a
x
(
i
)
+
b
\hat{y}^{(i)}=a x^{(i)}+b
y^(i)=ax(i)+b。
找到这个最佳的a和b,等价于使
∣
y
(
i
)
−
y
^
(
i
)
∣
\left|y^{(i)}-\hat{y}^{(i)}\right|
∣∣y(i)−y^(i)∣∣最小,即,使
(
y
(
i
)
−
y
^
(
i
)
)
2
\left(y^{(i)}-\hat{y}^{(i)}\right)^{2}
(y(i)−y^(i))2最小,等价于使
∑
i
=
1
m
(
y
i
−
y
^
i
)
2
\sum_{i=1}^{m}\left(y^{i}-\hat{y}^{i}\right)^{2}
∑i=1m(yi−y^i)2最小,等价于使
∑
i
=
1
m
(
y
i
−
a
x
i
−
b
)
2
\sum_{i=1}^{m}\left(y^{i}-a x^{i}-b\right)^{2}
∑i=1m(yi−axi−b)2最小。
这里我们称
∑
i
=
1
m
(
y
i
−
a
x
i
−
b
)
2
\sum_{i=1}^{m}\left(y^{i}-a x^{i}-b\right)^{2}
∑i=1m(yi−axi−b)2这个函数为损失函数。这里的损失可以理解为误差或错误。损失越小也就是错误越小。
2、使用最小二乘法求a,b
令损失函数为:
J
(
a
,
b
)
=
∑
i
=
1
m
(
y
(
i
)
−
a
x
(
i
)
−
b
)
2
J(a, b)=\sum_{i=1}^{m}\left(y^{(i)}-a x^{(i)}-b\right)^{2}
J(a,b)=∑i=1m(y(i)−ax(i)−b)2
问题转化为了求损失函数的最值问题。高等数学的求偏导:
令: ∂ J ( a , b ) ∂ a = 0 \frac{\partial J(a, b)}{\partial a}=0 ∂a∂J(a,b)=0 和 ∂ J ( a , b ) ∂ b = 0 \frac{\partial J(a, b)}{\partial b}=0 ∂b∂J(a,b)=0
对b求偏导:
因为:
∂
J
(
a
,
b
)
∂
b
\frac{\partial J(a, b)}{\partial b}
∂b∂J(a,b)=
2
∑
i
=
1
m
(
y
i
−
a
x
i
−
b
)
(
−
1
)
2 \sum_{i=1}^{m}\left(y^{i}-a x^{i}-b\right)(-1)
2∑i=1m(yi−axi−b)(−1)
所以: 2 ∑ i = 1 m ( y i − a x i − b ) ( − 1 ) = 0 2 \sum_{i=1}^{m}\left(y^{i}-a x^{i}-b\right)(-1)=0 2∑i=1m(yi−axi−b)(−1)=0
即: ∑ i = 1 m ( y i − a x i − b ) = 0 \sum_{i=1}^{m}\left(y^{i}-a x^{i}-b\right)=0 ∑i=1m(yi−axi−b)=0
即: ∑ i = 1 m y i − ∑ i = 1 m a x i − m b = 0 \sum_{i=1}^{m} y^{i}-\sum_{i=1}^{m} a x^{i}-m b=0 ∑i=1myi−∑i=1maxi−mb=0
即: m b = ∑ i = 1 m y i − ∑ i = 1 m a x i m b=\sum_{i=1}^{m} y^{i}-\sum_{i=1}^{m} a x^{i} mb=∑i=1myi−∑i=1maxi
即:
b
=
y
ˉ
−
a
x
ˉ
b=\bar{y}-a \bar{x}
b=yˉ−axˉ
对a求偏导:
因为:
∂
J
(
a
,
b
)
∂
a
\frac{\partial J(a, b)}{\partial a}
∂a∂J(a,b)=
2
∑
i
=
1
m
(
y
i
−
a
x
i
−
b
)
(
−
x
i
)
2 \sum_{i=1}^{m}\left(y^{i}-a x^{i}-b\right)\left(-x^{i}\right)
2∑i=1m(yi−axi−b)(−xi)
所以: ∑ i = 1 m ( y ( i ) − a x ( i ) − b ) x ( i ) = 0 \sum_{i=1}^{m}\left(y^{(i)}-a x^{(i)}-b\right) x^{(i)}=0 ∑i=1m(y(i)−ax(i)−b)x(i)=0
即: ∑ i = 1 m ( y ( i ) − a x ( i ) − y ˉ + a x ˉ ) x ( i ) = 0 \sum_{i=1}^{m}\left(y^{(i)}-a x^{(i)}-\bar{y}+a \bar{x}\right) x^{(i)}=0 ∑i=1m(y(i)−ax(i)−yˉ+axˉ)x(i)=0
此时上式仅有一个未知数a。整理后有:
∑ i = 1 m ( x ( i ) y ( i ) − x ( i ) y ˉ − a ( x ( i ) ) 2 + a x ˉ x ( i ) ) = 0 \sum_{i=1}^{m}\left(x^{(i)} y^{(i)}-x^{(i)} \bar{y}-a\left(x^{(i)}\right)^{2}+a \bar{x} x^{(i)}\right)=0 ∑i=1m(x(i)y(i)−x(i)yˉ−a(x(i))2+axˉx(i))=0
∑ i = 1 m ( x ( i ) y ( i ) − x ( i ) y ˉ ) − ∑ i = 1 m ( a ( x ( i ) ) 2 − a x ˉ x ( i ) ) = 0 \sum_{i=1}^{m}\left(x^{(i)} y^{(i)}-x^{(i)} \bar{y}\right)-\sum_{i=1}^{m}\left(a\left(x^{(i)}\right)^{2}-a \bar{x} x^{(i)}\right)=0 ∑i=1m(x(i)y(i)−x(i)yˉ)−∑i=1m(a(x(i))2−axˉx(i))=0
即:
a
=
∑
i
=
1
m
(
x
i
y
i
−
x
i
y
ˉ
)
∑
i
=
1
m
(
(
x
i
)
2
−
x
ˉ
x
i
)
a=\frac{\sum_{i=1}^{m}\left(x^{i} y^{i}-x^{i} \bar{y}\right)}{\sum_{i=1}^{m}\left(\left(x^{i}\right)^{2}-\bar{x} x^{i}\right)}
a=∑i=1m((xi)2−xˉxi)∑i=1m(xiyi−xiyˉ)
其中: ∑ i = 1 m x i y ˉ = y ˉ ∑ i = 1 m x i = m y ˉ 1 m ∑ i = 1 m x i = m y ˉ ⋅ x ˉ \sum_{i=1}^{m} x^{i} \bar{y}=\bar{y} \sum_{i=1}^{m} x^{i}=m \bar{y} \frac{1}{m} \sum_{i=1}^{m} x^{i}=m \bar{y} \cdot \bar{x} ∑i=1mxiyˉ=yˉ∑i=1mxi=myˉm1∑i=1mxi=myˉ⋅xˉ
则: ∑ i = 1 m x i y ˉ = m y ˉ ⋅ x ˉ = x ˉ ∑ i = 1 m y i = ∑ i = 1 m x ˉ y i \sum_{i=1}^{m} x^{i} \bar{y}=m \bar{y} \cdot \bar{x}=\bar{x} \sum_{i=1}^{m} y^{i}=\sum_{i=1}^{m} \bar{x} y^{i} ∑i=1mxiyˉ=myˉ⋅xˉ=xˉ∑i=1myi=∑i=1mxˉyi
即:
∑
i
=
1
m
x
i
y
ˉ
=
∑
i
=
1
m
x
ˉ
y
i
\sum_{i=1}^{m} x^{i} \bar{y}=\sum_{i=1}^{m} \bar{x} y^{i}
∑i=1mxiyˉ=∑i=1mxˉyi
那么:
a
=
∑
i
=
1
m
(
x
i
y
i
−
x
i
y
ˉ
)
∑
i
=
1
m
(
(
x
i
)
2
−
x
ˉ
x
i
)
=
∑
i
=
1
m
(
x
i
y
i
−
x
i
y
ˉ
−
x
ˉ
y
i
+
x
ˉ
⋅
y
ˉ
)
∑
i
=
1
m
(
(
x
i
)
2
−
x
ˉ
x
i
−
x
ˉ
x
i
+
x
ˉ
2
)
\begin{aligned} a &=\frac{\sum_{i=1}^{m}\left(x^{i} y^{i}-x^{i} \bar{y}\right)}{\sum_{i=1}^{m}\left(\left(x^{i}\right)^{2}-\bar{x} x^{i}\right)} \\ &=\frac{\sum_{i=1}^{m}\left(x^{i} y^{i}-x^{i} \bar{y}-\bar{x} y^{i}+\bar{x} \cdot \bar{y}\right)}{\sum_{i=1}^{m}\left(\left(x^{i}\right)^{2}-\bar{x} x^{i}-\bar{x} x^{i}+\bar{x}^{2}\right)} \end{aligned}
a=∑i=1m((xi)2−xˉxi)∑i=1m(xiyi−xiyˉ)=∑i=1m((xi)2−xˉxi−xˉxi+xˉ2)∑i=1m(xiyi−xiyˉ−xˉyi+xˉ⋅yˉ)
则:
a
=
∑
i
=
1
m
(
x
i
y
i
−
x
2
y
ˉ
)
∑
i
=
1
m
(
(
x
i
)
2
−
x
ˉ
x
i
)
=
∑
i
=
1
m
(
x
i
y
i
−
x
i
y
ˉ
−
x
ˉ
y
i
+
x
ˉ
⋅
y
ˉ
)
∑
i
=
1
m
(
(
x
i
)
2
−
2
x
ˉ
x
i
+
x
ˉ
2
)
=
∑
i
=
1
m
(
x
i
−
x
ˉ
)
(
y
i
−
y
ˉ
)
(
x
i
−
x
ˉ
)
2
\begin{aligned} a &=\frac{\sum_{i=1}^{m}\left(x^{i} y^{i}-x^{2} \bar{y}\right)}{\sum_{i=1}^{m}\left(\left(x^{i}\right)^{2}-\bar{x} x^{i}\right)} \\ &=\frac{\sum_{i=1}^{m}\left(x^{i} y^{i}-x^{i} \bar{y}-\bar{x} y^{i}+\bar{x} \cdot \bar{y}\right)}{\sum_{i=1}^{m}\left(\left(x^{i}\right)^{2}-2 \bar{x} x^{i}+\bar{x}^{2}\right)} \\ &=\frac{\sum_{i=1}^{m}\left(x^{i}-\bar{x}\right)\left(y^{i}-\bar{y}\right)}{\left(x^{i}-\bar{x}\right)^{2}} \end{aligned}
a=∑i=1m((xi)2−xˉxi)∑i=1m(xiyi−x2yˉ)=∑i=1m((xi)2−2xˉxi+xˉ2)∑i=1m(xiyi−xiyˉ−xˉyi+xˉ⋅yˉ)=(xi−xˉ)2∑i=1m(xi−xˉ)(yi−yˉ)
那么问题更加明确了:找到a、b的取值,使损失函数取最小值。这里:
a = ∑ i = 1 m ( x i − x ˉ ) ( y i − y ˉ ) ∑ i = 1 m ( x i − x ˉ ) 2 a=\frac{\sum_{i=1}^{m}\left(x^{i}-\bar{x}\right)\left(y^{i}-\bar{y}\right)}{\sum_{i=1}^{m}\left(x^{i}-\bar{x}\right)^{2}} a=∑i=1m(xi−xˉ)2∑i=1m(xi−xˉ)(yi−yˉ)
b
=
y
ˉ
−
a
x
ˉ
b=\bar{y}-a \bar{x}
b=yˉ−axˉ
代码实现:
# 最小二乘法实现的
import numpy as np
class SimpleLinearRegression1:
def __init__(self):
self.a_ = None
self.b_ = None
def fit(self,x_train,y_train):
x_mean = np.mean(x_train)
y_mean = np.mean(y_train)
num = 0.0
d = 0.0
# 先通过循环来求解
for x,y in zip (x_train,y_train):
num += (x - x_mean) * (y - y_mean)
d += (x - x_mean) ** 2
self.a_ = num / d
self.b_ = y_mean - self.a_ * x_mean
return self
def predict(self,x_predict):
return np.array([self._predict(x) for x in x_predict])
def _predict(self,x_single):
return self.a_ * x_single + self.b_
def __repr__(self):
return "SimpleLinearRegression1(a=%s,b=%s)" % (self.a_,self.b_)
其中加入for循环来计算效率太低了,其实可以使用numpy提供的向量运算来实现。
2.1 向量化运算
代码实现:
# 最小二乘法实现的
import numpy as np
class SimpleLinearRegression2:
def __init__(self):
self.a_ = None
self.b_ = None
def fit(self,x_train,y_train):
x_mean = np.mean(x_train)
y_mean = np.mean(y_train)
# 通过向量化运算
num = (x_train - x_mean).dot(y_train - y_mean) # 返回的是标量
d = (x_train - x_mean).dot(x_train - x_mean)
self.a_ = num / d
self.b_ = y_mean - self.a_ * x_mean
return self
def predict(self,x_predict):
return np.array([self._predict(x) for x in x_predict])
def _predict(self,x_single):
return self.a_ * x_single + self.b_
def __repr__(self):
return "SimpleLinearRegression2(a=%s,b=%s)" % (self.a_,self.b_)
3、评估线性回归的指标
在分类问题中可以使用准确度衡量,那在回归问题中使用什么衡量呢。
将数据分为训练数据和测试数据。利用训练集得到的a,b,将测试集的样本数据带入训练模型,也就是带入拟合后的方程,得到一个结果,对比模型得出的结果与真是的结果的差值。
在训练的时候:使
∑
i
=
1
m
(
y
train
(
i
)
−
a
x
train
(
i
)
−
b
)
2
\sum_{i=1}^{m}\left(y_{\text {train}}^{(i)}-a x_{\text {train}}^{(i)}-b\right)^{2}
∑i=1m(ytrain(i)−axtrain(i)−b)2尽可能小。训练结束后,得到a,b。即:得到相应的
y
^
test
(
i
)
=
a
x
test
(
i
)
+
b
\hat{y}_{\text {test}}^{(i)}=a x_{\text {test}}^{(i)}+b
y^test(i)=axtest(i)+b
3.1、均方误差(MSE,Mean Squared Error)
1 m ∑ i = 1 m ( y test ( i ) − y ^ test ( i ) ) 2 \frac{1}{m} \sum_{i=1}^{m}\left(y_{\text {test }}^{(i)}-\hat{y}_{\text {test }}^{(i)}\right)^{2} m1∑i=1m(ytest (i)−y^test (i))2
3.2、均方根误差(RMSE,Root Mean Squared Error)
1 m ∑ i = 1 m ( y test ( i ) − y ^ test ( i ) ) 2 = M S E test \sqrt{\frac{1}{m} \sum_{i=1}^{m}\left(y_{\text {test }}^{(i)}-\hat{y}_{\text {test }}^{(i)}\right)^{2}}=\sqrt{M S E_{\text {test }}} m1∑i=1m(ytest (i)−y^test (i))2=MSEtest
3.3、平均绝对误差(MAE,Mean Absolute Error)
1 m ∑ i = 1 m ∣ y test ( i ) − y ^ test ( i ) \frac{1}{m} \sum_{i=1}^{m} \mid y_{\text {test}}^{(i)}-\hat{y}_{\text {test}}^{(i)} m1∑i=1m∣ytest(i)−y^test(i)
3.4 判定系数 R-square
上面我们介绍了MSE,RMSE和MAE这三个指标。其实这三个指标还有一定的问题。我们评价分类问题的指标是分类的准确度:1表示最好,0表示最差。即使我们的分类问题不同,我们也可以用这个指标来表示不同分类算分的优劣。但是上面的这三个指标是没有这个性质的。这个问题的解决方法是使用判定系数R Square)这个指标。
对于分子来说,我们可以理解为使用我们的模型预测产生的错误;对于分母来说,可以理解为使用均值预测产生的错误。如果一个模型使用均值来进行预测,那么这个模型叫基准模型(Baseline Model),就是说我们的模型至少要比基准模型要好。如果整个式子的结果大于1,说明我们的模型甚至比不上基准模型。 结合整个式子来考虑的话,
R
2
{R^2}
R2越大越好,最大为 1 ,因为最好的模型就是使得分子为 0 ,整个式子为 1 。
这个式子还可以这样变形:
R
2
=
1
−
∑
i
(
y
^
i
−
y
i
)
2
∑
i
(
y
ˉ
−
y
i
)
2
=
1
−
(
∑
i
=
1
m
(
y
^
i
−
y
i
)
2
)
/
m
(
∑
i
=
1
m
(
y
i
−
y
ˉ
)
2
)
/
m
R^{2}=1-\frac{\sum_{i}\left(\hat{y}^{i}-y^{i}\right)^{2}}{\sum_{i}\left(\bar{y}-y^{i}\right)^{2}}=1-\frac{\left(\sum_{i=1}^{m}\left(\hat{y}^{i}-y^{i}\right)^{2}\right) / m}{\left(\sum_{i=1}^{m}\left(y^{i}-\bar{y}\right)^{2}\right) / m}
R2=1−∑i(yˉ−yi)2∑i(y^i−yi)2=1−(∑i=1m(yi−yˉ)2)/m(∑i=1m(y^i−yi)2)/m
则:
R
2
=
1
−
M
S
E
(
y
^
,
y
)
/
V
a
r
(
y
)
{R^2}=1-MSE(\hat{y},y)/Var(y)
R2=1−MSE(y^,y)/Var(y),分子变成了MSE,分母变成了y的方差。
代码实现:
def mean_squared_error(y_true,y_predict):
return np.sum((y_true - y_predict)**2)/ len(y_true)
def root_mean_squared_error(y_true,y_predict):
return sqrt(mean_squared_error(y_true,y_predict))
def mean_absolute_error(y_true,y_predict):
return np.sum(np.absolute(y_true - y_predict)) / len(y_true)
def r_squire(y_test,y_predict):
return 1-mean_squared_error(y_test,y_predict) / np.var(y_test)
最后,在自定义的简单线性回归中添加评分功能:
# 最小二乘法实现的
import numpy as np
class SimpleLinearRegression2:
def __init__(self):
self.a_ = None
self.b_ = None
def fit(self,x_train,y_train):
x_mean = np.mean(x_train)
y_mean = np.mean(y_train)
# 通过向量化运算
num = (x_train - x_mean).dot(y_train - y_mean) #dot是向量点乘,返回的是标量
d = (x_train - x_mean).dot(x_train - x_mean)
self.a_ = num / d
self.b_ = y_mean - self.a_ * x_mean
return self
def predict(self,x_predict):
return np.array([self._predict(x) for x in x_predict])
def _predict(self,x_single):
return self.a_ * x_single + self.b_
def score(self,x_test,y_test):
y_predict = self.predict(x_test)
return mean_squared_error(y_test,y_predict) / np.var(y_test)
def __repr__(self):
return "SimpleLinearRegression2(a=%s,b=%s)" % (self.a_,self.b_)
4、多元线性回归
其中
θ
0
\theta_{0}
θ0是截距,其余
θ
i
\theta_{i}
θi为各个系数。
问题转化为:求得一个
θ
{\theta}
θ,使目标尽可能小。下边的式子称为多元线性回归的正规方程解(Normal Equation)
代码实现:
import numpy as np
from sklearn.metrics import r2_score
class LinearRegression:
def __init(self):
self.coef_ = None # 系数
self.interception_ = None # 截距
self._theta = None
def fit_normal(self, X_train, y_train):
X_b = np.hstack([np.ones((len(X_train), 1)), X_train]) # 构造X_b X_train加上 虚构的都等于1的列
self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train) # 通过正规方程解求得theta
# 分开保存
self.interception_ = self._theta[0]
self.coef_ = self._theta[1:]
return self
def predict(self, X_predict):
X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
return X_b.dot(self._theta)
def score(self, X_test, y_test):
y_predict = self.predict(X_test)
return r2_score(y_test, y_predict)
def __repr__(self):
return "LinearRegression(coef_=%s,interception_=%s)" % (self.coef_, self.interception_)
5、总结
线性回归算法只能解决回归问题,它是一个典型的参数学习问题。线性回归算法的特点就是对数据有很强的解释性。
我们在使用线性回归算法时,是有一个假设的。假设就是数据和最终的输出结果之间有一定的线性关系。
通过正规方程解来实现的线性回归算法, 如果我们的样本数量或特征数量巨大的话,使用这种方法是不合适的。梯度下降法是求解最优模型的通用方法。