本文已参与「新人创作礼」活动，一起开启掘金创作之路
推荐b站视频：www.bilibili.com/video/BV1Ht…
回归分析(Regression toward mediocrity)-趋中回归
有人可能会好奇，为什么叫“回归”这个名称，它有什么具体含义？实际上，回归这种现象最早由英国生物统计学家高尔顿在研究父母亲和子女的遗传特性时所发现的一种有趣的现象：
身高这种遗传特性表现出“高个子父母，其子代身高也高于平均身高；但不见得比其父母更高，到一定程度后会往平均身高方向发生‘回归’”。 这种效应被称为“趋中回归”。
回归分析最早是19世纪末期高尔顿（Sir Francis Galton）所发展。高尔顿是生物统计学派的奠基人，他的表哥达尔文的巨著《物种起源》问世以后，触动他用统计方法研究智力进化问题，统计学上的“相关”和“回归”的概念也是高尔顿第一次使用的。
1855年，他发表了一篇“遗传的身高向平均数方向的回归”文章，分析儿童身高与父母身高之间的关系，发现父母的身高可以预测子女的身高，当父母越高或越矮时，子女的身高会比一般儿童高或矮，他将儿子与父母身高的这种现象拟合出一种线形关系。但是有趣的是：通过观察他注意到，尽管这是一种拟合较好的线形关系，但仍然存在例外现象：矮个的人的儿子比其父要高，身材较高的父母所生子女的身高将回降到人的平均身高。换句话说，当父母身高走向极端（或者非常高，或者非常矮）的人的子女，子女的身高不会象父母身高那样极端化，其身高要比父母们的身高更接近平均身高。高尔顿选用“回归”一词，把这一现象叫做“向平均数方向的 回归”（regression toward mediocrity）。
而关于父辈身高与子代身高的具体关系是如何的，高尔顿和他的学生K·Pearson观察了1078对夫妇，以每对夫妇的平均身高作为自变量，取他们的一个成年儿子的身高作为因变量，结果发现两者近乎一条直线，其回归 直线方程为：y^=33.73+0.516x ，这种趋势及回归方程表明父母身高每增加一个单位时，其成年儿子的身高平均增加0.516个单位。这样当然极端值就会向中心靠拢。
根据换算公式1英寸=0.0254米进行单位换算后可得： y=0.8567 + 0.516x
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
Test = pd.DataFrame({'父亲身高':[1.51, 1.64, 1.6, 1.73, 1.82, 1.87],
'孩子身高':[1.63, 1.7, 1.71, 1.72, 1.76, 1.86]})
Test = Test.sort_values('父亲身高')
通过观察我们发现随着父亲身高上涨孩子身高也随着上涨
plt.scatter(Test['父亲身高'], Test['孩子身高'])
我们可以通过公式构建一个身高的回归方程
from sklearn.linear_model import LinearRegression
ans = []
Y = Test['孩子身高']
X = Test['父亲身高'].values.reshape(-1,1)
model = LinearRegression()
model.fit(X, Y)
#模型预测
test_pred = model.predict(X)
ans.append(model.intercept_)
ans.append(model.coef_[0])
print('使得经验误差函数 RD(h) 取最小值的参数为：{}'.format(ans))
plt.scatter(Test['父亲身高'], Test['孩子身高'],label='真实值')
plt.plot(Test['父亲身高'], Test['孩子身高'])
plt.scatter(Test['父亲身高'], test_pred,label='预测值')
plt.plot(Test['父亲身高'], test_pred)
plt.legend()
使得经验误差函数 RD(h) 取最小值的参数为：[0.8585439999999991, 0.5141333333333337]
<matplotlib.legend.Legend at 0x11b04270588>
线性回归的数学推导主要涉及到以下几个知识点。
- 利用矩阵的知识对线性公式进行整合
- 误差项的分析
- 似然函数的理解
- 矩阵求偏导
- 线性回归的最终求解
利用矩阵的知识对线性公式进行整合
当只有一个变量是一元线性函数，但是日常生活中不止有一个变量，可能是多个变量起作用，所以我们需要将多个变量组合在一起，其中b为截距，为了方便计算可以去掉,也可以进行转化，将b转化为，其中为1，将多项求和的式子，转换成矩阵的乘法的表达式，因为矩阵计算快且方便。
注：行列式计算
为求和符号
为转置
二：误差项的分析
简单线性回归模型定义：
实际得到的可贷款金额和预估的可贷款金额之间是有一定误差的，这个就是误差项
思考: 当最小时是不是函数模拟效果最好？
所以我们构造目标函数：,并且需要取绝对值，不然当数据正负为0时，也是最小值， 所以转变为 也就是
我们只需要将y与x带入，就可以得到 与 的函数，如同所示
为截距
为斜率
凸函数意味着画出来看上去像山谷。且具有最小值，通过求偏导
from sympy import symbols, diff, solve
import numpy as np
# 数据集 D
X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])
# 构造经验误差函数
w, b = symbols('w b', real=True)
RDh = 0
for (xi, yi) in zip(X, y):
RDh += (yi - (xi*w + b))**2
RDh *= 1/len(X)
# 对 w 和 b 求偏导
eRDhw = diff(RDh, w)
eRDhb = diff(RDh, b)
# 求解方程组
ans = solve((eRDhw, eRDhb), (w, b))
print('使得经验误差函数 RD(h) 取最小值的参数为：{}'.format(ans))
使得经验误差函数 RD(h) 取最小值的参数为：{w: 0.514133333333440, b: 0.858543999999819}
注: w为 b为
Y = Test['孩子身高']
X = Test['父亲身高'].values.reshape(-1,1)
#模型预测
test_pred_jingyan = X*0.514133333333440 + 0.858543999999819
plt.scatter(Test['父亲身高'], Test['孩子身高'],label='真实值')
plt.plot(Test['父亲身高'], Test['孩子身高'])
plt.scatter(Test['父亲身高'], test_pred,label='预测值')
plt.plot(Test['父亲身高'], test_pred)
plt.scatter(Test['父亲身高'], test_pred_jingyan,label='预测值_经验')
plt.plot(Test['父亲身高'], test_pred)
plt.legend()
<matplotlib.legend.Legend at 0x11b043716d8>
因为线性回归的经验误差函数是平方之和，所以本节介绍的求解该经验误差函数的最小值的方法被称为 最小平方法 国内各种教材中也常称为 最小二乘法 。
三 ：似然函数的理解
补充：以概率的角度去解决问题
我们对勒让德的猜测，即最小二乘法，仍然抱有怀疑，万一这个猜测是错误的怎么办？
数学王子高斯（1777－1855）也像我们一样心存怀疑。 高斯换了一个思考框架，通过概率统计那一套来思考
每次的预测值与真实值x之间存在一个误差：
这些误差最终会形成一个概率分布，只是现在不知道误差的概率分布是什么。假设概率密度函数为：
再假设一个联合概率密度函数，这样方便把所有的所有数据利用起来
现在是不是需要将最大就是最优结果，就会想到最大似然估计。
最大似然估计
我们假设硬币有两面，一面是“花”，一面是“字”。 一般来说，我们都觉得硬币是公平的，也就是“花”和“字”出现的概率是差不多的。 如果我扔了100次硬币，100次出现的都是“花”。 在这样的事实下，我觉得似乎硬币的参数不是公平的。你硬要说是公平的，那就是侮辱我的智商。 这种通过事实，反过来猜测硬币的情况，就是似然。
通过事实，推断出最有可能的硬币情况，就是最大似然估计。
因为是关于的函数，并且也是一个概率密度函数,根据极大似然估计的思想，概率最大的最应该出现
当下面这个式子成立时，取得最大值：
然后高斯想，最小二乘法给出的答案是：
如果最小二乘法是对的，那么时应该取得最大值，即：
好，现在可以来解这个微分方程了。最终得到：
这是什么？这就是正态分布啊。 并且这还是一个充要条件：
也就是说，如果误差的分布是正态分布，那么最小二乘法得到的就是最有可能的值。
高斯分布的条件是独立且具有相同的分布，并且服从均值为0方差为的高斯分布
- 独立：张三的父亲是张三的父亲，李四的父亲是李四的父亲
- 同分布:衡量单位一样
- 高斯分布 绝大多数的情况下，在一个的空间内浮动不大 ，身高有一个范围
第二步，已经知道误差项是符合高斯分布的，所以误差项的概率值就是下面的式子。
其中
这只是一个误差，需要将所有误差结合起来，组成联合概率，联合概率等于各个概率乘，但是乘积的最大值无法获取，我们直接使用多个数相乘，转化成多个数相加的形式。取对数
注：
因为似然函数是越大越好，似然函数的值和对数似然函数的值是成正比的，对值求对数，并不会影响到最后求极限的值。所以才敢进行对数处理。
对上面的式子进行整合，得到
通过上面一系列推导，就把式子转化为最小二乘法的相关知识了
之后可以进行求偏导也可以使用梯度下降算法，寻找最优值
梯度下降算法
场景假设
个人被困在山上，需要从山上下来(找到山的最低点，也就是山谷)。但此时山上的浓雾很大，导致可视度很低；因此，下山的路径就无法确定，必须利用自己周围的信息一步一步地找到下山的路。这个时候，便可利用梯度下降算法来帮助自己下山。怎么做呢，首先以他当前的所处的位置为基准，寻找这个位置最陡峭的地方，然后朝着下降方向走一步，然后又继续以当前位置为基准，再找最陡峭的地方，再走直到最后到达最低处；同理上山也是如此，只是这时候就变成梯度上升算法了
首先，我们有一个可微分的函数。这个函数就代表着一座山。我们的目标就是找到这个函数的最小值，也就是山底。根据之前的场景假设，最快的下山的方式就是找到当前位置最陡峭的方向，然后沿着此方向向下走，对应到函数中，就是找到给定点的梯度 ，然后朝着梯度相反的方向，就能让函数值下降的最快！因为梯度的方向就是函数之变化最快的方向
当在最小值左边时，斜率为负数，向正方向移动，靠近最小值
当在最小值右边时，斜率为正数，向负方向移动，靠近最小值
注：
- 学习率 过大时，找不到最小值
- 学习率 过小时，速度太慢
- 当初始值 选取不合适时会陷入局部最优解
梯度下降法
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
二维问题的梯度下降法示例
"""
import math
import numpy as np
def func_2d(x):
"""
目标函数
:param x: 自变量，二维向量
:return: 因变量，标量
"""
return - math.exp(-(x[0] ** 2 + x[1] ** 2))
def grad_2d(x):
"""
目标函数的梯度
:param x: 自变量，二维向量
:return: 因变量，二维向量
"""
deriv0 = 2 * x[0] * math.exp(-(x[0] ** 2 + x[1] ** 2))
deriv1 = 2 * x[1] * math.exp(-(x[0] ** 2 + x[1] ** 2))
return np.array([deriv0, deriv1])
def gradient_descent_2d(grad, cur_x=np.array([0.1, 0.1]), learning_rate=0.01, precision=0.0001, max_iters=10000):
"""
二维问题的梯度下降法
:param grad: 目标函数的梯度
:param cur_x: 当前 x 值，通过参数可以提供初始值
:param learning_rate: 学习率，也相当于设置的步长
:param precision: 设置收敛精度
:param max_iters: 最大迭代次数
:return: 局部最小值 x*
"""
print(f"{cur_x} 作为初始值开始迭代...")
for i in range(max_iters):
grad_cur = grad(cur_x)
if np.linalg.norm(grad_cur, ord=2) < precision:
break # 当梯度趋近为 0 时，视为收敛
cur_x = cur_x - grad_cur * learning_rate
print("第", i, "次迭代：x 值为 ", cur_x)
print("局部最小值 x =", cur_x)
return cur_x
if __name__ == '__main__':
gradient_descent_2d(grad_2d, cur_x=np.array([1, -1]), learning_rate=0.2, precision=0.000001, max_iters=10000)
[ 1 -1] 作为初始值开始迭代...
第 0 次迭代：x 值为 [ 0.94586589 -0.94586589]
第 1 次迭代：x 值为 [ 0.88265443 -0.88265443]
第 2 次迭代：x 值为 [ 0.80832661 -0.80832661]
第 3 次迭代：x 值为 [ 0.72080448 -0.72080448]
第 4 次迭代：x 值为 [ 0.61880589 -0.61880589]
第 5 次迭代：x 值为 [ 0.50372222 -0.50372222]
第 6 次迭代：x 值为 [ 0.3824228 -0.3824228]
第 7 次迭代：x 值为 [ 0.26824673 -0.26824673]
第 8 次迭代：x 值为 [ 0.17532999 -0.17532999]
第 9 次迭代：x 值为 [ 0.10937992 -0.10937992]
第 10 次迭代：x 值为 [ 0.06666242 -0.06666242]
第 11 次迭代：x 值为 [ 0.04023339 -0.04023339]
第 12 次迭代：x 值为 [ 0.02419205 -0.02419205]
第 13 次迭代：x 值为 [ 0.01452655 -0.01452655]
第 14 次迭代：x 值为 [ 0.00871838 -0.00871838]
第 15 次迭代：x 值为 [ 0.00523156 -0.00523156]
第 16 次迭代：x 值为 [ 0.00313905 -0.00313905]
第 17 次迭代：x 值为 [ 0.00188346 -0.00188346]
第 18 次迭代：x 值为 [ 0.00113008 -0.00113008]
第 19 次迭代：x 值为 [ 0.00067805 -0.00067805]
第 20 次迭代：x 值为 [ 0.00040683 -0.00040683]
第 21 次迭代：x 值为 [ 0.0002441 -0.0002441]
第 22 次迭代：x 值为 [ 0.00014646 -0.00014646]
第 23 次迭代：x 值为 [ 8.78751305e-05 -8.78751305e-05]
第 24 次迭代：x 值为 [ 5.27250788e-05 -5.27250788e-05]
第 25 次迭代：x 值为 [ 3.16350474e-05 -3.16350474e-05]
第 26 次迭代：x 值为 [ 1.89810285e-05 -1.89810285e-05]
第 27 次迭代：x 值为 [ 1.13886171e-05 -1.13886171e-05]
第 28 次迭代：x 值为 [ 6.83317026e-06 -6.83317026e-06]
第 29 次迭代：x 值为 [ 4.09990215e-06 -4.09990215e-06]
第 30 次迭代：x 值为 [ 2.45994129e-06 -2.45994129e-06]
第 31 次迭代：x 值为 [ 1.47596478e-06 -1.47596478e-06]
第 32 次迭代：x 值为 [ 8.85578865e-07 -8.85578865e-07]
第 33 次迭代：x 值为 [ 5.31347319e-07 -5.31347319e-07]
第 34 次迭代：x 值为 [ 3.18808392e-07 -3.18808392e-07]
局部最小值 x = [ 3.18808392e-07 -3.18808392e-07]
线性回归代码实现
import sys
import os
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
#对数据集中的样本属性进行分割，制作X和Y矩阵
def feature_label_split(pd_data):
#行数、列数
row_cnt, column_cnt = pd_data.shape
#生成新的X、Y矩阵
X = np.empty([row_cnt, column_cnt-1]) #生成两个随机未初始化的矩阵
Y = np.empty([row_cnt, 1])
for i in range(0, row_cnt):
row_array = redwine_data.iloc[i, ]
X[i] = np.array(row_array[0:-1])
Y[i] = np.array(row_array[-1])
return X, Y
#把特征数据进行标准化为均匀分布
def uniform_norm(X_in):
X_max = X_in.max(axis=0)
X_min = X_in.min(axis=0)
X = (X_in-X_min)/(X_max-X_min)
return X
#线性回归模型
class linear_regression():
def fit(self, train_X_in, train_Y, learning_rate=0.03, lamda=0.03, regularization="l2"):
#样本个数、样本的属性个数
case_cnt, feature_cnt = train_X_in.shape
#X矩阵添加X0向量
train_X = np.c_[train_X_in, np.ones(case_cnt, )]
#初始化待调参数theta
self.theta = np.zeros([feature_cnt+1, 1])
max_iter_num = sys.maxsize #最多迭代次数
step = 0 #当前已经迭代的次数
pre_step = 0 #上一次得到较好学习误差的迭代学习次数
last_error_J = sys.maxsize #上一次得到较好学习误差的误差函数值
threshold_value = 1e-6 #定义在得到较好学习误差之后截止学习的阈值
stay_threshold_times = 10 #定义在得到较好学习误差之后截止学习之前的学习次数
for step in range(0, max_iter_num):
#预测值
pred = train_X.dot(self.theta)
#损失函数
J_theta = sum((pred-train_Y)**2) / (2*case_cnt)
#更新参数theta
self.theta -= learning_rate*(lamda*self.theta + (train_X.T.dot(pred-train_Y))/case_cnt)
#检测损失函数的变化值，提前结束迭代
if J_theta < last_error_J - threshold_value:
last_error_J = J_theta
pre_step = step
elif step - pre_step > stay_threshold_times:
break
#定期打印，方便用户观察变化
if step % 100 == 0:
print("step %s: %.6f" % (step, J_theta))
def predict(self, X_in):
case_cnt = X_in.shape[0]
X = np.c_[X_in, np.ones(case_cnt, )]
pred = X.dot(self.theta)
return pred
#主函数
if __name__ == "__main__":
#读取样本数据
redwine_data = pd.read_csv("winequality-red.csv", sep=";")
#样本数据进行X、Y矩阵分离
X, Y = feature_label_split(redwine_data)
#对X矩阵进行归一化
unif_X = uniform_norm(X)
#对样本数据进行训练集和测试集的划分
unif_trainX, unif_testX, train_Y, test_Y = train_test_split(unif_X, Y, test_size=0.3, random_state=0)
#模型训练
model = linear_regression()
model.fit(unif_trainX, train_Y, learning_rate=0.1)
#模型预测
test_pred = model.predict(unif_testX)
test_pred_error = sum((test_pred-test_Y)**2) / (2*unif_testX.shape[0])
print("Test error is %.6f" % (test_pred_error))
step 0: 16.308311
step 100: 0.330867
step 200: 0.306297
step 300: 0.295023
step 400: 0.289310
step 500: 0.286208
step 600: 0.284439
step 700: 0.283390
step 800: 0.282748
step 900: 0.282346
step 1000: 0.282087
step 1100: 0.281918
step 1200: 0.281806
step 1300: 0.281730
step 1400: 0.281679
step 1500: 0.281645
step 1600: 0.281621
step 1700: 0.281604
step 1800: 0.281593
Test error is 0.239596
官方文档解释
sklearn.linear_model.LinearRegression(*, fit_intercept=True, normalize=False, copy_X=True, n_jobs=None, positive=False)
fit_intercept
默认为True,当为True时计算截距,Flase为不计算,一般不修改
normalize
当FIT_CHECTECTITS设置为False时，忽略此参数。如果为True，则在回归之前，通过减去均值并用L2-范数除以回归变量X将被规范化。如果您希望标准化，请在调用FIT之前使用StandardScaler，并使用Normalize=false。一般不修改
copy_X
如果为True，将复制X；否则为X。 否则，它可能会被覆盖。一般不修改
n_jobs
None表示1，-1表示使用所有处理器 一般不修改
positive
设置为True时，强制系数为正
案例
from sklearn.linear_model import LinearRegression
from sklearn.metrics import explained_variance_score
redwine_data = pd.read_csv("winequality-red.csv", sep=";")
redwine_data.head()
| fixed acidity | volatile acidity | citric acid | residual sugar | chlorides | free sulfur dioxide | total sulfur dioxide | density | pH | sulphates | alcohol | quality | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 7.4 | 0.70 | 0.00 | 1.9 | 0.076 | 11.0 | 34.0 | 0.9978 | 3.51 | 0.56 | 9.4 | 5 |
| 1 | 7.8 | 0.88 | 0.00 | 2.6 | 0.098 | 25.0 | 67.0 | 0.9968 | 3.20 | 0.68 | 9.8 | 5 |
| 2 | 7.8 | 0.76 | 0.04 | 2.3 | 0.092 | 15.0 | 54.0 | 0.9970 | 3.26 | 0.65 | 9.8 | 5 |
| 3 | 11.2 | 0.28 | 0.56 | 1.9 | 0.075 | 17.0 | 60.0 | 0.9980 | 3.16 | 0.58 | 9.8 | 6 |
| 4 | 7.4 | 0.70 | 0.00 | 1.9 | 0.076 | 11.0 | 34.0 | 0.9978 | 3.51 | 0.56 | 9.4 | 5 |
Y = redwine_data['quality']
del redwine_data['quality']
unif_X = redwine_data
## 训练集划分
unif_trainX, unif_testX, train_Y, test_Y = train_test_split(unif_X, Y, test_size=0.3, random_state=0)
## 模型校验
model = LinearRegression()
model.fit(unif_trainX, train_Y)
#模型预测
test_pred = model.predict(unif_testX)
均误差方（MSE）
指标解释：所有样本的样本误差的平方的均值
指标解读：均误差方越接近0，模型越准确
平均绝对误差（MAE）
指标解释：所有样本的样本误差的绝对值的均值
指标解读：平均绝对误差的单位与因变量单位一致，越接近0，模型越准确
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
print(f'mean squared error is: {mean_squared_error(test_Y,test_pred)}')
print(f'mean absolute error is: {mean_absolute_error(test_Y,test_pred)}')
mean squared error is: 0.4007803663750049
mean absolute error is: 0.4871262164592796
补充
