5.4多元线性回归实践
5.4.1多元线性回归之房屋价格预测
在本文开始前,笔者给出的一个房价预测是一元线性的房价预测,接下来笔者要讲的也是房价预测,不在只是和面积有关,还和房间数有关,就变成了一个多元的问题。关于这部分的理论在5.2节笔者已经讲过了,这里就不在讲了。直接来代码吧。
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import axes3d
import math
from matplotlib import cm
import matplotlib.ticker as mtick
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
"""
函数说明:读取数据
Parameters:
filename - 文件名
Returns:
xArr - x数据集
yArr - y数据集
"""
def loadDataSet(filename):
'''加载数据'''
X = []
Y = []
with open(filename, 'rb') as f:
for idx, line in enumerate(f):
line = line.decode('utf-8').strip()
if not line:
continue
eles = line.split()
if idx == 0:
numFea = len(eles)
eles = list(map(float, eles))
X.append(eles[:-1])
Y.append([eles[-1]])
# 将X,Y列表转化成矩阵
xArr = np.array(X)
yArr = np.array(Y)
return xArr,yArr
"""
函数说明:特征标准化处理
Parameters:
X - 样本集
Returns:
X, values - 标准后的样本集
"""
def standarize(X):
m, n = X.shape
values = {} # 保存每一列的mean和std,便于对预测数据进行标准化
for j in range(n):
features = X[:,j]
meanVal = features.mean(axis=0)
stdVal = features.std(axis=0)
values[j] = [meanVal, stdVal]
if stdVal != 0:
X[:,j] = (features - meanVal) / stdVal
else:
X[:,j] = 0
return X, values
"""
函数说明:代价函数
Parameters:
theta, X, Y
Returns:
"""
def J(theta, X, Y):
m = len(X)
return np.sum(np.dot((predict(theta, X) - Y).T, (predict(theta, X) - Y)) / (2 * m))
"""
函数说明:梯度下降求解参数
Parameters:
alpha, X, Y, maxloop, epsilon
Returns:
theta, costs, thetas
"""
def bgd(alpha, X, Y, maxloop, epsilon):
m, n = X.shape
theta = np.zeros((n,1)) #初始化参数为0
count = 0 # 记录迭代次数
converged = False # 是否已收敛的标志
cost = np.inf # 初始化代价值为无穷大
costs = [J(theta, X, Y),] # 记录每一次的代价值
thetas = {} # 记录每一次参数的更新
for i in range(n):
thetas[i] = [theta[i,0],]
while count <= maxloop:
if converged:
break
count += 1
# n个参数计算,并存入thetas中(单独计算)
#for j in range(n):
# deriv = np.sum(np.dot(X[:,j].T, (h(theta, X) - Y))) / m
# thetas[j].append(theta[j,0] - alpha*deriv)
# n个参数在当前theta中更新
#for j in range(n):
# theta[j,0] = thetas[j][-1]
#n个参数同时计算更新值
theta = theta - alpha * 1.0 / m * np.dot(X.T, (predict(theta, X) - Y))
#添加到thetas中
for j in range(n):
thetas[j].append(theta[j,0])
# 记录当前参数的函数代价,并存入costs
cost = J(theta, X, Y)
costs.append(cost)
if abs(costs[-1] - costs[-2]) < epsilon:
converged = True
return theta, thetas, costs
"""
函数说明:定义模型函数
Parameters:
theta, X
Returns:
返回预测后的结果
"""
def predict(theta, X):
return np.dot(X, theta) # 此时的X为处理后的X
"""
函数说明:绘制回归曲线和数据点
Parameters
X, Y, theta
Returns:
无
"""
def plotRegression(X, Y, theta):
# 打印拟合平面
fittingFig = plt.figure(figsize=(16, 12))
##title = 'bgd: rate=%.3f, maxloop=%d, epsilon=%.3f \n'%(alpha,maxloop,epsilon)
ax=fittingFig.gca(projection='3d')
xx = np.linspace(0,200,25)
yy = np.linspace(0,5,25)
zz = np.zeros((25,25))
for i in range(25):
for j in range(25):
normalizedSize = (xx[i]-values[0][0])/values[0][1]
normalizedBr = (yy[j]-values[1][0])/values[1][1]
x = np.matrix([[1,normalizedSize, normalizedBr]])
zz[i,j] =predict(theta, x)
xx, yy = np.meshgrid(xx,yy)
ax.zaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
ax.plot_surface(xx, yy, zz, rstride=1, cstride=1, cmap=cm.rainbow, alpha=0.1, antialiased=True)
xs = X[:, 0].flatten()
ys = X[:, 1].flatten()
zs = Y[:, 0].flatten()
ax.scatter(xs, ys, zs, c='b', marker='o')
ax.set_xlabel(u'面积')
ax.set_ylabel(u'卧室数')
ax.set_zlabel(u'估价')
"""
函数说明:绘制代价函数曲线
Parameters
costs
Returns:
无
"""
def plotinline(costs):
errorsFig = plt.figure()
ax = errorsFig.add_subplot(111)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
ax.plot(range(len(costs)), costs)
ax.set_xlabel(u'迭代次数')
ax.set_ylabel(u'代价函数')
"""
函数说明:计算相关系数
Parameters
X, Y
Returns:
相关系数
"""
def computeCorrelation(X, Y):
xBar = np.mean(X)
yBar = np.mean(Y)
SSR = 0
varX = 0
varY = 0
for i in range(0 , len(X)):
diffXXBar = X[i] - xBar
diffYYBar = Y[i] - yBar
SSR += (diffXXBar * diffYYBar)
varX += diffXXBar**2
varY += diffYYBar**2
SST = math.sqrt(varX * varY)
return SSR / SST
#测试
if __name__ == '__main__':
## Step 1: load data...
print("Step 1: load data...")
originX, Y = loadDataSet('houses.txt')
#print(originX.shape, Y.shape)
# 对特征X增加x0列
m,n = originX.shape
X, values = standarize(originX.copy())
X = np.concatenate((np.ones((m,1)), X), axis=1)
#print(X.shape, Y.shape)
#print(X[:3],values)
## Step 2: training...
print("Step 2: training...")
alpha = 1 # 学习率
maxloop = 5000 # 最大迭代次数
epsilon = 0.000001
# 最优参数保存在theta中,costs保存每次迭代的代价值,thetas保存每次迭代更新的theta值
theta , thetas, costs= bgd(alpha, X, Y, maxloop, epsilon)
#print(theta , thetas, costs)
print(theta)
## Step 4: testing...
print("Step 4: testing...")
normalizedSize = (70-values[0][0])/values[0][1]
normalizedBr = (2-values[1][0])/values[1][1]
predicateX = np.matrix([[1, normalizedSize, normalizedBr]])
# 到此,参数学习出来了,模型也就定下来了,若要预测新的实例,进行以下即可
price = predict(theta, predicateX)
## Step 5: show the result...
print("Step 5: show the result...")
print('70㎡两居估价: ¥%.4f万元'%price)
## Step 6: show Regression plot...
print("Step 6: show Regression plot...")
plotRegression(X, Y, theta)
## Step 7: show plotinline...
print("Step 7: show plotinline...")
plotinline(costs)
## Step 8: math Pearson...
print("Step 8: math Pearson...")
xMat = np.mat(X) #创建xMat矩阵
yMat = np.mat(Y) #创建yMat矩阵
yHat = predict(theta, xMat)
Pearson = computeCorrelation(yHat, yMat)
print(Pearson)
结果如下所示。
【注】代码和一元回归的实例稍微有些不同,请读者朋友对比查看。
【完整代码参考附件3.MLR\MLR_v1\MLR_v1.0.py】
调用sklearn库实现多元线性回归
直接看代码咯。
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import math
from matplotlib import cm
import matplotlib.ticker as mtick
from sklearn import linear_model
from mpl_toolkits.mplot3d import axes3d
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
"""
函数说明:读取数据
Parameters:
filename - 文件名
Returns:
xArr - x数据集
yArr - y数据集
"""
def loadDataSet(filename):
'''加载数据'''
X = []
Y = []
with open(filename, 'rb') as f:
for idx, line in enumerate(f):
line = line.decode('utf-8').strip()
if not line:
continue
eles = line.split()
if idx == 0:
numFea = len(eles)
eles = list(map(float, eles))
X.append(eles[:-1])
Y.append([eles[-1]])
# 将X,Y列表转化成矩阵
xArr = np.array(X)
yArr = np.array(Y)
return xArr,yArr
"""
函数说明:特征标准化处理
Parameters:
X - 样本集
Returns:
X, values - 标准后的样本集
"""
def standarize(X):
m, n = X.shape
values = {} # 保存每一列的mean和std,便于对预测数据进行标准化
for j in range(n):
features = X[:,j]
meanVal = features.mean(axis=0)
stdVal = features.std(axis=0)
values[j] = [meanVal, stdVal]
if stdVal != 0:
X[:,j] = (features - meanVal) / stdVal
else:
X[:,j] = 0
return X, values
"""
函数说明:代价函数
Parameters:
theta, X, Y
Returns:
"""
def J(theta, X, Y):
m = len(X)
return np.sum(np.dot((predict(theta, X) - Y).T, (predict(theta, X) - Y)) / (2 * m))
"""
函数说明:梯度下降求解参数
Parameters:
alpha, X, Y, maxloop, epsilon
Returns:
theta, costs, thetas
"""
def bgd(alpha, X, Y, maxloop, epsilon):
m, n = X.shape
theta = np.zeros((n,1)) #初始化参数为0
count = 0 # 记录迭代次数
converged = False # 是否已收敛的标志
cost = np.inf # 初始化代价值为无穷大
costs = [J(theta, X, Y),] # 记录每一次的代价值
thetas = {} # 记录每一次参数的更新
for i in range(n):
thetas[i] = [theta[i,0],]
while count <= maxloop:
if converged:
break
count += 1
# n个参数计算,并存入thetas中(单独计算)
#for j in range(n):
# deriv = np.sum(np.dot(X[:,j].T, (h(theta, X) - Y))) / m
# thetas[j].append(theta[j,0] - alpha*deriv)
# n个参数在当前theta中更新
#for j in range(n):
# theta[j,0] = thetas[j][-1]
#n个参数同时计算更新值
theta = theta - alpha * 1.0 / m * np.dot(X.T, (predict(theta, X) - Y))
#添加到thetas中
for j in range(n):
thetas[j].append(theta[j,0])
# 记录当前参数的函数代价,并存入costs
cost = J(theta, X, Y)
costs.append(cost)
if abs(costs[-1] - costs[-2]) < epsilon:
converged = True
return theta, thetas, costs
"""
函数说明:定义模型函数
Parameters:
theta, X
Returns:
返回预测后的结果
"""
def predict(theta, X):
return np.dot(X, theta) # 此时的X为处理后的X
"""
函数说明:绘制回归曲线和数据点
Parameters
X, Y, theta
Returns:
无
"""
def plotRegression(X, Y, theta):
# 打印拟合平面
fittingFig = plt.figure(figsize=(16, 12))
##title = 'bgd: rate=%.3f, maxloop=%d, epsilon=%.3f \n'%(alpha,maxloop,epsilon)
ax=fittingFig.gca(projection='3d')
xx = np.linspace(0,200,25)
yy = np.linspace(0,5,25)
zz = np.zeros((25,25))
xx, yy = np.meshgrid(xx,yy)
ax.zaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
ax.plot_surface(xx, yy, zz, rstride=1, cstride=1, cmap=cm.rainbow, alpha=0.1, antialiased=True)
xs = X[:, 0].flatten()
ys = X[:, 1].flatten()
zs = Y[:, 0].flatten()
ax.scatter(xs, ys, zs, c='b', marker='o')
ax.set_xlabel(u'面积')
ax.set_ylabel(u'卧室数')
ax.set_zlabel(u'估价')
"""
函数说明:绘制代价函数曲线
Parameters
costs
Returns:
无
"""
def plotinline(costs):
errorsFig = plt.figure()
ax = errorsFig.add_subplot(111)
ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.2e'))
ax.plot(range(len(costs)), costs)
ax.set_xlabel(u'迭代次数')
ax.set_ylabel(u'代价函数')
"""
函数说明:计算相关系数
Parameters
X, Y
Returns:
相关系数
"""
def computeCorrelation(X, Y):
xBar = np.mean(X)
yBar = np.mean(Y)
SSR = 0
varX = 0
varY = 0
for i in range(0 , len(X)):
diffXXBar = X[i] - xBar
diffYYBar = Y[i] - yBar
SSR += (diffXXBar * diffYYBar)
varX += diffXXBar**2
varY += diffYYBar**2
SST = math.sqrt(varX * varY)
return SSR / SST
#测试
if __name__ == '__main__':
## Step 1: load data...
print("Step 1: load data...")
originX, Y = loadDataSet('houses.txt')
#print(originX.shape, Y.shape)
# 对特征X增加x0列
m,n = originX.shape
X = np.concatenate((np.ones((m,1)), originX), axis=1)
#print(X.shape, Y.shape)
#print(X[:3],values)
## Step 2: Create linear regression object...
print("Step 2: Create linear regression object...")
regr = linear_model.LinearRegression()
## Step 3: training...
print("Step 3: training...")
regr.fit(X, Y)
## Step 4: testing...
print("Step 4: testing...")
X_predict = [1, 70, 2]
# 到此,参数学习出来了,模型也就定下来了,若要预测新的实例,进行以下即可
Y_predict = regr.predict([X_predict])
## Step 5: show the result...
print("Step 5: show the result...")
#定义向量 w = (w_1,..., w_p) 作为 coef_ ,定义 w_0 作为 intercept_ 。
print("coef_:"+str(regr.coef_))
print("intercept_:"+str(regr.intercept_))
print("Y_predict:"+str(Y_predict))
## Step 6: show the result...
print("Step 6: show the result...")
print('70㎡两居估价: ¥%.4f万元'%Y_predict)
## Step 6: show Regression plot...
print("Step 6: show Regression plot...")
w0 = np.array(regr.intercept_)
w1 = np.array(regr.coef_.T[1:3])
print(w0)
print(w1)
#合并为一个矩阵
theta = np.vstack((w0,w1))
print(theta)
plotRegression(X, Y, theta)
## Step 7: math Pearson...
print("Step 7: math Pearson...")
xMat = np.mat(X) #创建xMat矩阵
yMat = np.mat(Y) #创建yMat矩阵
yHat = predict(theta, xMat)
Pearson = computeCorrelation(yHat, yMat)
print(Pearson)
【完整代码参考附件3.MLR\MLR-sklearn_v2\MLR-sklearn_v2.0.py】
【注】结果和上一个代码一样就不在贴出来了。
5.4.2多元线性回归之鲍鱼年龄预测
接下来,我们将回归用于真是数据。在abalone.txt文件中记录了鲍鱼(一种水生物【呵呵一笑,】)的年龄,这个数据来自UCI数据集合的数据。鲍鱼年龄可以从鲍鱼壳的层数推算得到。先看看数据集吧。
这个有点吓人,这也算是一比较好的例子吧。直接看代码吧。
# -*- coding:utf-8 -*-
import numpy as np
"""
函数说明:加载数据
Parameters:
fileName - 文件名
Returns:
xArr - x数据集
yArr - y数据集
"""
def loadDataSet(fileName):
numFeat = len(open(fileName).readline().split('\t')) - 1
xArr = []; yArr = []
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)
yArr.append(float(curLine[-1]))
return xArr, yArr
"""
函数说明:使用局部加权线性回归计算回归系数w
Parameters:
testPoint - 测试样本点
xArr - x数据集
yArr - y数据集
k - 高斯核的k,自定义参数
Returns:
ws - 回归系数
"""
def lwlr(testPoint, xArr, yArr, k = 1.0):
xMat = np.mat(xArr); yMat = np.mat(yArr).T
m = np.shape(xMat)[0]
weights = np.mat(np.eye((m))) #创建权重对角矩阵
for j in range(m): #遍历数据集计算每个样本的权重
diffMat = testPoint - xMat[j, :]
weights[j, j] = np.exp(diffMat * diffMat.T/(-2.0 * k**2))
xTx = xMat.T * (weights * xMat)
if np.linalg.det(xTx) == 0.0:
print("矩阵为奇异矩阵,不能求逆")
return
ws = xTx.I * (xMat.T * (weights * yMat)) #计算回归系数
return testPoint * ws
"""
函数说明:局部加权线性回归测试
Parameters:
testArr - 测试数据集,测试集
xArr - x数据集,训练集
yArr - y数据集,训练集
k - 高斯核的k,自定义参数
Returns:
ws - 回归系数
"""
def lwlrTest(testArr, xArr, yArr, k=1.0):
m = np.shape(testArr)[0] #计算测试数据集大小
yHat = np.zeros(m)
for i in range(m): #对每个样本点进行预测
yHat[i] = lwlr(testArr[i],xArr,yArr,k)
return yHat
"""
函数说明:计算回归系数w
Parameters:
xArr - x数据集
yArr - y数据集
Returns:
ws - 回归系数
"""
def standRegres(xArr,yArr):
xMat = np.mat(xArr); yMat = np.mat(yArr).T
xTx = xMat.T * xMat #根据文中推导的公示计算回归系数
if np.linalg.det(xTx) == 0.0:
print("矩阵为奇异矩阵,不能求逆")
return
ws = xTx.I * (xMat.T*yMat)
return ws
"""
函数说明:误差大小评价函数
Parameters:
yArr - 真实数据
yHatArr - 预测数据
Returns:
误差大小
"""
def rssError(yArr, yHatArr):
return ((yArr - yHatArr) **2).sum()
if __name__ == '__main__':
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次测试来比较结果。
本示例展示了如何使用局部加权线性回归来构建模型,可以得到比普通线性回归更好的效果。局部加权线性回归的问题在于,每次必须在整个数据集上运行。也就是说为了做出预测,必须保存所有的训练数据。
参考文献:
[1]机器学习基础,(英)罗杰斯,(英)吉罗拉 著,郭茂祖 等译
[2]机器学习,周志华著
[3]机器学习实战
本章参考代码
点击进入