线性回归

来源：互联网发布：淘宝晒图内裤爆光图片编辑：程序博客网时间：2024/05/22 03:33

1 线性回归理论

1.1 问题描述

对于只有一个自变量x和一个因变量y的数据集，我们假设(x,y)之间的关系可以用

h (x) = θ x + b

进行映射，那么我们就把求解这个函数的过程称为线性回归求解。可以在二维坐标系拟合出一条直线来代表(x,y)的趋势。很显然，不可能所有的点刚好在这条直线上，所以求解最优的标准就是，使点集(x,y)与拟合函数间的误差（平方误差）最小。

推广到多个自变量，假设一个数据集有n个自变量，也可以说成n个属性，x=(x1;x2;...;xn)，如果使用线性模型去拟合则有

h θ (x) = θ 1 x 1 + θ 2 x 2 + . . . + θ n x n + b

这种方式称为多元线性回归（multivariate linear regression）.需要求解的未知量有

θ1,θ2,...,θn,b.为了便于讨论，令

b=θ0x0，其中

x0=1，那么，原来的目标函数就可以表示为

h θ (x i) = θ 0 x 0 + θ 1 x 1 + θ 2 x 2 + . . . + θ n x n

用矩阵的方式表示为

h θ (x i) = θ T x i, i = 0, 1, . . ., n

即

y = θ T X

其中

X 是一个m*(n+1)的矩阵，前n维数据是数据集的数据，最后一维值为1，代表

X = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ x 11 x 21 ⋮ x m 1 x 12 x 22 ⋮ x m 2 \dots \dots ⋱ \dots x 1 n x 2 n ⋮ x m n 11 ⋮ 1 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟

所以说，求解线性模型的解，就是求解

θ。当

∑mi=1(hθ(x(i))−y(i))2最小，即平方误差最小时，有最优解。求解方法有梯度下降法（Gradient Descent）、Normal Equation等等。梯度下降有如下特点：需要预先选定步长a、需要多次迭代、特征值需要Scaling（统一到同一个尺度范围）。因此比较复杂，还有一种不需要迭代的求解方式-Normal Equation，简单、方便、不需要Feature Scaling。但是Normal Equation最小二乘法需要计算X的转置与逆矩阵，计算量很大，因此特征个数多时计算会很慢，只适用于特征个数小于100000时使用；当特征数量大于100000时使用梯度法。另外，当X不可逆时就有岭回归算法的用武之地了。

1.2 Normal Equation

Normal Equation最常用的就是最小二乘法。

对θ矩阵求导，得到其最优解为

θ^* = (X T X) - 1 X T y

1.3 Gradient Descent

梯度下降法（Gradient Descent）需要根据平方误差，定义该线性回归模型的损失函数（Cost Function）：

J (θ) = J (θ 0, θ 1, . . ., θ n) = 1 2 m \sum i = 1 m (h θ (x (i)) - y (i)) 2

线性回归的损耗函数的值与回归系数θ的关系是碗状的，只有一个最小点。线性回归的求解过程如同Logistic回归，区别在于学习模型函数hθ(x)不同.

1.4 局部加权线性回归

线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有最小均方误差的无偏估计。显而易见，如果模型欠拟合将不能取得最好的预测效果。所以有些方法允许在估计中引人一些偏差，从而降低预测的均方误差。其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression, LWLR )。在该算法中，我们给待预测点附近的每个点赋予一定的权重.于是公式变为:

θ^= (X T W X) - 1 X T W y

W是（m,m）矩阵，m表示样本数。

LWLR使用 “核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：

w (i, i) = e x p (| x ( i ) - x | - 2 k 2)

k需要优化选择.

局部加权线性回归也存在一个问题，即增加了计算量，因为它对每个点做预测时都必须使用整个数据集，而不是计算出回归系数得到回归方程后代入计算即可。因此该算法不被推荐。

2 简单的Python实例

根据线性回归的理论的Normal Equation方法，使用Python进行简单的求解，并把结果可视化。

import numpy as npimport pandas as pdimport seaborn as snsimport matplotlib.pyplot as plt%matplotlib inline

# 使用Pandas读取数据，并查看几条样例数据。# 第一列全是0，根据上面的理论知识，这是x0的占位符，不要太在意。我们需要拟合的是x1与y之间的线性关系。ex0 = pd.read_csv("ex0.txt", sep="\t", names=["x0", "x1", "y"])ex0.head()

x0 x1 y 0 1.0 0.067732 3.176513 1 1.0 0.427810 3.816464 2 1.0 0.995731 4.550095 3 1.0 0.738336 4.256571 4 1.0 0.981083 4.560815

2.1 整体拟合效果预览

# seaborn可以直接拟合出直线，我们先看一下最终的可视化效果。# 下图的圆点就是原始数据(x1,y),直线就是数据集的回归线sns.lmplot(x = 'x1', y='y', data=ex0)

<seaborn.axisgrid.FacetGrid at 0x7fa34ec26438>

seaborn线性拟合

2.2 逐步求解线性回归参数

下面根据理论知识，自己动手，逐步求解回归线的参数θ.

根据上面推导的多元线性回归最优解公式

θ^* = (X T X) - 1 X T y

其中矩阵X和矩阵y可从刚刚读取的ex0变量中获得。

# 获取矩阵X和矩阵y# 为了使用NumPy中矩阵计算的方法，将Pandas的DataFrame对象转换成NumPy中的ndarray对象# 如果读取的数据文件没有这里的x0这个数据，需要手动添加一列全为1的数据到xMat对象中xMat = ex0.as_matrix(['x1', 'x0'])yMat = ex0.as_matrix(['y'])

# 计算XTX，注意矩阵乘法是np.dot()，而不是np.mul()xTx = np.dot(xMat.T, xMat)# 由于有矩阵的逆的计算，所以需要判断行列式不为零np.linalg.det(xTx) == 0

False

# 完成最后的theta值的计算ws = np.dot(np.linalg.inv(xTx), xMat.T)ws = np.dot(np.linalg.inv(xTx), np.dot(xMat.T, yMat))ws

array([[ 1.69532264],       [ 3.00774324]])

于是，我们得到解

f (x) = 1.69532264 x 1 + 3.00774324

注意，这里ws变量中哪个是x0的系数，根据矩阵运算的过程，这个系数在ws中的位置和x0在xMat变量的位置一致，即最后一位。

根据这个解，画出X对应的最佳拟合直线图

# 直线上的y值，即预测值yHat = np.dot(xMat, ws)# 画原始值的散点图，即(x1,y)fig = plt.figure()ax = fig.add_subplot(111)ax.scatter(ex0['x1'], yMat)# 画拟合的直线ax.plot(xMat[:,0], yHat)plt.show()

线性拟合效果

2.3 评估拟合效果

可以发现，不论什么样的数据集，都可以使用线性回归来拟合，那么拟合的效果如何呢。可以使用皮尔森相关系数（Pearson correlation coefficient）也称皮尔森积矩相关系数(Pearson product-moment correlation coefficient)来衡量，值在[-1,1],表示数据的线性相关程度。

# 计算y的估计值和实际值的相关系数，可以反应这套数据使用线性回归模型拟合的效果。# 两个向量都必须是行向量np.corrcoef(yHat.T, yMat.T)

array([[ 1.        ,  0.98647356],       [ 0.98647356,  1.        ]])

从结果可以看出，变量和自己的相关度都是1，表示完全呈线性相关，这是显然的。估计值和实际值线性相关度为0.98，说明(x1,y)呈线性正相关，线性拟合效果不错。

2.4 线性回归工具

在开发中，我们往往不必自己从零实现线性回归求解过程，可以使用SciPy中的linear_model包，或者StatsModels包，下面以Scipy为例。

from sklearn.linear_model import LinearRegressionregr = LinearRegression()# 线性拟合，fix(x,y)，xy都必须是列向量regr.fit(xMat[:, 0].reshape(200,1), yMat)predictions={}# 根据拟合的模型计算预测值predictions['value'] = regr.predict(xMat[:, 0].reshape(200,1))predictions['coefficient'] = regr.coef_predictions['intercept'] = regr.intercept_predictions

{'coefficient': array([[ 1.69532264]]), 'intercept': array([ 3.00774324]), 'value': array([[ 3.12257084],        [ 3.73301922],        [ 4.69582855],        ...        [ 3.20467701]])}

可以看出coefficient就是上面我们自己求解的系数x1，intercept就是θ0，value是使用求解的函数预测的y值。结果和上面自己实现的方法一致。

阅读全文

0 0