Machine Learning in Action(7) 回归算法
按照《机器学习实战》的主线，结束有监督学习中关于分类的机器学习方法，进入回归部分。所谓回归就是数据进行曲线拟合，回归一般用来做预测，涵盖线性回归（经典最小二乘法）、局部加权线性回归、岭回归和逐步线性回归。先来看下线性回归，即经典最小二乘法，说到最小二乘法就不得说下线性代数，因为一般说线性回归只通过计算一个公式就可以得到答案，如（公式一）所示：
（公式一）
其中X是表示样本特征组成的矩阵，Y表示对应的值，比如房价，股票走势等，（公式一）是直接通过对（公式二）求导得到的，因为（公式二）是凸函数，导数等于零的点就是最小点。
（公式二）
不过并不是所有的码农能从（公式二）求导得到（公式一）的解，因此这里给出另外一个直观的解，直观理解建立起来后，后续几个回归就简单类推咯。从初中的投影点说起，如（图一）所示：
（图一）
在（图一）中直线a上离点b最近的点是点b在其上的投影，即垂直于a的交点p。p是b在a上的投影点。试想一下，如果我们把WX看成多维的a,即空间中的一个超面来代替二维空间中的直线，而y看成b,那现在要使得（公式二）最小是不是就是寻找（图一）中的e,即垂直于WX的垂线，因为只有垂直时e才最小。下面来看看如何通过寻找垂线并最终得到W。要寻找垂线，先从（图二）中的夹角theta 说起吧，因为当cos(theta)=0时，他们也就垂直了。下面来分析下直线或者向量之间的夹角，如（图二）所示：
（图二）
在（图二）中，表示三角形
的斜边，那么：
角beta也可以得到同样的计算公式，接着利用三角形和差公式得到（公式三）：
（公式三）
（公式三）表示的是两直线或者两向量之间的夹角公式，很多同学都学过。再仔细看下，发现分子其实是向量a,b之间的内积（点积），因此公式三变为简洁的（公式四）的样子：
（公式四）
接下来继续分析（图一）中的投影，为了方便观看，增加了一些提示如（图三）所示：
（图三）
在（图三）中，假设向量b在向量a中的投影为p（注意，这里都上升为向量空间，不再使用直线，因为（公式四）是通用的）。投影p和a 在同一方向上（也可以反方向），因此我们可以用一个系数乘上a来表示p,比如（图三）中的，有了投影向量p，那么我们就可以表示向量e,因为根据向量法则，
，有因为a和e垂直，因此
，展开求得系数x,如（公式五）所示：
（公式五）
（公式五）是不是很像（公式一）？只不过（公式一）的分母写成了另外的形式，不过别急，现在的系数只是一个标量数字，因为a,b都是一个向量，我们要扩展一下，把a从向量扩展到子空间，因为（公式一）中的X是样本矩阵，矩阵有列空间和行空间，如（图四）所示：
（图四）
（图四）中的A表示样本矩阵X，假设它有两个列a1和a2,我们要找一些线性组合系数来找一个和（图三）一样的接受b 投影的向量，而这个向量通过矩阵列和系数的线性组合表示。求解的这个系数的思路和上面完全一样，就是寻找投影所在的向量和垂线e的垂直关系，得到系数，如（公式六）所示：
（公式六）
这下（公式六）和（公式一）完全一样了，基于最小二乘法的线性回归也就推导完成了，而局部加权回归其实只是相当于对不同样本之间的关系给出了一个权重，所以叫局部加权，如（公式七）所示：
（公式七）
而权重的计算可通过高斯核（高斯公式）来完成，核的作用就是做权重衰减，很多地方都要用到，表示样本的重要程度，一般离目标进的重要程度大些，高斯核可以很好的描述这种关系。如（公式八）所示，其中K是个超参数，根据情况灵活设置：
（公式八）
（图五）是当K分别为1.0, 0.01,0.003时的局部加权线性回归的样子，可以看出当K=1.0时，和线性回归没区别：
（图五）
而岭回归的样子如（公式九）所示：
（公式九）
岭回归主要是解决的问题就是当XX’无法求逆时，比如当特征很多，样本很少，矩阵X不是满秩矩阵，此时求逆会出错，但是通过加上一个对角为常量lambda的矩阵，就可以很巧妙的避免这个计算问题，因此会多一个参数lambda,lambda的最优选择由交叉验证（cross-validation）来决定，加上一个对角不为0的矩阵很形象的在对角上抬高了，因此称为岭。不同的lambda会使得系数缩减，如（图六）所示：
（图六）
说到系数缩减大家可能会觉得有奇怪，感觉有点类似于正则，但是这里只是相当于在（公式六）中增大分母，进而缩小系数，另外还有一些系数缩减的方法，比如直接增加一些约束，如（公式十）和（公式十一）所示：
（公式十）
（公式十一）
当线性回归增加了（公式十）的约束变得和桥回归差不多，系数缩减了，而如果增加了（公式十一）的约束时就是稀疏回归咯，（我自己造的名词，sorry）,系数有一些0。
有了约束后，求解起来就不像上面那样直接计算个矩阵运算就行了，回顾第五节说中支持向量机原理，需要使用二次规划求解，不过仍然有一些像SMO算法一样的简化求解算法，比如前向逐步回归方法：
前向逐步回归的伪代码如（图七）所示，也不难，仔细阅读代码就可以理解：
（图七）
下面直接给出上面四种回归的代码：
1 from numpy import * 2 3 def loadDataSet(fileName): #general function to parse tab -delimited floats 4 numFeat = len(open(fileName).readline().split('\t')) - 1 #get number of fields 5 dataMat = []; labelMat = [] 6 fr = open(fileName) 7 for line in fr.readlines(): 8 lineArr =[] 9 curLine = line.strip().split('\t') 10 for i in range(numFeat): 11 lineArr.append(float(curLine[i])) 12 dataMat.append(lineArr) 13 labelMat.append(float(curLine[-1])) 14 return dataMat,labelMat 15 16 def standRegres(xArr,yArr): 17 xMat = mat(xArr); yMat = mat(yArr).T 18 xTx = xMat.T*xMat 19 if linalg.det(xTx) == 0.0: 20 print "This matrix is singular, cannot do inverse" 21 return 22 ws = xTx.I * (xMat.T*yMat) 23 return ws 24 25 def lwlr(testPoint,xArr,yArr,k=1.0): 26 xMat = mat(xArr); yMat = mat(yArr).T 27 m = shape(xMat)[0] 28 weights = mat(eye((m))) 29 for j in range(m): #next 2 lines create weights matrix 30 diffMat = testPoint - xMat[j,:] # 31 weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2)) 32 xTx = xMat.T * (weights * xMat) 33 if linalg.det(xTx) == 0.0: 34 print "This matrix is singular, cannot do inverse" 35 return 36 ws = xTx.I * (xMat.T * (weights * yMat)) 37 return testPoint * ws 38 39 def lwlrTest(testArr,xArr,yArr,k=1.0): #loops over all the data points and applies lwlr to each one 40 m = shape(testArr)[0] 41 yHat = zeros(m) 42 for i in range(m): 43 yHat[i] = lwlr(testArr[i],xArr,yArr,k) 44 return yHat 45 46 def lwlrTestPlot(xArr,yArr,k=1.0): #same thing as lwlrTest except it sorts X first 47 yHat = zeros(shape(yArr)) #easier for plotting 48 xCopy = mat(xArr) 49 xCopy.sort(0) 50 for i in range(shape(xArr)[0]): 51 yHat[i] = lwlr(xCopy[i],xArr,yArr,k) 52 return yHat,xCopy 53 54 def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays 55 return ((yArr-yHatArr)**2).sum() 56 57 def ridgeRegres(xMat,yMat,lam=0.2): 58 xTx = xMat.T*xMat 59 denom = xTx + eye(shape(xMat)[1])*lam 60 if linalg.det(denom) == 0.0: 61 print "This matrix is singular, cannot do inverse" 62 return 63 ws = denom.I * (xMat.T*yMat) 64 return ws 65 66 def ridgeTest(xArr,yArr): 67 xMat = mat(xArr); yMat=mat(yArr).T 68 yMean = mean(yMat,0) 69 yMat = yMat - yMean #to eliminate X0 take mean off of Y 70 #regularize X's 71 xMeans = mean(xMat,0) #calc mean then subtract it off 72 xVar = var(xMat,0) #calc variance of Xi then divide by it 73 xMat = (xMat - xMeans)/xVar 74 numTestPts = 30 75 wMat = zeros((numTestPts,shape(xMat)[1])) 76 for i in range(numTestPts): 77 ws = ridgeRegres(xMat,yMat,exp(i-10)) 78 wMat[i,:]=ws.T 79 return wMat 80 81 def regularize(xMat):#regularize by columns 82 inMat = xMat.copy() 83 inMeans = mean(inMat,0) #calc mean then subtract it off 84 inVar = var(inMat,0) #calc variance of Xi then divide by it 85 inMat = (inMat - inMeans)/inVar 86 return inMat 87 88 def stageWise(xArr,yArr,eps=0.01,numIt=100): 89 xMat = mat(xArr); yMat=mat(yArr).T 90 yMean = mean(yMat,0) 91 yMat = yMat - yMean #can also regularize ys but will get smaller coef 92 xMat = regularize(xMat) 93 m,n=shape(xMat) 94 #returnMat = zeros((numIt,n)) #testing code remove 95 ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy() 96 for i in range(numIt): 97 print ws.T 98 lowestError = inf; 99 for j in range(n): 100 for sign in [-1,1]: 101 wsTest = ws.copy() 102 wsTest[j] += eps*sign 103 yTest = xMat*wsTest 104 rssE = rssError(yMat.A,yTest.A) 105 if rssE < lowestError: 106 lowestError = rssE 107 wsMax = wsTest 108 ws = wsMax.copy() 109 #returnMat[i,:]=ws.T 110 #return returnMat 111 112 #def scrapePage(inFile,outFile,yr,numPce,origPrc): 113 # from BeautifulSoup import BeautifulSoup 114 # fr = open(inFile); fw=open(outFile,'a') #a is append mode writing 115 # soup = BeautifulSoup(fr.read()) 116 # i=1 117 # currentRow = soup.findAll('table', r="%d" % i) 118 # while(len(currentRow)!=0): 119 # title = currentRow[0].findAll('a')[1].text 120 # lwrTitle = title.lower() 121 # if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): 122 # newFlag = 1.0 123 # else: 124 # newFlag = 0.0 125 # soldUnicde = currentRow[0].findAll('td')[3].findAll('span') 126 # if len(soldUnicde)==0: 127 # print "item #%d did not sell" % i 128 # else: 129 # soldPrice = currentRow[0].findAll('td')[4] 130 # priceStr = soldPrice.text 131 # priceStr = priceStr.replace('$','') #strips out $ 132 # priceStr = priceStr.replace(',','') #strips out , 133 # if len(soldPrice)>1: 134 # priceStr = priceStr.replace('Free shipping', '') #strips out Free Shipping 135 # print "%s\t%d\t%s" % (priceStr,newFlag,title) 136 # fw.write("%d\t%d\t%d\t%f\t%s\n" % (yr,numPce,newFlag,origPrc,priceStr)) 137 # i += 1 138 # currentRow = soup.findAll('table', r="%d" % i) 139 # fw.close() 140 141 from time import sleep 142 import json 143 import urllib2 144 def searchForSet(retX, retY, setNum, yr, numPce, origPrc): 145 sleep(10) 146 myAPIstr = 'AIzaSyD2cR2KFyx12hXu6PFU-wrWot3NXvko8vY' 147 searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum) 148 pg = urllib2.urlopen(searchURL) 149 retDict = json.loads(pg.read()) 150 for i in range(len(retDict['items'])): 151 try: 152 currItem = retDict['items'][i] 153 if currItem['product']['condition'] == 'new': 154 newFlag = 1 155 else: newFlag = 0 156 listOfInv = currItem['product']['inventories'] 157 for item in listOfInv: 158 sellingPrice = item['price'] 159 if sellingPrice > origPrc * 0.5: 160 print "%d\t%d\t%d\t%f\t%f" % (yr,numPce,newFlag,origPrc, sellingPrice) 161 retX.append([yr, numPce, newFlag, origPrc]) 162 retY.append(sellingPrice) 163 except: print 'problem with item %d' % i 164 165 def setDataCollect(retX, retY): 166 searchForSet(retX, retY, 8288, 2006, 800, 49.99) 167 searchForSet(retX, retY, 10030, 2002, 3096, 269.99) 168 searchForSet(retX, retY, 10179, 2007, 5195, 499.99) 169 searchForSet(retX, retY, 10181, 2007, 3428, 199.99) 170 searchForSet(retX, retY, 10189, 2008, 5922, 299.99) 171 searchForSet(retX, retY, 10196, 2009, 3263, 249.99) 172 173 def crossValidation(xArr,yArr,numVal=10): 174 m = len(yArr) 175 indexList = range(m) 176 errorMat = zeros((numVal,30))#create error mat 30columns numVal rows 177 for i in range(numVal): 178 trainX=[]; trainY=[] 179 testX = []; testY = [] 180 random.shuffle(indexList) 181 for j in range(m):#create training set based on first 90% of values in indexList 182 if j < m*0.9: 183 trainX.append(xArr[indexList[j]]) 184 trainY.append(yArr[indexList[j]]) 185 else: 186 testX.append(xArr[indexList[j]]) 187 testY.append(yArr[indexList[j]]) 188 wMat = ridgeTest(trainX,trainY) #get 30 weight vectors from ridge 189 for k in range(30):#loop over all of the ridge estimates 190 matTestX = mat(testX); matTrainX=mat(trainX) 191 meanTrain = mean(matTrainX,0) 192 varTrain = var(matTrainX,0) 193 matTestX = (matTestX-meanTrain)/varTrain #regularize test with training params 194 yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)#test ridge results and store 195 errorMat[i,k]=rssError(yEst.T.A,array(testY)) 196 #print errorMat[i,k] 197 meanErrors = mean(errorMat,0)#calc avg performance of the different ridge weight vectors 198 minMean = float(min(meanErrors)) 199 bestWeights = wMat[nonzero(meanErrors==minMean)] 200 #can unregularize to get model 201 #when we regularized we wrote Xreg = (x-meanX)/var(x) 202 #we can now write in terms of x not Xreg: x*w/var(x) - meanX/var(x) +meanY 203 xMat = mat(xArr); yMat=mat(yArr).T 204 meanX = mean(xMat,0); varX = var(xMat,0) 205 unReg = bestWeights/varX 206 print "the best model from Ridge Regression is:\n",unReg 207 print "with constant term: ",-1*sum(multiply(meanX,unReg)) + mean(yMat)
以上各种回归方法没有考虑实际数据的噪声，如果噪声很多，直接用上述的回归不是太好，因此需要加上正则，然后迭代更新权重
以上内容转载至朋友博客:http://blog.csdn.net/marvin521/article/details/9387305
