《人工智能》机器学习 - 第5章 线性回归(三 多元线性回归实践)

102 阅读11分钟

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]机器学习实战

本章参考代码
点击进入