本文传送机
用线性回归找到最佳拟合直线
线性回归
优点:结果易于理解,计算上不复杂
缺点:对非线性的数据拟合不好
适用数据类型:数值型和标称型数据
回归方程为:,其中为回归系数
回归的一般方法
- 收集数据:采用任意方法收集数据
- 准备数据:回归需要数值型数据,标称型数据将被转换为二值型数据
- 分析数据:绘制出数据的可视化二维图将有助于对数据做出理解和分析,在采用压缩减法求得新回归系数之后,可以将新拟合线在图上作为对比
- 训练算法:找到回归系数
- 测试算法:使用或者预测值和数据的拟合度来分析模型的效果
- 使用算法:使用回归,可以在给定输入的时候预测出一个数值,这是对分类方法的提升,因为这样可以预测连续型数值而不仅仅是离散的类别标签。
现在的问题是,给定和对应的,怎样才能找到呢?一个常用的方法就是找出使误差最小的 。这里的误差是指预测值和真实值之间的差值,使用丐武训哈的简单累加将使得正差值和负差值相互抵消,所以我们采用平方误差。
平方误差为:
写成矩阵运算的形式为:
为求的极小值,可对进行求导,的极值点一定出现在的位置:
可得
这里的便是回归系数的最优估计向量。
代码:
import numpy as np
import matplotlib.pyplot as plt
def load_dataset(filename):
feature_number = len(open(filename).readline().split('\t'))-1
data_mat = []
label_mat = []
file = open(filename)
for line in file.readlines():
line_array = []
curr_line = line.strip().split('\t')
for i in range(feature_number):
line_array.append(float(curr_line[i]))
data_mat.append(line_array)
label_mat.append(curr_line[-1])
return data_mat, label_mat
def linear_regress(x_array, y_array):
x_mat = np.mat(x_array, dtype=np.float64)
y_mat = np.mat(y_array, dtype=np.float64).T
xTx = x_mat.T * x_mat
if np.linalg.det(xTx) == 0.0:
print("This matrix is singular, cannot do inverse")
return
w_hat = xTx.I * (x_mat.T * y_mat)
return w_hat
if __name__ == '__main__':
x_array, y_array = load_dataset('ex0.txt')
w_hat = linear_regress(x_array, y_array)
x_mat = np.mat(x_array, dtype=np.float64)
y_mat = np.mat(y_array, dtype=np.float64)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x_mat[:, 1].flatten().A[0], y_mat.T[:, 0].flatten().A[0])
x_copy = x_mat.copy()
x_copy.sort(0)
y_hat = x_copy * w_hat
ax.plot(x_copy[:, 1], y_hat, c="red")
plt.show()
结果:
那么,如何判断这些模型的好坏呢?
有种方法可以计算预测序列和真实值之间的匹配程度,那就是计算这两个序列的相关系数(numpy中提供了corrcoef方法):
print(np.corrcoef(y_hat.T, y_mat))
结果:
[[1. 0.13653777]
[0.13653777 1. ]]
结果的含义是:和自己的匹配最完美,而和的相关系数为0.1365
局部加权线性回归
线性回归的一个问题是有可能出现欠拟合现象,因为它求得是具有最小均方误差得无偏估计。显而易见,如果模型欠拟合将不能取得最好得预测结果。所以有些方法允许在估计中引入一些偏差,从而降低预测得均方误差。
其中的一个方法是局部加权线性回归(Locally Weighted Linear Regression, LWLR),在该算法中,我们给待测点附近的每一个点赋予一定的权重,然后在这个子集上基于最小均方差来进行普通的回归。与kNN一样,这种算法每次预测均需要事先选取出对应的数据子集。该算法决出的回归系数的形式如下:
其中是一个矩阵,用来给每个数据点赋予权重。
LWLR(局部加权线性回归)使用“核”(与支持向量机中的核类似)来对附近点赋予更高的权重。核的类型可以自由选择,最常用的核就是高斯核,高斯核对应的权重如下:
这样就构建了一个只含有对角元素的权重矩阵,并且点与越近,将会越大。其中是一个需要用户指定的参数,它决定了对附近的点赋予多大的权重,这也是LWLR唯一需要考虑的参数。
代码:
import numpy as np
import matplotlib.pyplot as plt
def load_dataset(filename):
feature_number = len(open(filename).readline().split('\t'))-1
data_mat = []
label_mat = []
file = open(filename)
for line in file.readlines():
line_array = []
curr_line = line.strip().split('\t')
for i in range(feature_number):
line_array.append(float(curr_line[i]))
data_mat.append(line_array)
label_mat.append(curr_line[-1])
return data_mat, label_mat
def lwlr(test_point, x_array, y_array, k=1.0):
x_mat = np.mat(x_array, dtype=np.float64)
y_mat = np.mat(y_array, dtype=np.float64).T
m = np.shape(x_mat)[0]
weights = np.mat(np.eye(m))
for j in range(m):
diff_mat = test_point - x_mat[j, :]
weights[j, j] = np.exp(diff_mat * diff_mat.T / (-2.0 * k ** 2))
xTx = x_mat.T * (weights * x_mat)
if np.linalg.det(xTx) == 0.0:
print("This matrix is singular, cannot do inverse")
return
w_hat = xTx.I * (x_mat.T * weights * y_mat)
return test_point * w_hat
def lwlr_test(test_array, x_array, y_array, k=1.0):
m = np.shape(test_array)[0]
y_hat = np.zeros(m)
for i in range(m):
y_hat[i] = lwlr(test_array[i], x_array, y_array, k)
return y_hat
if __name__ == '__main__':
x_array, y_array = load_dataset('ex0.txt')
x_array = np.array(x_array, dtype=np.float64)
y_array = np.array(y_array, dtype=np.float64)
y_hat = lwlr_test(x_array, x_array, y_array, 0.003)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x_array[:, 1], y_array)
sorted_index = x_array[:, 1].argsort(0)
sorted_x_array = x_array[sorted_index]
sorted_y_hat = y_hat[sorted_index]
ax.plot(sorted_x_array[:, 1], sorted_y_hat, c="red")
plt.show()
结果:
通过缩减系数来“理解”数据
如果特征比样本点还多(n>m),也就是说输入数据的矩阵不是满秩矩阵。非满秩矩阵在求逆时会出现问题。
为了解决这个问题,统计学家引入了岭回归(ridge regression)的概念,同时还有lasso法(这种方法效果好但计算复杂);第二种缩减的方法为前向逐步法,可以得到和lasso差不多的效果,且更容易实现。
岭回归
简单来说,岭回归就是在矩阵上加一个从而使得矩阵变得非奇异,这样便能对求逆,其中是一个用户自定义参数。在这种情况下,回归系数的计算公式为:
岭回归最先用来处理特征数多余样本数的情况,现在也用于在估计中加入偏差,从而得到更好的估计。这里通过引入来限制所有的之和,通过引入该惩罚项,能够减少不重要的参数,这个技术在统计学中称为缩减(shrinkage)
缩减方法可以去掉不重要的参数,因此能更好得理解数据。此外,与简单的线性回归相比,缩减法能够取得更好的预测效果。
def ridge_regress(x_mat, y_mat, lam=0.2):
xTx = x_mat.T * x_mat
denom = xTx + np.eye(np.shape(x_mat)[1]) * lam
if np.linalg.det(denom) == 0.0:
print("This matrix is singular, cannot do inverse")
return
w_hat = denom.T * (x_mat.T * y_mat)
return w_hat
lasso
不难证明,在增加如下约束时,普通的最小二乘法回归会得到与岭回归的一样的公式:
上式限定了所有回归系数的平方和不能大于。使用普通的最小二乘法回归在当两个或更多的特征相关时,可能会得出一个很大的正系数和很大的负系数。正是因为上述限制条件的存在,使用岭回归可以避免这个问题。
与岭回归类似,另一个缩减方法lasso也对回归系数做出了限定,对应的约束条件如下:
唯一的不同点在于,这个约束条件使用绝对值取代了平方和。虽然约束形式只是稍作变化,结果却大项径庭:在足够小的时候,一些系数会因此被迫缩减到0,这个特性可以帮助我们更好地理解数据,这两个约束条件在公式上看起来相差无几,但细微得变化却极大地增加了计算复杂度(为了在这个新得约束条件下解出回归系数,需要使用二次规划算法)。下面将介绍一个更为简单的方法来得到结果,该方法叫做前向逐步回归。
前向逐步回归
前向逐步回归算法可以得到与lasso差不多的效果,但更加简单。它属于一种贪心算法,即每一步都尽可能减小误差。一开始,所有权重都设为1,然后每一步所做的决策是对某个权重增加或减少一个很小的值。
前向逐步回归伪代码
- 数据标准化,使其分布满足0均值和单位方差
- 在每轮迭代过程中:
- 设置当前最小误差lowest_error为正无穷
- 对每个特征:
- 增大或缩小:
- 改变一个系数得到一个新的
- 计算新下的误差error
- 如果误差error小于当前最小误差lowest_error:设置等于当前的
- 将设置为新的
代码:
import numpy as np
def load_dataset(filename):
feature_number = len(open(filename).readline().split('\t'))-1
data_array = []
label_array = []
file = open(filename)
for line in file.readlines():
line_array = []
curr_line = line.strip().split('\t')
for i in range(feature_number):
line_array.append(float(curr_line[i]))
data_array.append(line_array)
label_array.append(curr_line[-1])
return data_array, label_array
def regularize(x_mat):
in_mat = x_mat.copy()
in_means = np.mean(in_mat, 0)
in_var = np.var(in_mat, 0)
in_mat = (in_mat - in_means) / in_var
return in_mat
def calc_rss_error(y_array, y_hat_array):
return ((y_array - y_hat_array) ** 2).sum()
def stepwise(x_array, y_array, eps=0.01, iter_num=100):
x_mat = np.mat(x_array, dtype=np.float64)
y_mat = np.mat(y_array, dtype=np.float64).T
# 将特征值按照均值为0,方差为1进行标准化处理
y_mean = np.mean(y_mat, 0)
y_mat = y_mat - y_mean
x_mat = regularize(x_mat)
m, n = np.shape(x_mat)
return_mat = np.zeros((iter_num, n))
w_hat = np.zeros((n, 1))
w_hat_test = w_hat.copy()
w_hat_max = w_hat.copy()
for i in range(iter_num):
print("[iteration %d]" % i, w_hat.T)
lowest_error = np.inf
for j in range(n):
for sign in [-1, 1]:
w_hat_test = w_hat.copy()
w_hat_test[j] += eps * sign
y_test = x_mat * w_hat_test
rss_error = calc_rss_error(y_mat.A, y_test.A)
if rss_error < lowest_error:
lowest_error = rss_error
w_hat_max = w_hat_test
w_hat = w_hat_max.copy()
return_mat[i, :] = w_hat.T
return return_mat
if __name__ == '__main__':
x_array, y_array = load_dataset('abalone.txt')
w_hat = stepwise(x_array, y_array, 0.01, 200)
print("w_hat: ", w_hat)
CONTACT INFORMATION
E-mail: birdguan@seu.edu.cn
QQ: 46611253