线性回归与岭回归python代码实现

发布时间:2019-08-09 10:34:14编辑:auto阅读(1589)

    一、标准线性回归

     

    在线性回归中我们要求的参数为:

    详细的推导可以参见:http://blog.csdn.net/weiyongle1996/article/details/73727505

    所以代码实现主要就是实现上式,python代码如下:

    import numpy as np
    import pandas as pd
    import matplotlib.pyplot as plt
    
    # implement stand regress
    def standRegress(xArr,yArr):
        # 将数组转换为矩阵
        xMat = np.mat(xArr)
        yMat = np.mat(yArr)
        xTx = xMat.T * xMat # 计算xTx的
        if np.linalg.det(xTx) == 0.0:
            print('xTx不能求逆矩阵')
            return
        theta = xTx.I * (xMat.T * yMat)
        yHat = xMat*theta
        return yHat
    
    # import data
    ex0 = np.loadtxt('ex0.txt',delimiter='\t')
    
    # deal with data
    xArr = []
    yArr = []
    for data in ex0:
        # print(data)
        xTmp = []
        yTmp = []
        xTmp.append(data[0])
        xTmp.append(data[1])
        yTmp.append(data[2])
        xArr.append(xTmp)
        yArr.append(yTmp)
    
    # print(ex0)
    # print(xArr[0:2])
    print(yArr)
    # ws = standRegress(xArr,yArr)
    # print(ws)
    yHat = standRegress(xArr,yArr)
    xMat = np.mat(xArr)
    yMat = np.mat(yArr)
    # print(yMat.T[0,:].flatten().A[0])
    plt.scatter(xMat[:,1].flatten().A[0],yMat.T[0,:].flatten().A[0]) # real data
    plt.plot(xMat[:,1],yHat,'r-') # predict data
    plt.show()
    

    运行结果如下:


    二、局部加权线性回归


    局部加权线性回归是在线性回归的基础上增加权值,以更好的拟合弯曲的线段(详细参见:http://blog.csdn.net/weiyongle1996/article/details/74537567

    使用局部加权解出的回归系数为:

    python代码如下:

    import numpy as np
    import matplotlib.pyplot as plt
    
    def lwlr(testPoint,xArr,yArr,k=1.0):
        xMat = np.mat(xArr)
        yMat = np.mat(yArr)
        m = np.shape(xMat)[0] #shape 读取矩阵的长度 shape[0]获得矩阵第一维的长度
        # print(m)
        weights = np.mat(np.eye(m)) # 创建对角矩阵
        # print(weights)
        for j in range(m):                      #next 2 lines create weights matrix
            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("This matrix is singular, cannot do inverse")
            return
        ws = xTx.I * (xMat.T * (weights * yMat))
        return testPoint * 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
    
    # import data
    ex0 = np.loadtxt('ex0.txt',delimiter='\t')
    
    # deal with data
    xArr = []
    yArr = []
    for data in ex0:
        # print(data)
        xTmp = []
        yTmp = []
        xTmp.append(data[0])
        xTmp.append(data[1])
        yTmp.append(data[2])
        xArr.append(xTmp)
        yArr.append(yTmp)
    
    # 对单点估计
    # yHat = lwlr(xArr[0],xArr,yArr,1.0)
    # print(yHat)
    # 得到所有点的估计
    yHat = lwlrTest(xArr,xArr,yArr,0.02)
    xMat = np.mat(xArr)
    yMat = np.mat(yArr)
    # print(xMat)
    strInd = xMat[:,1].argsort(0) # argsort返回数组值从小到大排列后各元素对应的索引值
    # print(strInd)
    xSort = xMat[strInd][:,0,:] # 排序
    # print(xSort)
    
    plt.scatter(xMat[:,1].flatten().A[0],yMat.T[0,:].flatten().A[0]) # real data
    plt.plot(xSort[:,1],yHat[strInd],'r-') # predict data
    plt.show()
    
    运行结果如下:


    更改k的值会获得不同的曲线,k越小,对真实数据拟合的越好(但可能过拟合),k越大,越趋向于标准的线性回归。


    三、岭回归


    岭回归就是在矩阵xTx上增加一项使得矩阵非奇异,从而能够对其求逆。从上面两端代码我们可以看到,在之前对xTx求逆时都需要先判断xTx是否可以求逆,而岭回归就是解决这个问题的。岭回归的回归系数计算公式为:


    实现代码如下:

    import numpy as np
    import matplotlib.pyplot as plt
    
    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:
            print("This matrix is singular, cannot do inverse")
            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) # 数据标准化
        # print(yMean)
        yMat = yMat - yMean
        # print(xMat)
        #regularize X's
        xMeans = np.mean(xMat,0)
        xVar = np.var(xMat,0)
        xMat = (xMat - xMeans) / xVar #(特征-均值)/方差
        numTestPts = 30
        wMat = np.zeros((numTestPts,np.shape(xMat)[1]))
        for i in range(numTestPts): # 测试不同的lambda取值,获得系数
            ws = ridgeRegres(xMat,yMat,np.exp(i-10))
            wMat[i,:]=ws.T
        return wMat
    
    
    # import data
    ex0 = np.loadtxt('abalone.txt',delimiter='\t')
    xArr = ex0[:,0:-1]
    yArr = ex0[:,-1]
    # print(xArr,yArr)
    ridgeWeights = ridgeTest(xArr,yArr)
    # print(ridgeWeights)
    plt.plot(ridgeWeights)
    plt.show()
    运行结果如图:

    纵坐标为回归系数,横坐标为log(lambda),在最左边,回归系数与线性回归一致,最右边系数全部缩减为0.

    其中间某部分可以得到最好的预测结果,为了定量进行寻找最佳参数,还需要进行交叉验证。


    以上代码python环境均为python3.6

    代码参考:

    《机器学习实战》

    数据取自《机器学习实战》附带数据

关键字