你對樂高（LEGO）品牌的玩具瞭解嗎？樂高公司生產拼裝類玩具，由很多大小不同的塑料插塊組成。這些塑料插塊的設計非常出色，不需要任何粘合劑就可以隨意拼裝起來。除了簡單玩具之外，樂高玩具在一些成人中也很流行。一般來說，這些插塊都成套出售，它們可以拼裝成很多不同的東西，如船、城堡、一些著名建築，等等。樂高公司每個套裝包含的部件數目從10件到5000件不等。
一種樂高套裝基本上在幾年後就會停產，但樂高的收藏者之間仍會在停產後彼此交易。Dangler喜歡為樂高套裝估價，下面將用本章的回歸技術幫助他建立一個預測模型。
示例：用回歸法預測樂高套裝的價格
- 收集數據：用Google Shopping的API收集數據。
- 準備數據：從返回的JSON數據中抽取價格。
- 分析數據：可視化並觀察數據。
- 訓練算法：構建不同的模型，採用逐步線性回歸和直接的線性回歸模型。
- 測試算法：使用交叉驗證來測試不同的模型，分析哪個效果最好。
- 使用算法：這次練習的目標就是生成數據模型。
在這個例子中，我們將從不同的數據集上獲取價格，然後對這些數據建立回歸模型。需要做的第一件事就是如何獲取數據。
8.6.1 收集數據：使用Google購物的API
Google已經為我們提供了一套購物的API來抓取價格。在使用API之前，需要註冊一個Google賬號，然後訪問Google API的控制台來確保購物API可用。完成之後就可以發送HTTP請求，API將以JSON格式返回所需的產品信息。Python提供了JSON解析模塊，我們可以從返回的JSON格式裡整理出所需數據。詳細的API介紹可以參見：http://code.google.com/apis/shopping/search/v1/getting_started.html。
打開regression.py文件並加入如下代碼。
程序清單8-5 購物信息的獲取函數
from time import sleep
import json
import urllib2
def searchForSet(retX, retY, setNum, yr, numPce, origPrc):
sleep(10)
myAPIstr = \'get from code.google.com\'
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 \"%dt%dt%dt%ft%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)
上述程序清單中的第一個函數是searchForSet，它調用Google購物API並保證數據抽取的正確性。這裡需要導入新的模塊：time.sleep、json和urllib2。但是一開始要休眠10秒鐘，這是為了防止短時間內有過多的API調用。接下來，我們拼接查詢的URL字符串，添加API的key和待查詢的套裝信息，打開和解析操作通過json.loads方法實現。完成後我們將得到一部字典，下面需要做的是從中找出價格和其他信息。
部分返回結果的是一個產品的數組，我們將在這些產品上循環迭代，判斷該產品是否是新產品並抽取它的價格。我們知道，樂高套裝由很多小插件組成，有的二手套裝很可能會缺失其中一兩件。也就是說，賣家只出售套裝的若幹部件（不完整）。因為這種不完整的套裝也會通過檢索結果返回，所以我們需要將這些信息過濾掉（可以統計描述中的關鍵詞或者是用貝葉斯方法來判斷）。我在這裡僅使用了一個簡單的啟髮式方法：如果一個套裝的價格比原始價格低一半以上，則認為該套裝不完整。程序清單8-5在代碼❶處過濾掉了這些套裝，解析成功後的套裝將在屏幕上顯示出來並保存在list對像retX和retY中。
程序清單8-5的最後一個函數是 setDataCollect，它負責多次調用searchForSet。函數searchForSet的其他參數是從www.brickset.com收集來的，它們也一併輸出到文件中。
下面看一下執行結果，添加程序清單8-5中的代碼之後保存regression.py，在Python提示符下輸入如下命令：
>>> lgX = ; lgY =
>>> regression.setDataCollect(lgX, lgY)
2006 800 1 49.990000 549.990000
2006 800 1 49.990000 759.050000
2006 800 1 49.990000 316.990000
2002 3096 1 269.990000 499.990000
2002 3096 1 269.990000 289.990000
.
.
2009 3263 0 249.990000 524.990000
2009 3263 1 249.990000 672.000000
2009 3263 1 249.990000 580.000000
檢查一下lgX和lgY以確認一下它們非空。下節我們將使用這些數據來構建回歸方程並預測樂高玩具套裝的售價。
8.6.2 訓練算法：建立模型
上一節從網上收集到了一些真實的數據，下面將為這些數據構建一個模型。構建出的模型可以對售價做出預測，並幫助我們理解現有數據。看一下Python是如何完成這些工作的。
首先需要添加對應常數項的特徵X0（X0=1），為此創建一個全1的矩陣：
>>> shape(lgX) (58, 4) >>> lgX1=mat(ones((58,5)))
接下來，將原數據矩陣lgX複製到新數據矩陣lgX1的第1到第5列：
>>> lgX1[:,1:5]=mat(lgX)
確認一下數據複製的正確性：
>>> lgX[0]
[2006.0, 800.0, 0.0, 49.990000000000002]
>>> lgX1[0]
matrix([[ 1.00000000e+00, 2.00600000e+03, 8.00000000e+02,
0.00000000e+00, 4.99900000e+01]])
很顯然，後者除了在第0列加入1之外其他數據都一樣。最後在這個新數據集上進行回歸處理：
>>> ws=regression.standRegres(lgX1,lgY)
>>> ws
matrix([[ 5.53199701e+04],
[ -2.75928219e+01],
[ -2.68392234e-02],
[ -1.12208481e+01],
[ 2.57604055e+00]])
檢查一下結果，看看模型是否有效：
>>> lgX1[0]*ws
matrix([[ 76.07418853]])
>>> lgX1[-1]*ws
matrix([[ 431.17797672]])
>>> lgX1[43]*ws
matrix([[ 516.20733105]])
可以看到模型有效。下面看看具體的模型。該模型認為套裝的售價應該採用如下公式計算：
$55319.97-27.59*Year-0.00268*NumPieces-11.22*NewOrUsed+2.57*original price
這個模型的預測效果非常好，但模型本身並不能令人滿意。它對於數據擬合得很好，但看上去沒有什麼道理。從公式看，套裝裡零部件越多售價反而會越低。另外，該公式對新套裝也有一定的懲罰。
下面使用縮減法中一種，即嶺回歸再進行一次實驗。前面討論過如何對係數進行縮減，但這次將會看到如何用縮減法確定最佳回歸係數。打開regression.py並輸入下面的代碼。
程序清單8-6 交叉驗證測試嶺回歸
def crossValidation(xArr,yArr,numVal=10):
m = len(yArr)
indexList = range(m)
errorMat = zeros((numVal,30))
for i in range(numVal):
#❶（以下兩行）創建訓練集和測試集容器
trainX=; trainY=
testX = ; testY =
random.shuffle(indexList)
for j in range(m):
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)
for k in range(30):
#❸（以下三行）用訓練時的參數將測試數據標準化
matTestX = mat(testX); matTrainX=mat(trainX)
meanTrain = mean(matTrainX,0)
varTrain = var(matTrainX,0)
matTestX = (matTestX-meanTrain)/varTrain
yEst = matTestX * mat(wMat[k,:]).T + mean(trainY)
errorMat[i,k]=rssError(yEst.T.A,array(testY))
meanErrors = mean(errorMat,0)
minMean = float(min(meanErrors))
bestWeights = wMat[nonzero(meanErrors==minMean)]
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)
上述程序清單中的函數crossValidation有三個參數，前兩個參數lgX和lgY存有數據集中的X和Y值的list對象，默認lgX和lgY具有相同的長度。第三個參數numVal是算法中交叉驗證的次數，如果該值沒有指定，就取默認值10。函數crossValidation首先計算數據點的個數 m。創建好了訓練集和測試集容器❶，之後創建了一個list並使用Numpy提供的random.shuffle函數對其中的元素進行混洗（shuffle），因此可以實現訓練集或測試集數據點的隨機選取。❷處將數據集的90%分割成訓練集，其餘10%為測試集，並將二者分別放入對應容器中。
一旦對數據點進行混洗之後，就建立一個新的矩陣wMat來保存嶺回歸中的所有回歸係數。我們還記得在8.4.1節中，函數ridgeTest使用30個不同的λ值創建了30組不同的回歸係數。接下來我們也在上述測試集上用30組回歸係數來循環測試回歸效果。嶺回歸需要使用標準化後的數據，因此測試數據也需要用與測試集相同的參數來執行標準化。在❸處用函數rssError計算誤差，並將結果保存在errorMat中。
在所有交叉驗證完成後，errorMat保存了ridgeTest裡每個λ對應的多個誤差值。為了將得出的回歸係數與standRegres作對比，需要計算這些誤差估計值的均值1。有一點值得注意：嶺回歸使用了數據標準化，而standRegres則沒有，因此為了將上述比較可視化還需將數據還原。在❹處對數據做了還原並將最終結果展示。
1. 此處為10折，所以每個λ應該對應10個誤差。應選取使誤差的均值最低的回歸係數。----譯者注
來看一下整體的運行效果，在regression.py中輸入程序清單8-6中的代碼並保存，然後執行如下命令：
>>> regression.crossValidation(lgX,lgY,10)
The best model from Ridge Regression is:
[[ -2.96472902e+01 -1.34476433e-03 -3.38454756e+01 2.44420117e+00]]
with constant term: 59389.2069537
為了便於與常規的最小二乘法進行比較，下面給出當前的價格公式：
$59389.21-29.64*Year-0.00134*NumPieces-33.85*NewOrUsed+2.44*original price.
可以看到，該結果與最小二乘法沒有太大差異。我們本期望找到一個更易於理解的模型，顯然沒有達到預期效果。為了達到這一點，我們來看一下在縮減過程中回歸係數是如何變化的，輸入下面的命令：
>>> regression.ridgeTest(lgX,lgY)
array([[ -1.45288906e+02, -8.39360442e+03, -3.28682450e+00, 4.42362406e+04],
[ -1.46649725e+02, -1.89952152e+03, -2.80638599e+00,4.27891633e+04],
.
.
[ -4.91045279e-06, 5.01149871e-08, 2.40728171e-05,8.14042912e-07]])
這些係數是經過不同程度的縮減得到的。首先看第1行，第4項比第2項的係數大5倍，比第1項大57倍。這樣看來，如果只能選擇一個特徵來做預測的話，我們應該選擇第4個特徵，也就是原始價格。如果可以選擇2個特徵的話，應該選擇第4個和第2個特徵。
這種分析方法使得我們可以挖掘大量數據的內在規律。在僅有4個特徵時，該方法的效果也許並不明顯；但如果有100個以上的特徵，該方法就會變得十分有效：它可以指出哪些特徵是關鍵的，而哪些特徵是不重要的。