一、标准线性回归
在线性回归中我们要求的参数为：
详细的推导可以参见：http://blog.csdn.net/weiyongle1996/article/details/73727505
所以代码实现主要就是实现上式，python代码如下：
代码语言：javascript复制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代码如下：
代码语言：javascript复制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是否可以求逆，而岭回归就是解决这个问题的。岭回归的回归系数计算公式为：
实现代码如下：
代码语言：javascript复制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
代码参考：
《机器学习实战》
数据取自《机器学习实战》附带数据
