之前介紹的分類的目標變量都是標稱型數據，接下來我們將介紹連續型的數據並且作出預測，本篇介紹的是線性回歸，接下來引入局部平滑技術，能夠更好地擬合數據
本篇我們主要討論欠擬合情況下的縮減的技術，探討偏差和方差的概念。
優點：結構易於理解，計算上不復雜
缺點：對非線性的數據擬合不好
適合數值型和標稱型數據
有回歸方程，求回歸方程的回歸系數的過程就是回歸，一旦有了回歸系數，再給定了輸入，做預測就非常容易。具體做法就是回歸系數乘以輸入數據，再將結果全部加到一起，就得到預測值
機器學習算法的基本任務就是預測，預測目標按照數據類型可以分為兩類：一種是標稱型數據（通常表現為類標簽），另一種是連續型數據（例如房價或者銷售量等等）。針對標稱型數據的預測就是我們常說的分類，針對數值型數據的預測就是回歸了。這里有一個特殊的算法需要注意，邏輯回歸（logistic regression）是一種用來分類的算法，那為什么又叫“回歸”呢？這是因為邏輯回歸是通過擬合曲線來進行分類的。也就是說，邏輯回歸只不過在擬合曲線的過程中采用了回歸的思想，其本質上仍然是分類算法
這個簡單的式子就叫回歸方程，其中0.7和0.19稱為回歸系數，面積和房子的朝向稱為特征。有了這些概念，我們就可以說，回歸實際上就是求回歸系數的過程。在這里我們看到，房價和面積以及房子的朝向這兩個特征呈線性關系，這種情況我們稱之為線性回歸。當然還存在非線性回歸，在這種情況下會考慮特征之間出現非線性操作的可能性（比如相乘或者相除），由於情況有點復雜，不在這篇文章的討論范圍之內。
簡便起見，我們規定代表輸入數據的矩陣為XX （維度為m*n，m為樣本數，n為特征維度），回歸系數向量為 θθ（維度為n*1）。對於給定的數據矩陣XX ，其預測結果由：Y=XθY=Xθ 這個式子給出。我們手里有一些現成的x和y作為訓練集，那么如何根據訓練集找到合適的回歸系數向量θθ是我們要考慮的首要問題，一旦找到θθ，預測問題就迎刃而解了。在實際應用中，我們通常認為能帶來最小平方誤差的θθ就是我們所要尋找的回歸系數向量。平方誤差指的是預測值與真實值的差的平方。采用平方這種形式的目的在於規避正負誤差的互相抵消。所以，我們的目標函數如下所示：
minθ∑i=0m(yi−xTiθ)2minθ∑i=0m(yi−xiTθ)2
這里的m代表訓練樣本的總數。對這個函數的求解有很多方法，由於網絡上對於詳細解法的相關資料太少，下面展示一種利用正規方程組的解法： 
(1) 
∇AtrAB=BT∇AtrAB=BT(2)
針對上式不太清楚的朋友可以看我之前寫的這篇博文：http://blog.csdn.net/qrlhl/article/details/47758509。根據以上式子，解法如下： 
令其等於0，即可得：θ=(XTX)−1XTyθ=(XTX)−1XTy 。有一些需要說明的地方：第三步是根據實數的跡和等於本身這一事實推導出的（括號中的每一項都為實數），第四步是根據式（2）推導出來的。第五步是根據式（1）推導出來的，其中的C為單位矩陣II。這樣，我們就得到了根據訓練集求得回歸系數矩陣θθ的方程。這種方法的特點是簡明易懂，不過缺點也很明顯，就是XTXXTX 這一項不一定可以順利的求逆。由於只有滿秩才可以求逆，這對數據矩陣XX提出了一定的要求。有人也許會問XTXXTX不是滿秩的情況下怎么辦？這個時候就要用到嶺回歸（ridge regression）了，這一部分留到下次再講。
貼上代碼
from numpy import * def loadDataSet(fileName): #general function to parse tab -delimited floats numFeat = len(open(fileName).readline().split('\t')) - 1 #get number of fields dataMat = []; labelMat = [] fr = open(fileName) for line in fr.readlines(): lineArr =[] curLine = line.strip().split('\t') for i in range(numFeat): lineArr.append(float(curLine[i])) dataMat.append(lineArr) labelMat.append(float(curLine[-1])) return dataMat,labelMat def standRegres(xArr,yArr): xMat = mat(xArr); yMat = mat(yArr).T xTx = xMat.T*xMat if linalg.det(xTx) == 0.0: print "This matrix is singular, cannot do inverse" return ws = xTx.I * (xMat.T*yMat) return ws def lwlr(testPoint,xArr,yArr,k=1.0): xMat = mat(xArr); yMat = mat(yArr).T m = shape(xMat)[0] weights = mat(eye((m))) for j in range(m): #next 2 lines create weights matrix diffMat = testPoint - xMat[j,:] # weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2)) xTx = xMat.T * (weights * xMat) if 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): #loops over all the data points and applies lwlr to each one m = shape(testArr)[0] yHat = zeros(m) for i in range(m): yHat[i] = lwlr(testArr[i],xArr,yArr,k) return yHat def lwlrTestPlot(xArr,yArr,k=1.0): #same thing as lwlrTest except it sorts X first yHat = zeros(shape(yArr)) #easier for plotting xCopy = mat(xArr) xCopy.sort(0) for i in range(shape(xArr)[0]): yHat[i] = lwlr(xCopy[i],xArr,yArr,k) return yHat,xCopy def rssError(yArr,yHatArr): #yArr and yHatArr both need to be arrays return ((yArr-yHatArr)**2).sum() def ridgeRegres(xMat,yMat,lam=0.2): xTx = xMat.T*xMat denom = xTx + eye(shape(xMat)[1])*lam if 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 = mat(xArr); yMat=mat(yArr).T yMean = mean(yMat,0) yMat = yMat - yMean #to eliminate X0 take mean off of Y #regularize X's xMeans = mean(xMat,0) #calc mean then subtract it off xVar = var(xMat,0) #calc variance of Xi then divide by it xMat = (xMat - xMeans)/xVar numTestPts = 30 wMat = zeros((numTestPts,shape(xMat)[1])) for i in range(numTestPts): ws = ridgeRegres(xMat,yMat,exp(i-10)) wMat[i,:]=ws.T return wMat def regularize(xMat):#regularize by columns inMat = xMat.copy() inMeans = mean(inMat,0) #calc mean then subtract it off inVar = var(inMat,0) #calc variance of Xi then divide by it inMat = (inMat - inMeans)/inVar return inMat def stageWise(xArr,yArr,eps=0.01,numIt=100): xMat = mat(xArr); yMat=mat(yArr).T yMean = mean(yMat,0) yMat = yMat - yMean #can also regularize ys but will get smaller coef xMat = regularize(xMat) m,n=shape(xMat) #returnMat = zeros((numIt,n)) #testing code remove ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy() for i in range(numIt): print ws.T lowestError = inf; for j in range(n): for sign in [-1,1]: wsTest = ws.copy() wsTest[j] += eps*sign yTest = xMat*wsTest rssE = rssError(yMat.A,yTest.A) if rssE < lowestError: lowestError = rssE wsMax = wsTest ws = wsMax.copy() #returnMat[i,:]=ws.T #return returnMat #def scrapePage(inFile,outFile,yr,numPce,origPrc): # from BeautifulSoup import BeautifulSoup # fr = open(inFile); fw=open(outFile,'a') #a is append mode writing # soup = BeautifulSoup(fr.read()) # i=1 # currentRow = soup.findAll('table', r="%d" % i) # while(len(currentRow)!=0): # title = currentRow[0].findAll('a')[1].text # lwrTitle = title.lower() # if (lwrTitle.find('new') > -1) or (lwrTitle.find('nisb') > -1): # newFlag = 1.0 # else: # newFlag = 0.0 # soldUnicde = currentRow[0].findAll('td')[3].findAll('span') # if len(soldUnicde)==0: # print "item #%d did not sell" % i # else: # soldPrice = currentRow[0].findAll('td')[4] # priceStr = soldPrice.text # priceStr = priceStr.replace('$','') #strips out $ # priceStr = priceStr.replace(',','') #strips out , # if len(soldPrice)>1: # priceStr = priceStr.replace('Free shipping', '') #strips out Free Shipping # print "%s\t%d\t%s" % (priceStr,newFlag,title) # fw.write("%d\t%d\t%d\t%f\t%s\n" % (yr,numPce,newFlag,origPrc,priceStr)) # i += 1 # currentRow = soup.findAll('table', r="%d" % i) # fw.close() from time import sleep import json import urllib2 def searchForSet(retX, retY, setNum, yr, numPce, origPrc): sleep(10) myAPIstr = 'AIzaSyD2cR2KFyx12hXu6PFU-wrWot3NXvko8vY' searchURL = 'https://www.googleapis.com/shopping/search/v1/public/products?key=%s&country=US&q=lego+%d&alt=json' % (myAPIstr, setNum) pg = urllib2.urlopen(searchURL) retDict = json.loads(pg.read()) for i in range(len(retDict['items'])): try: currItem = retDict['items'][i] if currItem['product']['condition'] == 'new': newFlag = 1 else: newFlag = 0 listOfInv = currItem['product']['inventories'] for item in listOfInv: sellingPrice = item['price'] if sellingPrice > origPrc * 0.5: print "%d\t%d\t%d\t%f\t%f" % (yr,numPce,newFlag,origPrc, sellingPrice) retX.append([yr, numPce, newFlag, origPrc]) retY.append(sellingPrice) except: print 'problem with item %d' % i def setDataCollect(retX, retY): searchForSet(retX, retY, 8288, 2006, 800, 49.99) searchForSet(retX, retY, 10030, 2002, 3096, 269.99) searchForSet(retX, retY, 10179, 2007, 5195, 499.99) searchForSet(retX, retY, 10181, 2007, 3428, 199.99) searchForSet(retX, retY, 10189, 2008, 5922, 299.99) searchForSet(retX, retY, 10196, 2009, 3263, 249.99) def crossValidation(xArr,yArr,numVal=10): m = len(yArr) indexList = range(m) errorMat = zeros((numVal,30))#create error mat 30columns numVal rows for i in range(numVal): trainX=[]; trainY=[] testX = []; testY = [] random.shuffle(indexList) for j in range(m):#create training set based on first 90% of values in indexList if j < m*0.9: trainX.append(xArr[indexList[j]]) trainY.append(yArr[indexList[j]]) else: testX.append(xArr[indexList[j]]) testY.append(yArr[indexList[j]]) wMat = ridgeTest(trainX,trainY) #get 30 weight vectors from ridge for k in range(30):#loop over all of the ridge estimates matTestX = mat(testX); matTrainX=mat(trainX) meanTrain = mean(matTrainX,0) varTrain = var(matTrainX,0) matTestX = (matTestX-meanTrain)/varTrain #regularize test with training params yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)#test ridge results and store errorMat[i,k]=rssError(yEst.T.A,array(testY)) #print errorMat[i,k] meanErrors = mean(errorMat,0)#calc avg performance of the different ridge weight vectors minMean = float(min(meanErrors)) bestWeights = wMat[nonzero(meanErrors==minMean)] #can unregularize to get model #when we regularized we wrote Xreg = (x-meanX)/var(x) #we can now write in terms of x not Xreg: x*w/var(x) - meanX/var(x) +meanY xMat = mat(xArr); yMat=mat(yArr).T meanX = mean(xMat,0); varX = var(xMat,0) unReg = bestWeights/varX print "the best model from Ridge Regression is:\n",unReg print "with constant term: ",-1*sum(multiply(meanX,unReg)) + mean(yMat)