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

py3study

发布于 2020-01-13 15:36:14

1.5K0

发布于 2020-01-13 15:36:14

文章被收录于专栏：python3python3

一、标准线性回归

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

详细的推导可以参见：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()

运行结果如图：