1 from numpy import *
  2 import matplotlib.pyplot as plt
  3 
  4 # 计算普通线性回归，回归系数最优解
  5 def standRegres(xArr,yArr):
  6     xMat=mat(xArr)
  7     yMat=mat(yArr).T
  8     xTx=xMat.T*xMat
  9     if linalg.det(xTx)==0.0:
 10         print("this matrix is singular,cannot do inverse")
 11         return
 12     ws=xTx.I*(xMat.T*yMat)
 13     return ws
 14 
 15 # 岭回归
 16 # 计算回归系数
 17 def ridgeRegres(xMat,yMat,lam=0.2):
 18     xTx=xMat.T*xMat
 19     denom=xTx+eye(shape(xMat)[1])*lam
 20     if linalg.det(denom)==0.0:
 21         print("this matrix is singular,cannot do inverse")
 22         return
 23     ws=denom.I*(xMat.T*yMat)
 24     return ws
 25 # 在一组lambda上测试效果
 26 def ridgeTest(xArr,yArr):
 27     xMat=mat(xArr)
 28     yMat=mat(yArr).T
 29     yMean=mean(yMat,0)
 30     yMat=yMat-yMean
 31     xMeans=mean(xMat,0)
 32     xVar=var(xMat,0)
 33     xMat=(xMat-xMeans)/xVar
 34     numTestPts=30
 35     wMat=zeros((numTestPts,shape(xMat)[1]))
 36     for i in range(numTestPts):
 37         ws=ridgeRegres(xMat,yMat,exp(i-10))
 38         wMat[i,:]=ws.T
 39     return wMat
 40 
 41 # 2.交叉验证岭回归
 42 def crossValidation(xArr,yArr,numVal=10):
 43     m=len(yArr)
 44     indexList=list(range(m))  # 创建一组包含0-76的顺序列表
 45     errorMat=zeros((numVal,30))  # 用于存放错误率的矩阵
 46     # 10折交叉验证
 47     for i in range(numVal):
 48         trainX=[]
 49         trainY=[]
 50         testX=[]
 51         testY=[]
 52         # 打乱列表，以使每一折的训练集和测试集都不一样
 53         random.shuffle(indexList)
 54         # 分配训练集和测试集，前90%放入训练集，后10%放入测试集
 55         for j in range(m):
 56             if j<m*0.9:
 57                 trainX.append(xArr[indexList[j]])
 58                 trainY.append(yArr[indexList[j]])
 59             else:
 60                 testX.append(xArr[indexList[j]])
 61                 testY.append(yArr[indexList[j]])
 62             # 使用岭回归得到回归系数矩阵，用了30个lambda创建了三十组回归系数
 63         wMat=ridgeTest(trainX,trainY)
 64         for k in range(30):
 65             matTestX=mat(testX)
 66             matTrainX=mat(trainX)
 67             meanTrain=mean(matTrainX,0)
 68             varTrain=var(matTrainX,0)
 69             # 测试集使用与训练集相同的参数进行标准化
 70             matTestX=(matTestX-meanTrain)/varTrain
 71             # 预测测试集的价格
 72             yEst=matTestX*mat(wMat[k,:]).T+mean(trainY)
 73             # 计算预测偏差
 74             errorMat[i,k]=rssError(yEst.T.A,array(testY))
 75     # 对各列求平均值
 76     meanErrors=mean(errorMat,0)
 77     minMean = float(min(meanErrors))
 78     # nonzero 返回不为0元素下标
 79     # 平均误差最小的组的回归系数即为所求最佳
 80     bestWeights = wMat[nonzero(meanErrors == minMean)]
 81     xMat = mat(xArr);
 82     yMat = mat(yArr).T
 83     meanX = mean(xMat, 0);
 84     varX = var(xMat, 0)
 85     # 数据还原标准化以前
 86     unReg = bestWeights / varX
 87     constant = -1 * sum(multiply(meanX, unReg)) + mean(yMat)  # 常数项
 88     print("the best model from Ridge Regression is:\n", unReg)
 89     print("with constant term: ", constant)
 90     return unReg, constant
 91 
 92 # 1.乐高玩具，数据依次是年份，零部件数，新旧，原始价格
 93 xArr,yArr = loadDataSet("legoAllData.txt")
 94 xMat=mat(xArr)
 95 lgX=xMat
 96 yMat=mat(yArr)
 97 lgX1=mat(ones((76,5)))
 98 lgX1[:,1:5]=xMat
 99 lgY=yMat
100 
101 # 2.普通线性回归
102 ws=standRegres(lgX1,lgY)
103 print(ws)
104 # 售价约等于=59349-29*年份-.0.03*零部件数+55*新旧+2*原始价格
105 
106 # 3.交叉验证岭回归
107 crossValidation(xArr,yArr,10)