线性回归理论参考
1.1、加载数据
# 加载数据集，
def loadDataSet(fileName): #general function to parse tab -delimited floats
numFeat = len(open(fileName).readline().split('\t')) - 1 #get number of fields
xArr = [] # x数据集
yArr = [] # y数据集
fr = open(fileName)
for line in fr.readlines():
lineArr =[]
curLine = line.strip().split('\t')
for i in range(numFeat):
lineArr.append(float(curLine[i]))
xArr.append(lineArr)
# 最后一列是y的值
yArr.append(float(curLine[-1]))
return xArr, yArr
1.1.1、样例数据中的第一列都是1.0, 即X0, 我们假定偏移量是一个常数， 第二列是下X1, 就是下图中的横坐标。
1.000000 0.067732 3.176513
1.000000 0.427810 3.816464
1.000000 0.995731 4.550095
1.000000 0.738336 4.256571
1.000000 0.981083 4.560815
1.000000 0.526171 3.929515
1.2、计算回归系数w
# 计算回归系数w
def standRegres(xArr,yArr):
'''
计算回归系数
:param xArr: x数据集
:param yArr: y数据集
:return: 回归系数
'''
xMat = mat(xArr)
yMat = mat(yArr).T # 由于yArr是一个列表， 而yMat需要的是一个列向量， 所以需要转置
xTx = xMat.T*xMat
# 前提条件， xTx不可逆
if linalg.det(xTx) == 0.0:
print("This matrix is singular, cannot do inverse")
return
ws = xTx.I * (xMat.T*yMat)
return ws
1.3、根据上文中推导的回归系数计算方法，求出回归系数向量，并根据回归系数向量绘制回归曲线，编写代码如下：
def plotRegression(xArr, yArr, ws):
"""
函数说明:绘制回归曲线和数据点
"""
xMat = np.mat(xArr) #创建xMat矩阵
yMat = np.mat(yArr) #创建yMat矩阵
xCopy = xMat.copy() #深拷贝xMat矩阵
xCopy.sort(0) #排序 如果直线的数据点次序混乱，绘图的时候会出现问题。所以先将点按照升序排列
yHat = xCopy * ws #计算对应的y值
fig = plt.figure()
ax = fig.add_subplot(111) #添加subplot
ax.plot(xCopy[:, 1], yHat, c = 'red') #绘制回归曲线
ax.scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue',alpha = .5) #绘制样本点
plt.title('DataSet') #绘制title
plt.xlabel('X')
plt.show()
if __name__ == '__main__':
# 加载数据集
xArr, yArr = loadDataSet('ex0.txt')
# 计算回归系数
ws = standRegres(xArr, yArr)
# 绘制回归曲线
plotRegression(xArr, yArr, ws)
1.4、如何判断拟合曲线的拟合效果的如何呢？当然，我们可以根据自己的经验进行观察，除此之外，我们还可以使用corrcoef方法，来比较预测值和真实值的相关性。
if __name__ == '__main__':
# plotDataSet()
# 加载数据集
xArr, yArr = loadDataSet('ex0.txt')
# 计算回归系数
ws = standRegres(xArr, yArr)
# 绘制回归曲线
plotRegression(xArr, yArr, ws)
# 使用corrcoef方法，来比较预测值和真实值的相关性。
xMat = np.mat(xArr) #创建xMat矩阵
yMat = np.mat(yArr) #创建yMat矩阵
yHat = xMat * ws
# 计算相关系数(需要保证两个向量都是行向量)
print(np.corrcoef(yHat.T, yMat))
结果如下：
可以看到，对角线上的数据是1.0，因为yMat和自己的匹配是完美的，而YHat和yMat的相关系数为0.98。
最佳拟合直线方法将数据视为直线进行建模，具有十分不错的表现。数据当中似乎还存在其他的潜在模式。那么如何才能利用这些模式呢？我们可以根据数据来局部调整预测，下面就会介绍这种方法。
二、 局部加权线性回归
线性回归的一个问题是有可能出现欠拟合现象，因为它求的是具有小均方误差的无偏估 计。显而易见，如果模型欠拟合将不能取得好的预测效果。所以有些方法允许在估计中引入一 些偏差，从而降低预测的均方误差。
其中的一个方法是局部加权线性回归（Locally Weighted Linear Regression，LWLR）。在该方法中，我们给待预测点附近的每个点赋予一定的权重。与kNN一样，这种算法每次预测均需要事先选取出对应的数据子集。该算法解除回归系数W的形式如下：
其中W是一个矩阵，这个公式跟我们上面推导的公式的区别就在于W，它用来给每个点赋予权重。
LWLR使用”核”（与支持向量机中的核类似）来对附近的点赋予更高的权重。核的类型可以自由选择，最常用的核就是高斯核，高斯核对应的权重如下：
2.1、局部加权回归系数计算
#局部加权线性回归
def lwlr(testPoint,xArr,yArr,k=1.0):
'''
:param testPoint: 测试样本点
:param xArr: x数据集
:param yArr: y数据集
:param k: 高斯核的k,自定义参数
:return: 回归系数
'''
xMat = mat(xArr); yMat = mat(yArr).T
m = shape(xMat)[0]
weights = mat(eye((m))) #创建权重对角矩阵， 是一个方阵， 阶数等于样本点个数
#遍历数据集 计算每个样本的权重
for j in range(m): # next 2 lines create weights matrix
diffMat = testPoint - xMat[j, :]
# 每个点高斯核对应的权重
weights[j, j] = exp(diffMat * diffMat.T / (-2.0 * k ** 2)) # 权重值大小以指数级别衰减
xTx = xMat.T * (weights * xMat)
if linalg.det(xTx) == 0.0:
print ("This matrix is singular, cannot do inverse")
return
ws = xTx.I * (xMat.T * (weights * yMat))
return testPoint * ws
2.2、 局部加权线性回归测试
def lwlrTest(testArr,xArr,yArr,k=1.0): #loops over all the data points and applies lwlr to each one
'''
:param testArr: 测试数据集x
:param xArr: x数据集
:param yArr: y数据集
:param k: 高斯核的k, 自定义参数
:return:
'''
m = shape(testArr)[0] # 测试数据集的大小
yHat = zeros(m)
for i in range(m):
yHat[i] = lwlr(testArr[i],xArr,yArr,k)
return yHat
2.3、绘制多条局部加权回归曲线
def plotlwlrRegression():
"""
函数说明:绘制多条局部加权回归曲线
"""
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
xArr, yArr = loadDataSet('ex0.txt') #加载数据集
yHat_1 = lwlrTest(xArr, xArr, yArr, 1.0) #根据局部加权线性回归计算yHat
yHat_2 = lwlrTest(xArr, xArr, yArr, 0.01) #根据局部加权线性回归计算yHat
yHat_3 = lwlrTest(xArr, xArr, yArr, 0.003) #根据局部加权线性回归计算yHat
xMat = np.mat(xArr) #创建xMat矩阵
yMat = np.mat(yArr) #创建yMat矩阵
srtInd = xMat[:, 1].argsort(0) #排序，返回索引值
xSort = xMat[srtInd][:,0,:]
fig, axs = plt.subplots(nrows=3, ncols=1,sharex=False, sharey=False, figsize=(10,8))
axs[0].plot(xSort[:, 1], yHat_1[srtInd], c = 'red') #绘制回归曲线
axs[1].plot(xSort[:, 1], yHat_2[srtInd], c = 'red') #绘制回归曲线
axs[2].plot(xSort[:, 1], yHat_3[srtInd], c = 'red') #绘制回归曲线
axs[0].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5) #绘制样本点
axs[1].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5) #绘制样本点
axs[2].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 20, c = 'blue', alpha = .5) #绘制样本点
#设置标题,x轴label,y轴label
axs0_title_text = axs[0].set_title(u'局部加权回归曲线,k=1.0',FontProperties=font)
axs1_title_text = axs[1].set_title(u'局部加权回归曲线,k=0.01',FontProperties=font)
axs2_title_text = axs[2].set_title(u'局部加权回归曲线,k=0.003',FontProperties=font)
plt.setp(axs0_title_text, size=8, weight='bold', color='red')
plt.setp(axs1_title_text, size=8, weight='bold', color='red')
plt.setp(axs2_title_text, size=8, weight='bold', color='red')
plt.xlabel('X')
plt.show()
if __name__ == '__main__':
plotlwlrRegression()
可以看到，当k越小，拟合效果越好。但是当k过小，会出现过拟合的情况，例如k等于0.003的时候。
三、 预测鲍鱼的年龄
接下来，我们将回归用于真实数据。在abalone.txt文件中记录了鲍鱼（一种水生物→__→）的年龄，鲍鱼年龄可以从鲍鱼壳的层数推算得到。
数据集是多维的，所以我们很难画出它的分布情况。每个维度数据的代表的含义没有给出，不过没有关系，我们只要知道最后一列的数据是y值就可以了，最后一列代表的是鲍鱼的真实年龄，前面几列的数据是一些鲍鱼的特征，例如鲍鱼壳的层数等。我们不做数据清理，直接用上所有特征，测试下我们的局部加权回归。
3.1、首先介绍下resError()， 用于评价误差大小
# 误差大小评价函数
def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays
'''
:param yArr: 真实数据
:param yHatArr: 预测数据
:return: 误差大小
'''
return ((yArr-yHatArr)**2).sum()
3.2、测试
def testAbaloneAge():
abX, abY = loadDataSet('abalone.txt')
print('训练集与测试集相同:局部加权线性回归,核k的大小对预测的影响:')
yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)
yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)
yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)
print('k=0.1时,误差大小为:', rssError(abY[0:99], yHat01.T))
print('k=1 时,误差大小为:', rssError(abY[0:99], yHat1.T))
print('k=10 时,误差大小为:', rssError(abY[0:99], yHat10.T))
print('')
print('训练集与测试集不同:局部加权线性回归,核k的大小是越小越好吗？更换数据集,测试结果如下:')
yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)
yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)
yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)
print('k=0.1时,误差大小为:', rssError(abY[100:199], yHat01.T))
print('k=1 时,误差大小为:', rssError(abY[100:199], yHat1.T))
print('k=10 时,误差大小为:', rssError(abY[100:199], yHat10.T))
print('')
print('训练集与测试集不同:简单的线性归回与k=1时的局部加权线性回归对比:')
print('k=1时,误差大小为:', rssError(abY[100:199], yHat1.T))
ws = standRegres(abX[0:99], abY[0:99])
yHat = np.mat(abX[100:199]) * ws
print('简单的线性回归误差大小:', rssError(abY[100:199], yHat.T.A))
可以看到，当k=0.1时，训练集误差小，但是应用于新的数据集之后，误差反而变大了。这就是经常说道的过拟合现象。我们训练的模型，我们要保证测试集准确率高，这样训练出的模型才可以应用于新的数据，也就是要加强模型的普适性。可以看到，当k=1时，局部加权线性回归和简单的线性回归得到的效果差不多。这也表明一点，必须在未知数据上比较效果才能选取到最佳模型。那么最佳的核大小是10吗？或许是，但如果想得到更好的效果，应该用10个不同的样本集做10次测试来比较结果。
本示例展示了如何使用局部加权线性回归来构建模型，可以得到比普通线性回归更好的效果。局部加权线性回归的问题在于，每次必须在整个数据集上运行。也就是说为了做出预测，必须保存所有的训练数据。
四、缩减系数来"理解"数据
如果数据样本的特征维度大于样本的数量，此时我们还能采取上面的线性回归方法求出最佳拟合参数么？显然不可能，因为当样本特征维度大于样本数时，数据矩阵显然是非满秩矩阵，那么对非满秩矩阵求逆运算会出现错误。
为了解决这个问题，科学家提出了岭回归(ridge regression)的概念，
缩减方法：
- 岭回归
- lasso法 效果很好，但是计算复杂。
此外还有一种称为"前向逐步回归"的算法，该算法可以取得很好的效果且计算相对容易。
4.1、 岭回归
简单而言，岭回归即是在矩阵xTx上加入一个λI从而使得矩阵非奇异，进而能对矩阵xTx+λI求逆。其中矩阵I是一个单位矩阵，即对角线上元素皆为1，其他均为0。这样，回归系数的计算公式变为：
公式中通过引入该惩罚项，从而减少不重要的参数，更好的理解和利用数据。此外，增加了相关约束：Σwi2<=λ，即回归系数向量中每个参数的平方和不能大于λ，这就避免了当两个或多个特征相关时，可能出现很大的正系数和很大的负系数。
上面的岭回归就是一种缩减方法，通过此方法可以去掉不重要的参数，更好的理解数据的重要性和非重要性，从而更好的预测数据。
在岭回归算法中，通过预测误差最小化来得到最优的λ值。数据获取之后，将数据随机分成两部分，一部分用于训练模型，另一部分则用于测试预测误差。为了找到最优的λ，可以选择多个不同λ值重复上述测试过程，最终得到一个使预测误差最小的λ。
岭的由来
4.2、python 实现
4.2.1、岭回归
若n阶方阵A的行列式不为零，即 |A|≠0，则称A为非奇异矩阵或满秩矩阵，否则称A为奇异矩阵或降秩矩阵。
# 岭回归
def ridgeRegres(xMat, yMat, lam=0.2):
xTx = xMat.T * xMat
# 添加 λI 使得矩阵非奇异 I是单位矩阵 即惩罚项
denom = xTx + eye(shape(xMat)[1]) * lam
if linalg.det(denom) == 0.0:
print("This matrix is singular, cannot do inverse")
return
# 岭回归的回归系数
ws = denom.I * (xMat.T * yMat)
return ws
4.2.2、特征需要标准化处理，使所有特征具有相同重要性(不考虑特征代表什么)
具体做法就是所有特征都减去各自的均值并除以方差
所有的回归系数输出到一个矩阵， 返回
# 在一组λ上测试结果
def ridgeTest(xArr, yArr):
xMat = mat(xArr); # 将列表转为矩阵
yMat = mat(yArr).T # y 需要的是一个列矩阵,mat(yArr)得到的是一个行矩阵， 所以需要转置
# 计算均值
yMean = mean(yMat, 0) # 压缩行，对各列求平均值
yMat = yMat - yMean # to eliminate X0 take mean off of Y
# regularize X's
xMeans = mean(xMat, 0) # 均值
xVar = var(xMat, 0) # 方差
np.seterr(divide='ignore', invalid='ignore') # xVar中存在0元素
# 特征标准化： (特征-均值)/方差
xMat = (xMat - xMeans) / xVar
numTestPts = 30
wMat = zeros((numTestPts, shape(xMat)[1]))
for i in range(numTestPts):
# λ以指数级变化 = exp(i - 10)
ws = ridgeRegres(xMat, yMat, exp(i - 10))
wMat[i, :] = ws.T
# 所有的回归系数输出到一个矩阵， 返回
return wMat
4.2.3、测试鲍鱼年龄
abX, abY = loadDataSet("abalone.txt")
ridgeWeights = ridgeTest(abX, abY)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(ridgeWeights)
plt.show()
需要说明的几点如下：
（1）代码中用到NumPy库中的eye()方法来生成单位矩阵。
（2）代码中仍保留了判断行列式是否为0的代码，原因是当λ取值为0时，又回到了普通的线性回归，那么矩阵很可能出现不可逆的情况
（3）岭回归中数据需要进行标准化处理，即数据的每一维度特征减去相应的特征均值之后，再除以特征的方差。
（4）这里，选择了30个不同的λ进行测试，并且这里的λ是按照指数级进行变化，从而可以看出当λ非常小和非常大的情况下对结果造成的影响
下图示出了**回归系数与log(λ)**之间的关系:
可以看出，当λ非常小时，系数与普通回归一样。而λ非常大时，所有回归系数缩减为0。这样，可以在中间的某处找到使得预测结果最好的λ。
4.3、lasso
五、前向逐步回归
逐步前向回归是一种贪心算法，即每一步都尽可能的减小误差。从一开始，所有的权重都设为1，然后每一步所做的决策是对某个权重增加或减少一个较小的数值。算法的伪代码为：
数据标准化，使其分布满足均值为0,和方差为1
在每轮的迭代中：
设置当前最小的误差为正无穷
对每个特征：
增大或减小：
改变一个系数得到一个新的w
计算新w下的误差
如果误差小于当前最小的误差：设置最小误差等于当前误差
将当前的w设置为最优的w
将本次迭代得到的预测误差最小的w存入矩阵中
返回多次迭代下的回归系数组成的矩阵
5.1、前向逐步回归代码实现
# 前向逐步回归
def stageWise(xArr,yArr,eps=0.01,numIt=100):
'''
:param xArr: x数据集
:param yArr: y数据集
:param eps: 每次迭代需要调整的步长
:param numIt: 迭代次数
:return:
'''
xMat = mat(xArr); yMat=mat(yArr).T
yMean = mean(yMat,0)
yMat = yMat - yMean #can also regularize ys but will get smaller coef
# 特征标准化
xMat = regularize(xMat)
m,n=shape(xMat)
# 将每次迭代中得到的回归系数存入矩阵
returnMat = zeros((numIt,n)) #testing code remove
ws = zeros((n,1)); # 初始化所有权重都是1，
wsTest = ws.copy(); wsMax = ws.copy()
for i in range(numIt): # 迭代次数
# 初始化最小误差为正无穷
lowestError = inf;
for j in range(n): # 遍历每个特征
for sign in [-1,1]:# 对每个特征的系数增加和减少eps*sign操作
wsTest = ws.copy()
# 改变系数
wsTest[j] += eps*sign
yTest = xMat*wsTest # 预测值
# 新的误差
rssE = rssError(yMat.A,yTest.A)
# 与所有误差比较后，取得最小误差
if rssE < lowestError:
lowestError = rssE
wsMax = wsTest
ws = wsMax.copy()
returnMat[i,:]=ws.T
return returnMat
5.1.1、使用逐步线性回归得到的系数与迭代次数关系
六、权衡偏差和方差
七、预测乐高玩具套装的价格
我们知道乐高玩具是一种拼装类玩具，由很多大小不同的塑料插件组成。一种乐高玩具套装基本上在几年后就会停产，但乐高收藏者之间仍会在停产后彼此交易。这样，我们可以拟合一个回归模型，从而对乐高套装进行估价。显然这样做十分有意义。
算法流程：
- 1 收集数据：用google shopping的api收集数据
- 2 准备数据：从返回的json数据中抽取价格
- 3 分析数据：可视化并观察数据
- 4 训练算法：构建不同的模型，采用岭回归和普通线性回归训练模型
- 5 测试算法：使用交叉验证来测试不同的模型，选择效果最好的模型
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
