基于pyhton3.6-机器学习实战-回归regression代码解释



本人是一名数学系研究生，于2017年底第一次接触python和机器学习，作为一名新手，欢迎与大家交流。

我主要给大家讲解代码，理论部分给大家推荐3本书：

《机器学习实战中文版》

《机器学习》周志华

《统计学习方法》李航

以上3本书，第一本是基于python2的代码实现；剩余两本主要作为第一本书理论省略部分的补充，理论大部分都讲得很细。

博客上关于机器学习实战理论解释都很多，参差不齐，好作品也大都借鉴了以上3本书，网上有很多电子版的书。

与其看看一些没用的博客，真心不如以上3本书有收获。

说实话，学习一定要静下心来，切忌浮躁。不懂可以每天看一点，每天你懂一点，天天积累就多了。

操作系统：windows8.1

python版本：python3.6

运行环境：spyder（anaconda）

# -*- coding: utf-8 -*-
"""
Created on Fri Mar  9 16:38:23 2018

@author: Loulch C.C
"""

import matplotlib.pyplot as plt
import numpy as np

def loadDataSet(fileName):
    """
    函数说明:数据导入函数
    Parameters:
        fileName - 文件名
    Returns:
        xArr - 数据矩阵
        yArr - 目标值，类似于分类中的类别标签
    """
    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

def standRegres(xArr,yArr):
    """
    函数说明:标准回归函数
    Parameters:
        xArr - 数据矩阵
        yArr - 目标值，类似于分类中的类别标签
    Returns:
        ws - 回归系数
    """
    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)
    print('标准回归函数计算的回归系数是：',ws)
    return ws

def plotRegression():
    """
    函数说明:绘制回归曲线和数据点
    Parameters:
        无
    Returns:
        无
    """
    xArr, yArr = loadDataSet('ex0.txt')                                    
    ws = standRegres(xArr, yArr)                                     
    xMat = np.mat(xArr)                                                   
    yMat = np.mat(yArr)                                                 
    xCopy = xMat.copy()  
    #print('xCopy:',xCopy)                                                  
    xCopy.sort(0) 
    #print('xCopy.sort(0):',xCopy)                                                       
    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.ylabel('Y')
    plt.show()
"""
if __name__ == '__main__':
    plotRegression()
#"""

#
def lwlr(testPoint, xArr, yArr, k = 1.0):
    """
    函数说明:局部加权线性回归函数（高斯核）
    Parameters:
        testPoint - 测试样本点
        xArr - 数据矩阵
        yArr - 目标值，类似于分类中的类别标签
        k - 高斯核的k,自定义参数
    Returns:
        ws - 回归系数
    """
    xMat = np.mat(xArr); yMat = np.mat(yArr).T
    m = np.shape(xMat)[0]
    weights = np.mat(np.eye((m)))                  #创建m*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)                   #m*m矩阵                     
    if np.linalg.det(xTx) == 0.0:
        print("矩阵为奇异矩阵,不能求逆")
        return
    ws = xTx.I * (xMat.T * (weights * yMat))          #计算回归系数，m*1矩阵
    return testPoint * ws
def lwlrTest(testArr, xArr, yArr, k=1.0):  
    """
    函数说明:局部加权线性回归测试函数
    Parameters:
        testArr - 测试数据矩阵
        xArr - 训练数据矩阵
        yArr - 目标值，类似于分类中的类别标签
        k - 高斯核的k,自定义参数
    Returns:
        ws - 回归系数
    """
    m = np.shape(testArr)[0]                                            #计算测试数据集大小
    yHat = np.zeros(m)    
    for i in range(m):                                                    #对每个样本点进行预测
        yHat[i] = lwlr(testArr[i],xArr,yArr,k)
    return yHat

from matplotlib.font_manager import FontProperties
#font_manager是matplotlib中的字体管理器

def plotlwlrRegression():
    """
    函数说明:绘制多条局部加权回归曲线
    Parameters:
        无
    Returns:
        无
    """
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)  
    #设置字体,simsun是宋体，TTC字体是TrueType字体集成文件
    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矩阵
    #print('xMat=',xMat)
    yMat = np.mat(yArr)                                                   #创建yMat矩阵
    srtInd = xMat[:, 1].argsort(0)        
    #首先xMat[:, 1]将xMat的第一列进行复制，返回一个行数*1的矩阵
    #对于矩阵来说，argsort(0)返回的是：矩阵对列排序后，值从小到大的索引值，返回一个行数*1的矩阵
    #print('xMat[:, 1]=',xMat[:, 1] )
    #print('srtInd',srtInd)
    #print(type(srtInd))
    xSort = xMat[srtInd][:,0,:]                           #对xMat按索引排序，并转换成正常矩阵
    #print('xMat[srtInd]',xMat[srtInd])                    #xMat[srtInd]只有两列       
    #print('xSort',xSort)
    #print('xSort的类型：',type(xSort))
    fig, axs = plt.subplots(nrows=3, ncols=1,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')                        #绘制回归曲线
    #print('xMat[:,1]=',xMat[:,1])
#    xMat[:,1]= [[ 0.067732]
#                [ 0.42781 ]
#                [ 0.995731]
#                ..., 
#                [ 0.070095]
#                [ 0.52707 ]
#                [ 0.116163]]
    #print('xMat[:,1].flatten()=',xMat[:,1].flatten())
#    xMat[:,1].flatten()= [[ 0.067732  0.42781   0.995731 ...,  0.070095  0.52707   0.116163]]
    #print('xMat[:,1].flatten().A=',xMat[:,1].flatten().A)
#    xMat[:,1].flatten().A= [[ 0.067732  0.42781   0.995731 ...,  0.070095  0.52707   0.116163]]
    #print('xMat[:,1].flatten().A[0]=',xMat[:,1].flatten().A[0])
#    xMat[:,1].flatten().A[0]= [ 0.067732  0.42781   0.995731 ...,  0.070095  0.52707   0.116163]
    axs[0].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 2, c = 'blue', alpha = .5)
    axs[1].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 2, c = 'blue', alpha = .5)
    axs[2].scatter(xMat[:,1].flatten().A[0], yMat.flatten().A[0], s = 2, c = 'blue', alpha = .5)
    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')  #bold粗体
    plt.setp(axs1_title_text, size=8, weight='bold', color='red')  
    plt.setp(axs2_title_text, size=8, weight='bold', color='red')  
    plt.show()
"""
if __name__ == '__main__':
    plotlwlrRegression()
#"""
#示例：预测鲍鱼年龄
def rssError(yArr, yHatArr):
    """
    平方误差函数
    Parameters:
        yArr - 真实数据
        yHatArr - 预测数据
    Returns:
        误差大小
    """
    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('更换数据集,测试结果如下:')
    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))
    #可以看到，k=10时测试误差最小，但是在训练集上的误差是最大的
    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))
#"""
#岭回归
def ridgeRegres(xMat,yMat,lam=0.2):
    xTx = xMat.T*xMat
    denom = xTx + np.eye(np.shape(xMat)[1])*lam
    if np.linalg.det(denom) == 0.0:    #若lam=0，矩阵行列式任可能为0，所以还得检查
        print( "该矩阵式奇异矩阵，不能求逆")
        return
    ws = denom.I * (xMat.T*yMat)
    return ws
    
def ridgeTest(xArr,yArr):
    xMat = np.mat(xArr); yMat=np.mat(yArr).T
    yMean = np.mean(yMat,0)       #求取y的均值
    """
    numpy.mean(a, axis=None)
    axis 不设置值，对 m*n 个数求均值，返回一个实数
    axis = 0：对各列求均值，返回 1* 列 矩阵，对于矩阵来说（对于数组，返回1*列）
    axis =1 ：对各行求均值，返回 行 *1 矩阵，对于矩阵来说（对于数组，返回1*行）
    """
    yMat = yMat - yMean           #计算y减去y的均值
    ############数据标准化：使每维特征具有相同重要性#############
    xMeans = np.mean(xMat,0)      #求取xMat每列的均值，返回一个1*列矩阵，对于矩阵来说
    xVar = np.var(xMat,0)         #求取xMat每列的方差，返回一个1*列矩阵,对于矩阵来说，方法类似于mean
    xMat = (xMat - xMeans)/xVar   #数据标准化
    ##########################################################
    numTestPts = 30
    wMat = np.zeros((numTestPts,np.shape(xMat)[1]))
    for i in range(numTestPts):   #在30个不同的lambda下调用ridgeRegres(),这里lambda以指数级变化
        ws = ridgeRegres(xMat,yMat,np.exp(i-10))         #这里lambda以指数级变化
        wMat[i,:]=ws.T
    return wMat
"""
if __name__ =='__main__':
    abX,abY=loadDataSet('abalone.txt')
    ridgeWeights=ridgeTest(abX,abY)
    #print(ridgeWeights)
    fig=plt.figure()
    ax=fig.add_subplot(111)
    ax.plot(ridgeWeights)     #ridgeWeights是一个30*8矩阵，plot按每列（特征）画图
    plt.xlabel('log(lambda)')
    plt.ylabel('regression coefficients of diffent feature')
    plt.show()
    #在最左边，即lambda最小时，所有系数的原始值与线性回归一致；而在最右边，系数全部缩减为0；
    #在中间部分的某值将可以取得最好的预测效果。为了定量地找到最佳参数，还需要进行交叉验证。
#"""
# 交叉验证
def crossValidation(xArr, yArr, numVal=10):
    m = len(yArr)
    indexList = list(range(m))
    errorMat = np.zeros((numVal, 30))  # 初始化一个10*30的错误数组
    for i in range(numVal):            #进行10次交叉验证
        trainX = []
        trainY = []
        testX = []
        testY = []
        ############随机取90%的数据作为训练集，剩下10%的数据作为测试集###############
        np.random.shuffle(indexList)    
        #对于一维数组，将indexList中数随即打乱，np.random.shuffle() 直接对原始矩阵进行修改（返回值为NoneType）
        for j in range(m):
            if j < m * 0.9:
                trainX.append(xArr[indexList[j]])
                trainY.append(yArr[indexList[j]])
            else:
                testX.append(xArr[indexList[j]])
                testY.append(yArr[indexList[j]])
        ##########################################################################
        
        #####################计算30个岭回归lambda参数下的平方误差#####################
        wMat = ridgeTest(trainX, trainY)  
        #调用ridgeTest(),在30个不同lambda下计算训练集的岭回归系数
        for k in range(30):  
            matTestX = np.mat(testX)
            matTrainX = np.mat(trainX)
            ######################数据标准化##############
            meanTrain = np.mean(matTrainX, 0)
            varTrain = np.var(matTrainX, 0)
            matTestX = (matTestX - meanTrain) / varTrain  #用训练集来数据标准化测试集
            ##############################################
            yEst = matTestX * np.mat(wMat[k, :]).T + np.mean(trainY) #计算测试结果yEst
            errorMat[i, k] = rssError(yEst.T.A, np.array(testY))     #计算平方误差
            #print(errorMat[i, k])     
            #errorMat是储存10个测试集（10行），每列是30个不同lambd下的平方误差
        ###########################################################################
        
    ###############取多次迭代的误差平均，使用平均误差最小的岭回归参数计算出的结果################
    meanErrors = np.mean(errorMat, 0) #对errorMat每列（即30个不同lambda下的平方误差）求均值
    minMean = float(min(meanErrors))
    bestWeights = wMat[np.nonzero(meanErrors == minMean)] 
    #nonzeros(a)返回数组a中非零元素的索引值数组
    xMat = np.mat(xArr)
    yMat = np.mat(yArr).T
    meanX = np.mean(xMat, 0)
    varX = np.var(xMat, 0)
    ##############################数据还原##################################
    unReg = bestWeights / varX     #想与标准线性回归做对比，所以进行数据还原
    print("最好的岭回归模型是:\n", unReg)
    print("常数项是: ", -1 * np.sum(np.multiply(meanX, unReg)) + np.mean(yMat))  
    #y=W^T*X+C, -1 * np.sum(np.multiply(meanX, unReg)) = -1 * meanX* unReg.T
    #######################################################################
    #########################################################################################
"""
if __name__ =='__main__':
    abX,abY=loadDataSet('abalone.txt')
    crossValidation(abX,abY)
    standRegres(abX,abY)
#"""

def regularize(xMat):
    xMeans = np.mean(xMat,0)      #求取xMat每列的均值，返回一个1*列矩阵，对于矩阵来说
    xVar = np.var(xMat,0)         #求取xMat每列的方差，返回一个1*列矩阵,对于矩阵来说，方法类似于mean
    xMat = (xMat - xMeans)/xVar   #数据标准化
    return xMat


"""
前向逐步线性回归（一种贪婪算法）:
数据标准化
每轮迭代： 
    设置当前最小误差lowestError为正无穷 
    对每个特征：
        增大或者缩小：（两次循环）
            改变一个系数得到新的w 
            计算新的w下的误差rss 
            如果新的误差rss小于当前最小误差：则设置wbest=当前的w 
        将当前w设置为最新的wbest
""" 

def stageWise(xArr, yArr, eps = 0.01, numIt = 100):  #eps步长;numIt迭代次数
    xMat = np.mat(xArr)  
    yMat =np. mat(yArr).T  
    yMean = np.mean(yMat, 0)  
    yMat = yMat - yMean  
    xMat = regularize(xMat)  
    m,n = np.shape(xMat)  
    returnMat = np.zeros((numIt,n))         #用来储存numIt次迭代误差最小的回归系数ws
    ws = np.zeros((n,1))                    #初始化回归系数ws
    wsTest = ws.copy()                      #用来储存测试中的回归系数ws
    wsMax = ws.copy()                       #用来储存误差最小的回归系数ws
    for i in range(numIt):                  #进入迭代
        print(ws.T)  
        lowestError = float('inf')          #设置当前最小误差lowestError为正无穷
        for j in range(n):                  #遍历每个特征
            for sign in [-1,1]:             #通过1，-1的切换来实现增大或缩小
                wsTest = ws.copy()  
                wsTest[j] += eps*sign       #改变一个系数，形成一个新的ws
                yTest = xMat*wsTest         #计算y的预测值
                rssE = rssError(yMat.A,yTest.A)  #计算新ws下的平方误差
                if rssE < lowestError:  
                    lowestError = rssE  
                    wsMax = wsTest  
        ws = wsMax.copy()  
        returnMat[i,:] = ws.T  
    return returnMat  

"""
if __name__ =='__main__':
    xArr,yArr=loadDataSet('abalone.txt')
    stageWise(xArr,yArr,0.01,200)
#"""