一、回歸介紹

監(jiān)督學(xué)習(xí)指的是有目標(biāo)變量或預(yù)測目標(biāo)的機器學(xué)習(xí)方法〉斓牛回歸與分類的不同，就在于其目標(biāo)變量是連續(xù)數(shù)值型嗡午，最直接的辦法是依據(jù)輸入寫出一個目標(biāo)值的計算公式囤躁。假如你想要預(yù)測姐姐男友汽車的功率大小，可能會這么計算：

HorsePower = 0.0015 * annualSalary - 0.99 * hoursListeningToPublicRadio

這就是所謂的回歸方程（regression equation）荔睹，其中的0.0015和-0.99稱作回歸系數(shù)（regression weights）狸演，求這些回歸系數(shù)的過程就是回歸。一旦有了這些回歸系數(shù)僻他，再給定輸入宵距，做預(yù)測就非常容易了。具體的做法是用回歸系數(shù)乘以輸入值吨拗，再將結(jié)果全部加在一起满哪，就得到了預(yù)測值。

說到回歸劝篷，一般都是指線性回歸（linear regression）哨鸭，所以本章里的回歸和線性回歸代表同一個意思。線性回歸意味著可以將輸入項分別乘以一些常量娇妓，再將結(jié)果加起來得到輸出像鸡。需要說明的是，存在另一種稱為非線性回歸的回歸模型哈恰，該模型不認(rèn)同上面的做法只估，比如認(rèn)為輸出可能是輸入的乘積志群。這樣，上面的功率計算公式也可以寫做：HorsePower = 0.0015 * annualSalary / hoursListeningToPublicRadio

二蛔钙、用線性回歸找到最佳擬合直線

全文代碼在最后第八點里锌云。

應(yīng)該怎么從一大堆數(shù)據(jù)里求出回歸方程呢？假定輸入數(shù)據(jù)存放在矩陣X中吁脱，結(jié)果存放在向量y中：

而回歸系數(shù)存放在向量w中：

那么對于給定的數(shù)據(jù)x1宾抓，即矩陣X的第一列數(shù)據(jù)，預(yù)測結(jié)果u1將會通過如下公式給出：

現(xiàn)在的問題是豫喧，手里有數(shù)據(jù)矩陣X和對應(yīng)的標(biāo)簽向量y石洗，怎么才能找到w呢？一個常用的方法就是找出使誤差最小的w紧显。這里的誤差是指預(yù)測u值和真實y值之間的差值讲衫，使用該誤差的簡單累加將使得正差值和負(fù)差值相互抵消，所以我們采用平方誤差孵班。

平方誤差和可以寫做：

用矩陣表示還可以寫做：

找到w涉兽，使平方誤差和最小。因為我們認(rèn)為平方誤差和越小篙程，說明線性回歸擬合效果越好枷畏。

現(xiàn)在，我們用矩陣表示的平方誤差和對w進行求導(dǎo)：

令上述公式等于0虱饿，得到：

w上方的小標(biāo)記表示拥诡，這是當(dāng)前可以估計出的w的最優(yōu)解。從現(xiàn)有數(shù)據(jù)上估計出的w可能并不是數(shù)據(jù)中的真實w值氮发，所以這里使用了一個"帽"符號來表示它僅是w的一個最佳估計渴肉。

值得注意的是，上述公式中包含逆矩陣爽冕，也就是說仇祭，這個方程只在逆矩陣存在的時候使用，也即是這個矩陣是一個方陣颈畸，并且其行列式不為0乌奇。

上述的最佳w求解是統(tǒng)計學(xué)中的常見問題，除了矩陣方法外還有很多其他方法可以解決眯娱。通過調(diào)用NumPy庫里的矩陣方法礁苗，我們可以僅使用幾行代碼就完成所需功能。該方法也稱作OLS困乒， 意思是“普通小二乘法”（ordinary least squares）寂屏。

現(xiàn)有一數(shù)據(jù)集：

第一列都為1.0贰谣，即x0娜搂。第二列為x1迁霎，即x軸數(shù)據(jù)。第三列為x2百宇，即y軸數(shù)據(jù)考廉。首先繪制下數(shù)據(jù)，看下數(shù)據(jù)分布携御。

如何判斷擬合曲線的擬合效果的如何呢昌粤？我們可以使用corrcoef方法計算這兩個序列的相關(guān)系數(shù)來比較預(yù)測值和真實值的相關(guān)性。

可以看到啄刹，對角線上的數(shù)據(jù)是1.0涮坐，因為yMat和自己的匹配是完美的，而YHat和yMat的相關(guān)系數(shù)為0.98誓军。

三袱讹、局部加權(quán)線性回歸

線性回歸的一個問題是有可能出現(xiàn)欠擬合現(xiàn)象，因為它求的是具有最小均方誤差的無偏估計昵时。顯而易見捷雕，如果模型欠擬合將不能取得最好的預(yù)測效果。所以有些方法允許在估計中引入一些偏差壹甥，從而降低預(yù)測的均方誤差救巷。

中的一個方法是局部加權(quán)線性回歸（Locally Weighted Linear Regression，LWLR）句柠。在該方法中浦译，我們給待預(yù)測點附近的每個點賦予一定的權(quán)重。與kNN一樣溯职，這種算法每次預(yù)測均需要事先選取出對應(yīng)的數(shù)據(jù)子集管怠。該算法解除回歸系數(shù)W的形式如下：

其中 w 是一個矩陣，用來給每個數(shù)據(jù)點賦予權(quán)重缸榄。LWLR使用“核”（與支持向量機中的核類似）來對附近的點賦予更高的權(quán)重渤弛。核的類型可以自由選擇，最常用的核就是高斯核甚带，高斯核對應(yīng)的權(quán)重如下：

這樣就構(gòu)建了一個只含對角元素的權(quán)重矩陣 w 她肯，并且點 x 與 x(i) 越近， w(i,i) 將會越大鹰贵。上述公式包含一個需要用戶指定的參數(shù) k 晴氨，它決定了對附近的點賦予多大的權(quán)重，這也是使用LWLR時唯一需要考慮的參數(shù)碉输。

當(dāng)k = 1.0時權(quán)重很大籽前，如同將所有的數(shù)據(jù)視為等權(quán)重，得出的最佳擬合直線與標(biāo)準(zhǔn)的回歸一致。使用k = 0.01得到了非常好的效果枝哄，抓住了數(shù)據(jù)的潛在模式肄梨。下圖使用k = 0.003納入了太多的噪聲點，擬合的直線與數(shù)據(jù)點過于貼近挠锥。所以众羡，圖中的最下圖是過擬合的一個例子，而最上圖則是欠擬合的一個例子蓖租。

四粱侣、示例：預(yù)測鮑魚的年齡

abalone.txt文件中記錄了鮑魚的年齡，這個數(shù)據(jù)來自UCI數(shù)據(jù)集合的數(shù)據(jù)蓖宦。鮑魚年齡可以從鮑魚殼的層數(shù)推算得到齐婴。

最后一列代表的是鮑魚的真實年齡，前面幾列的數(shù)據(jù)是一些鮑魚的特征稠茂，例如鮑魚殼的層數(shù)等尔店。

可以看到，當(dāng)k=0.1時主慰，訓(xùn)練集誤差小嚣州，但是應(yīng)用于新的數(shù)據(jù)集之后，誤差反而變大了共螺。這就是經(jīng)常說道的過擬合現(xiàn)象该肴。我們訓(xùn)練的模型，我們要保證測試集準(zhǔn)確率高藐不，這樣訓(xùn)練出的模型才可以應(yīng)用于新的數(shù)據(jù)匀哄，也就是要加強模型的普適性〕可以看到涎嚼，當(dāng)k=1時，局部加權(quán)線性回歸和簡單的線性回歸得到的效果差不多挑秉。這也表明一點法梯，必須在未知數(shù)據(jù)上比較效果才能選取到最佳模型。那么最佳的核大小是10嗎犀概？或許是立哑，但如果想得到更好的效果，應(yīng)該用10個不同的樣本集做10次測試來比較結(jié)果姻灶。

四铛绰、縮減系數(shù)來“理解”數(shù)據(jù)

如果數(shù)據(jù)的特征比樣本點還多，就不可以使用線性回歸和之前的方法來做預(yù)測了产喉，這是因為在計算 (X^T X)^ -1的時候會出錯捂掰。如果特征比樣本點還多（n > m）敢会，也就是說輸入數(shù)據(jù)的矩陣X不是滿秩矩陣。非滿秩矩陣在求逆時會出現(xiàn)問題这嚣。

為了解決這個問題鸥昏，統(tǒng)計學(xué)家引入了嶺回歸（ridge regression）、lasso法和前向逐步回歸疤苹。

4.1嶺回歸

嶺回歸即我們所說的L2正則線性回歸，在一般的線性回歸最小化均方誤差的基礎(chǔ)上增加了一個參數(shù)w的L2范數(shù)的罰項敛腌，從而最小化罰項殘差平方和：

簡單說來卧土，嶺回歸就是在矩陣X^TX上加一個λI 從而使得矩陣非奇異，進而能對X^TX + λI求逆像樊。其中矩陣 I 是一個m×m的單位矩陣尤莺，對角線上元素全為1，其他元素全為0生棍。而λ是一個用戶定義的數(shù)值颤霎，后面會做介紹。在這種情況下涂滴，回歸系數(shù)的計算公式將變成：

嶺回歸最先用來處理特征數(shù)多于樣本數(shù)的情況友酱，現(xiàn)在也用于在估計中加入偏差，從而得到更好的估計柔纵。這里通過引入λ來限制了所有w之和缔杉，通過引入該懲罰項，能夠減少不重要的參數(shù)搁料，這個技術(shù)在統(tǒng)計學(xué)中也可以叫做縮減（shrinkage）或详。

上圖繪制了回歸系數(shù)與log(λ)的關(guān)系。在最左邊郭计，即λ最小時霸琴，可以得到所有系數(shù)的原始值（與線性回歸一致）；而在右邊昭伸，系數(shù)全部縮減成0梧乘；在中間部分的某個位置，將會得到最好的預(yù)測結(jié)果庐杨。想要得到最佳的λ參數(shù)宋下，可以使用交叉驗證的方式獲得。

4.2前向逐步線性回歸

前向逐步線性回歸算法屬于一種貪心算法辑莫，即每一步都盡可能減少誤差学歧。我們計算回歸系數(shù)，不再是通過公式計算各吨，而是通過每次微調(diào)各個回歸系數(shù)枝笨，然后計算預(yù)測誤差袁铐。那個使誤差最小的一組回歸系數(shù)馋吗，就是我們需要的最佳回歸系數(shù)蛛勉。

上圖是反應(yīng)了迭代次數(shù)與回歸系數(shù)的關(guān)系曲線⊥兀可以看到徙融，有些系數(shù)從始至終都是約為0的洒缀，這說明它們不對目標(biāo)造成任何影響，也就是說這些特征很可能是不需要的欺冀。逐步線性回歸算法的優(yōu)點在于它可以幫助人們理解有的模型并做出改進树绩。當(dāng)構(gòu)建了一個模型后，可以運行該算法找出重要的特征隐轩，這樣就有可能及時停止對那些不重要特征的收集饺饭。

五、示例职车、預(yù)測樂高玩具套件的價格

樂高（LEGO）公司生產(chǎn)拼裝類玩具瘫俊，由很多大小不同的塑料插塊組成。一般來說悴灵，這些插塊都是成套出售扛芽，它們可以拼裝成很多不同的東西，如船积瞒、城堡胸哥、一些著名建筑等。樂高公司每個套裝包含的部件數(shù)目從10件到5000件不等赡鲜。

一種樂高套件基本上在幾年后就會停產(chǎn)空厌，但樂高的收藏者之間仍會在停產(chǎn)后彼此交易。本次實例银酬，就是使用回歸方法對收藏者之間的交易價格進行預(yù)測嘲更。

5.1、獲取數(shù)據(jù)
全部的數(shù)據(jù)集在各個網(wǎng)頁當(dāng)中揩瞪，需要利用爬蟲只是解析抓取出來赋朦。

抓取數(shù)據(jù)結(jié)果如下：

5.2、建立模型
抓取下來數(shù)據(jù)集后李破，接下來就是訓(xùn)練模型宠哄。首先需要添加全為0的特征X0列。因為線性回歸的第一列特征要求都是1.0嗤攻。然后使用最簡單的普通線性回歸進行建模毛嫉。訓(xùn)練出來的線性回歸模型如截圖最下面所示：

我們使用嶺回歸，通過交叉驗證妇菱，找到使誤差最小的λ對應(yīng)的回歸系數(shù)承粤。

這里隨機選取樣本暴区，因為其隨機性，所以每次運行的結(jié)果可能略有不同辛臊。不過整體如上圖所示仙粱，可以看出，它與常規(guī)的最小二乘法彻舰，即普通的線性回歸沒有太大差異伐割。我們本期望找到一個更易于理解的模型，顯然沒有達到預(yù)期效果刃唤。

現(xiàn)在隔心，我們看一下在縮減過程中回歸系數(shù)是如何變化的。

看運行結(jié)果的第一行透揣，可以看到最大的是第4項济炎，第二大的是第2項川抡。

因此辐真，如果只選擇一個特征來做預(yù)測的話，我們應(yīng)該選擇第4個特征崖堤，也就是原始加個侍咱。如果可以選擇2個特征的話，應(yīng)該選擇第4個和第2個特征密幔。

這種分析方法使得我們可以挖掘大量數(shù)據(jù)的內(nèi)在規(guī)律楔脯。在僅有4個特征時，該方法的效果也許并不明顯胯甩；但如果有100個以上的特征昧廷，該方法就會變得十分有效：它可以指出哪個特征是關(guān)鍵的，而哪些特征是不重要的偎箫。

六木柬、應(yīng)用scikit-learn構(gòu)建線性回歸模型

sklearn.linear_model提供了很多線性模型，包括嶺回歸淹办、貝葉斯回歸眉枕、Lasso等。現(xiàn)在我們用sklearn實現(xiàn)嶺回歸怜森。

我正則化項系數(shù)設(shè)為0.5速挑，其余參數(shù)使用默認(rèn)「惫瑁可以看到姥宝，獲得的結(jié)果與上面的結(jié)果類似。

七恐疲、小結(jié)：

與分類一樣伶授，回歸也是預(yù)測目標(biāo)值的過程断序。回歸與分類的不同點在于糜烹，前者預(yù)測連續(xù)類型變量违诗，而后者預(yù)測離散類型變量。

在回歸方程里疮蹦，求得特征對應(yīng)的最佳回歸系數(shù)的方法是最小化誤差的平方和诸迟。給定輸入矩陣X，如果X^TX的逆存在并可以求得的話愕乎，回歸法都可以直接使用阵苇。數(shù)據(jù)集上計算出的回歸方程并不一定意味著它是最佳的，可以使用預(yù)測值y_hat和原始值y的相關(guān)性來度量回歸方程的好壞感论。

當(dāng)數(shù)據(jù)的樣本數(shù)比特征數(shù)還少時候绅项，矩陣 X^TX的逆不能直接計算。即便當(dāng)樣本數(shù)比特征數(shù)多時比肄， X^TX的逆仍有可能無法直接計算快耿，這是因為特征有可能高度相關(guān)。這時可以考慮使用嶺回歸芳绩，因為當(dāng) X^TX的逆不能計算時掀亥，它仍保證能求得回歸參數(shù)。

嶺回歸是縮減法的一種妥色，相當(dāng)于對回歸系數(shù)的大小施加了限制搪花。另一種很好的縮減法是lasso。Lasso難以求解嘹害，但可以使用計算簡便的逐步線性回歸方法來求得近似結(jié)果撮竿。

縮減法還可以看做是對一個模型增加偏差的同時減少方差。

八笔呀、貼上全部代碼：

# -*- conding: utf-8 -*-
import matplotlib.pyplot as plt
from numpy import *
from matplotlib.font_manager import FontProperties
from bs4 import BeautifulSoup


def load_data_set(file_name):
    """
    Function:
        加載數(shù)據(jù)集
    Parameters:
        file_name - 文件名
    Returns:
        data_mat - 數(shù)據(jù)列表
        label_mat - 標(biāo)簽列表
    Modify:
        2018-10-22
    """
    num_feat = len(open(file_name).readline().split('\t')) - 1
    data_mat = []
    label_mat = []
    fr = open(file_name)
    for line in fr.readlines():
        line_arr = []
        cur_line = line.strip().split('\t')
        for i in range(num_feat):
            line_arr.append(float(cur_line[i]))
        data_mat.append(line_arr)
        label_mat.append(float(cur_line[-1]))
    return data_mat, label_mat


def plot_data_set():
    """
    Function:
        繪制數(shù)據(jù)集
    Parameters:
        無
    Returns:

    Modify:
        2018-10-22
    """
    data_mat, label_mat = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
    n = len(data_mat)
    xcord = []
    ycord = []
    for i in range(n):
        xcord.append(data_mat[i][1])
        ycord.append(label_mat[i])
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(xcord, ycord, s=20, c='blue', alpha=0.5)
    plt.title('data_set')
    plt.xlabel('X')
    plt.show()


def stand_regres(x_arr, y_arr):
    """
    Function:
        計算回歸系數(shù)w
    Parameters:
        x_arr - 數(shù)據(jù)列表
        y_arr - 標(biāo)簽列表
    Returns:
        ws - 回歸系數(shù)
    Modify:
        2018-10-22
    """
    x_mat = mat(x_arr)
    y_mat = mat(y_arr).T
    xTx = x_mat.T * x_mat
    if linalg.det(xTx) == 0.0:
        print('This matrix is singular, cannot do inverse')
        return
    ws = xTx.I * (x_mat.T * y_mat)
    return ws


def plot_regression():
    """
    Function:
        加載數(shù)據(jù)集
    Parameters:
        無
    Returns:
        dataMat - 數(shù)據(jù)列表
        labelMat - 標(biāo)簽列表
    Modify:
        2018-10-22
    """
    data_mat, label_mat = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
    ws = stand_regres(data_mat, label_mat)
    x_mat = mat(data_mat)
    y_mat = mat(label_mat)
    x_coyp = x_mat.copy()
    x_coyp.sort(0)
    y_hat = x_coyp * ws
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(x_coyp[:, 1], y_hat, c='red')
    ax.scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
    plt.title('data_set')
    plt.xlabel('X')
    plt.show()


def lwlr(test_point, x_arr, y_arr, k=1.0):
    """
    Function:
        使用局部加權(quán)線性回歸計算回歸系數(shù)w
    Parameters:
        test_point -  測試樣本點
        x_arr - x數(shù)據(jù)集
        y_arr - y數(shù)據(jù)集
    Returns:
        ws - 回歸系數(shù)
    Modify:
        2018-10-24
    """
    x_mat = mat(x_arr)
    y_mat = mat(y_arr).T
    m = shape(x_mat)[0]
    # 生成對角矩陣
    weights = mat(eye((m)))
    for j in range(m):
        diff_mat = test_point - x_mat[j, :]
        weights[j, j] = exp(diff_mat * diff_mat.T / (-2.0 * k ** 2))
    x_t_x = x_mat.T * (weights * x_mat)
    if linalg.det(x_t_x) == 0.0:
        print('矩陣為奇異矩陣,不能求逆')
        return
    # 計算回歸系數(shù)
    ws = x_t_x.I * (x_mat.T * (weights * y_mat))
    return test_point * ws


def lwlr_test(test_arr, x_arr, y_arr, k=1.0):
    """
    Function:
        局部加權(quán)線性回歸測試
    Parameters:
        test_arr -  測試數(shù)據(jù)集
        x_arr - x數(shù)據(jù)集
        y_arr - y數(shù)據(jù)集
    Returns:
        dataMat - 數(shù)據(jù)列表
        labelMat - 標(biāo)簽列表
    Modify:
        2018-10-24
    """
    m = shape(test_arr)[0]
    y_hat = zeros(m)
    for i in range(m):
        y_hat[i] = lwlr(test_arr[i], x_arr, y_arr, k)
    return y_hat


def plot_lwlr_regression():
    """
    Function:
        繪制多條局部加權(quán)回歸曲線
    Parameters:
        無
    Returns:
        無
    Modify:
        2018-10-24
    """
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    x_arr, y_arr = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
    y_hat_1 = lwlr_test(x_arr, x_arr, y_arr, 1.0)
    y_hat_2 = lwlr_test(x_arr, x_arr, y_arr, 0.01)
    y_hat_3 = lwlr_test(x_arr, x_arr, y_arr, 0.003)
    x_mat = mat(x_arr)
    y_mat = mat(y_arr)
    # 排序幢踏，返回索引值
    srt_ind = x_mat[:, 1].argsort(0)
    x_sort = x_mat[srt_ind][:, 0, :]
    fig, axs = plt.subplots(nrows=3, ncols=1, sharex=False, sharey=False, figsize=(10, 8))
    # 繪制回歸曲線
    axs[0].plot(x_sort[:, 1], y_hat_1[srt_ind], c='red')
    axs[1].plot(x_sort[:, 1], y_hat_2[srt_ind], c='red')
    axs[2].plot(x_sort[:, 1], y_hat_3[srt_ind], c='red')
    # 繪制原數(shù)據(jù)點
    axs[0].scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
    axs[1].scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
    axs[2].scatter(x_mat[:, 1].flatten().A[0], y_mat.flatten().A[0], s=20, c='blue', alpha=0.5)
    # 設(shè)置標(biāo)題,x軸label,y軸label
    axs0_title_text = axs[0].set_title(u'局部加權(quán)回歸曲線,k=1.0', FontProperties=font)
    axs1_title_text = axs[1].set_title(u'局部加權(quán)回歸曲線,k=0.01', FontProperties=font)
    axs2_title_text = axs[2].set_title(u'局部加權(quán)回歸曲線,k=0.003', FontProperties=font)
    # plt.setp():設(shè)置圖標(biāo)實例的屬性
    plt.setp(axs0_title_text, size=8, weight='bold', color='red')
    plt.setp(axs1_title_text, size=8, weight='bold', color='red')
    plt.setp(axs2_title_text, size=8, weight='bold', color='red')
    plt.xlabel('X')
    plt.show()


def rss_error(y_arr, y_hat_arr):
    """
    Function:
        誤差大小評價函數(shù)
    Parameters:
        y_arr - 真實數(shù)據(jù)
        y_hat_arr - 預(yù)測數(shù)據(jù)
    Returns:
        誤差大小
    Modify:
        2018-10-24
    """
    return ((y_arr - y_hat_arr) ** 2).sum()


def ridge_regress(x_mat, y_mat, lam=0.2):
    """
    Function:
        嶺回歸
    Parameters:
        x_mat - x數(shù)據(jù)集
        y_mat - y數(shù)據(jù)集
        lam - 縮減系數(shù)
    Returns:
         ws - 回歸系數(shù)
    Modify:
        2018-11-07
    """
    x_t_x = x_mat.T * x_mat
    denom = x_t_x + eye(shape(x_mat)[1]) * lam
    if linalg.det(denom) == 0.0:
        print('矩陣為奇異矩陣,不能轉(zhuǎn)置')
        return
    ws = denom.I * (x_mat.T * y_mat)
    return ws


def ridge_test(x_arr, y_arr):
    """
    Function:
        嶺回歸測試
    Parameters:
        x_mat - x數(shù)據(jù)集
        y_mat - y數(shù)據(jù)集
        lam - 縮減系數(shù)
    Returns:
         w_mat - 回歸系數(shù)矩陣
    Modify:
        2018-11-07
    """
    x_mat = mat(x_arr)
    y_mat = mat(y_arr).T
    # 行與行操作，求均值
    y_mean = mean(y_mat, axis=0)
    # 數(shù)據(jù)減去均值
    y_mat = y_mat - y_mean
    # 行與行操作凿可，求均值
    x_means = mean(x_mat, axis=0)
    # 行與行操作惑折，求方差
    x_var = var(x_mat, axis=0)
    # 數(shù)據(jù)減去均值除以方差實現(xiàn)標(biāo)準(zhǔn)化
    x_mat = (x_mat - x_means) / x_var
    # 30個不同的lambda測試
    num_test_pts = 30
    # 初始回歸系數(shù)矩陣
    w_mat = zeros((num_test_pts, shape(x_mat)[1]))
    for i in range(num_test_pts):
        # lambda以e的指數(shù)變化，最初是一個非常小的數(shù)
        ws = ridge_regress(x_mat, y_mat, exp(i - 10))
        w_mat[i, :] = ws.T
    return w_mat


def plot_ridge_regress_mat():
    """
    Function:
        繪制嶺回歸系數(shù)矩陣
    Parameters:
        無
    Returns:
        無
    Modify:
        2018-10-24
    """
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    ab_x, ab_y = load_data_set('../Machine/machinelearninginaction/Ch08/abalone.txt')
    redge_weights = ridge_test(ab_x, ab_y)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(redge_weights)
    ax_title_text = ax.set_title(u'log(lambada)與回歸系數(shù)的關(guān)系', FontProperties=font)
    ax_xlabel_text = ax.set_xlabel(u'log(lambada)', FontProperties=font)
    ax_ylabel_text = ax.set_ylabel(u'回歸系數(shù)', FontProperties=font)
    plt.setp(ax_title_text, size=20, weight='bold', color='red')
    plt.setp(ax_xlabel_text, size=10, weight='bold', color='black')
    plt.setp(ax_ylabel_text, size=10, weight='bold', color='black')
    plt.show()


def regularize(x_mat, y_mat):
    """
    Function:
        數(shù)據(jù)標(biāo)準(zhǔn)化
    Parameters:
        x_mat - x數(shù)據(jù)集
        y_mat - y數(shù)據(jù)集
    Returns:
        in_x_mat - 標(biāo)準(zhǔn)化后的x數(shù)據(jù)集
        in_y_mat - 標(biāo)準(zhǔn)化后的y數(shù)據(jù)集
    Modify:
        2018-11-11
    """
    in_x_mat = x_mat.copy()
    in_y_mat = y_mat.copy()
    in_y_mean = mean(y_mat, 0)
    in_y_mat = y_mat - in_y_mean
    in_means = mean(in_x_mat, 0)
    in_var = var(in_x_mat, 0)
    in_x_mat = (in_x_mat - in_means) / in_var
    return in_x_mat, in_y_mat


def stage_wise(x_arr, y_arr, eps=0.01, num_it=100):
    """
    Function:
        前向逐步線性回歸
    Parameters:
        x_arr - x輸入數(shù)據(jù)
        y_arr - y預(yù)測數(shù)據(jù)
        eps - 每次迭代需要調(diào)整的步長
        num_it - 迭代次數(shù)
    Returns:
        return_mat - 迭代次數(shù)
    Modify:
        2018-10-24
    """
    x_mat = mat(x_arr)
    y_mat = mat(y_arr).T
    x_mat, y_mat = regularize(x_mat, y_mat)
    m, n = shape(x_mat)
    # 初始化num_it次迭代的回歸系數(shù)矩陣
    return_mat = zeros((num_it, n))
    # 初始化回歸系數(shù)矩陣
    ws = zeros((n, 1))
    ws_test = ws.copy()
    ws_max = ws.copy()
    # 迭代num_it次
    for i in range(num_it):
        # 打印當(dāng)前回歸系數(shù)矩陣
        print(ws.T)
        lowest_error = float('inf')
        # 遍歷每個特征的回歸系數(shù)
        for j in range(n):
            for sign in [-1, 1]:
                ws_test = ws.copy()
                # 微調(diào)回歸系數(shù)
                ws_test[j] += eps * sign
                # 計算預(yù)測值
                y_test = x_mat * ws_test
                # 計算平方誤差
                rss_e = rss_error(y_mat.A, y_test.A)
                # 如果誤差更小枯跑，則更新當(dāng)前的最佳回歸系數(shù)
                if rss_e < lowest_error:
                    lowest_error = rss_e
                    ws_max = ws_test
        ws = ws_max.copy()
        # 記錄numIt次迭代的回歸系數(shù)矩陣
        return_mat[i, :] = ws.T
    return return_mat


def plot_stage_wise():
    """
    Function:
        繪制前向逐步線性回歸
    Parameters:
        無
    Returns:
        無
    Modify:
        2018-10-24
    """
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    ab_x, ab_y = load_data_set('../Machine/machinelearninginaction/Ch08/abalone.txt')
    stage_weights = stage_wise(ab_x, ab_y, 0.005, 1000)
    print(stage_weights)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(stage_weights)
    ax_title_text = ax.set_title(u'前向逐步回歸:迭代次數(shù)與回歸系數(shù)的關(guān)系', FontProperties=font)
    ax_xlabel_text = ax.set_xlabel(u'迭代次數(shù)', FontProperties=font)
    ax_ylabel_text = ax.set_ylabel(u'回歸系數(shù)', FontProperties=font)
    plt.setp(ax_title_text, size=15, weight='bold', color='red')
    plt.setp(ax_xlabel_text, size=10, weight='bold', color='black')
    plt.setp(ax_ylabel_text, size=10, weight='bold', color='black')
    plt.show()


def scrape_page(ret_x, ret_y, in_file, yr, num_pce, orig_prc):
    """
    Function:
        從html頁面讀取數(shù)據(jù)惨驶，生成ret_x和ret_y列表
    Parameters:
        ret_x - 數(shù)據(jù)x
        ret_y - 數(shù)據(jù)y
        in_file - HTML文件
        yr - 年份
        num_pce - 樂高部件數(shù)目
        orig_prc - 原價
    Returns:
        無
    Modify:
        2018-10-24
    """
    # 打開并讀取HTML文件
    with open(in_file, encoding='utf-8') as f:
        html = f.read()
    soup = BeautifulSoup(html, 'lxml')
    i = 1
    # 根據(jù)HTML頁面結(jié)構(gòu)進行解析
    current_row = soup.find_all('table', r='%d' % i)
    while (len(current_row) != 0):
        current_row = soup.find_all('table', r='%d' % i)
        title = current_row[0].find_all('a')[1].text
        lwr_title = title.lower()
        # 查找是否有全新標(biāo)簽
        if (lwr_title.find('new') > -1) or (lwr_title.find('nisb') > -1):
            new_flag = 1.0
        else:
            new_flag = 0.0
        # 查找是否已經(jīng)標(biāo)志出售，我們只收集已出售的數(shù)據(jù)
        sold_unicde = current_row[0].find_all('td')[3].find_all('span')
        if len(sold_unicde) == 0:
            print('商品 #%d 沒有出售' % i)
        else:
            sold_price = current_row[0].find_all('td')[4]
            price_str = sold_price.text
            price_str = price_str.replace('$', '')
            price_str = price_str.replace(',', '')
            if len(sold_price) > 1:
                price_str = price_str.replace('Free shipping', '')
            selling_price = float(price_str)
            if selling_price > orig_prc * 0.5:
                print("%d\t%d\t%d\t%f\t%f" % (yr, num_pce, new_flag, orig_prc, selling_price))
                ret_x.append([yr, num_pce, new_flag, orig_prc])
                ret_y.append(selling_price)
        i += 1
        current_row = soup.find_all('table', r='%d' % i)


def set_data_collect(ret_x, ret_y):
    """
    Function:
        依次讀取六種樂高套裝的數(shù)據(jù)敛助，并生成數(shù)據(jù)矩陣
    Parameters:
        ret_x - 數(shù)據(jù)x
        ret_y - 數(shù)據(jù)y
    Returns:
        無
    Modify:
        2018-10-24
    """
    file_path = './Lego/'
    scrape_page(ret_x, ret_y, file_path + 'lego8288.html', 2006, 800, 49.99)
    scrape_page(ret_x, ret_y, file_path + 'lego10030.html', 2002, 3096, 269.99)
    scrape_page(ret_x, ret_y, file_path + 'lego10179.html', 2007, 5195, 499.99)
    scrape_page(ret_x, ret_y, file_path + 'lego10181.html', 2007, 3428, 99.99)
    scrape_page(ret_x, ret_y, file_path + 'lego10189.html', 2008, 5922, 299.99)
    scrape_page(ret_x, ret_y, file_path + 'lego10196.html', 2009, 3263, 249.99)


def use_stand_regres():
    """
    Function:
        使用簡單的線性回歸進行建模
    Parameters:
        ret_x - 數(shù)據(jù)x
        ret_y - 數(shù)據(jù)y
    Returns:
        無
    Modify:
        2018-10-24
    """
    lg_x = []
    lg_y = []
    set_data_collect(lg_x, lg_y)
    # lg_y = mat(lg_y).T
    data_num, feature_num = shape(lg_x)
    lg_x1 = mat(ones((data_num, feature_num + 1)))
    lg_x1[:, 1:5] = mat(lg_x)
    ws = stand_regres(lg_x1, lg_y)
    print('%f%+f*年份%+f*部件數(shù)量%+f*是否為全新%+f*原價' % (ws[0], ws[1], ws[2], ws[3], ws[4]))


def cross_validation(x_arr, y_arr, num_val=10):
    """
    Function:
        交叉驗證嶺回歸
    Parameters:
        x_arr - x數(shù)據(jù)集
        y_arr - y數(shù)據(jù)集
        num_val - 交叉驗證次數(shù)
    Returns:
         w_mat - 回歸系數(shù)矩陣
    Modify:
        2018-10-24
    """
    # 統(tǒng)計樣本個數(shù)
    m = len(y_arr)
    # 生成索引值列表
    index_list = list(range(m))
    # create error mat 30columns num_val rows
    error_mat = zeros((num_val, 30))
    # 交叉驗證num_val次
    for i in range(num_val):
        # 訓(xùn)練集
        train_x = []
        train_y = []
        # 測試集
        test_x = []
        test_y = []
        # 打亂次序
        random.shuffle(index_list)
        # 劃分?jǐn)?shù)據(jù)集:90%訓(xùn)練集粗卜，10%測試集
        for j in range(m):
            if j < m * 0.9:
                train_x.append(x_arr[index_list[j]])
                train_y.append(y_arr[index_list[j]])
            else:
                test_x.append(x_arr[index_list[j]])
                test_y.append(y_arr[index_list[j]])
        # 獲得30個不同lambda下的嶺回歸系數(shù)
        w_mat = ridge_test(train_x, train_y)
        # 遍歷所有的嶺回歸系數(shù)
        for k in range(30):
            # 測試集
            mat_test_x = mat(test_x)
            mat_train_x = mat(train_x)
            # 測試集均值
            mean_train = mean(mat_train_x, 0)
            # 測試集方差
            var_train = var(mat_train_x, 0)
            # 測試集標(biāo)準(zhǔn)化
            mat_test_x = (mat_test_x - mean_train) / var_train
            # 根據(jù)ws預(yù)測y值
            y_est = mat_test_x * mat(w_mat[k, :]).T + mean(train_y)
            # 統(tǒng)計誤差
            error_mat[i, k] = rss_error(y_est.T.A, array(test_y))
    # 計算每次交叉驗證的平均誤差
    mean_errors = mean(error_mat, 0)
    # 找到最小誤差
    min_mean = float(min(mean_errors))
    # 找到最佳回歸系數(shù)
    best_weights = w_mat[nonzero(mean_errors == min_mean)]
    x_mat = mat(x_arr)
    y_mat = mat(y_arr).T
    mean_x = mean(x_mat, 0)
    var_x = var(x_mat, 0)
    # 數(shù)據(jù)經(jīng)過標(biāo)準(zhǔn)化，因此需要還原
    un_reg = best_weights / var_x
    print('%f%+f*年份%+f*部件數(shù)量%+f*是否為全新%+f*原價' % (
        (-1 * sum(multiply(mean_x, un_reg)) + mean(y_mat)), un_reg[0, 0], un_reg[0, 1], un_reg[0, 2], un_reg[0, 3]))


if __name__ == '__main__':
    # plot_data_set()

    # plot_regression()

    # data_mat, label_mat = load_data_set('../Machine/machinelearninginaction/Ch08/ex0.txt')
    # ws = stand_regres(data_mat, label_mat)
    # x_mat = mat(data_mat)
    # y_mat = mat(label_mat)
    # y_hat = x_mat * ws
    # print(corrcoef(y_hat.T, y_mat))

    # plot_lwlr_regression()

    # ab_x, ab_y = load_data_set('../Machine/machinelearninginaction/Ch08/abalone.txt')
    # print('訓(xùn)練集與測試集相同:局部加權(quán)線性回歸,核k的大小對預(yù)測的影響:')
    # y_hat_01 = lwlr_test(ab_x[0:99], ab_x[0:99], ab_y[0:99], 0.1)
    # y_hat_1 = lwlr_test(ab_x[0:99], ab_x[0:99], ab_y[0:99], 1)
    # y_hat_10 = lwlr_test(ab_x[0:99], ab_x[0:99], ab_y[0:99], 10)
    # print('k=0.1時,誤差大小為:', rss_error(ab_y[0:99], y_hat_01.T))
    # print('k=1  時,誤差大小為:', rss_error(ab_y[0:99], y_hat_1.T))
    # print('k=10 時,誤差大小為:', rss_error(ab_y[0:99], y_hat_10.T))
    # print('')
    # print('訓(xùn)練集與測試集不同:局部加權(quán)線性回歸,核k的大小是越小越好嗎纳击？更換數(shù)據(jù)集,測試結(jié)果如下:')
    # y_hat_1 = lwlr_test(ab_x[100:199], ab_x[0:99], ab_y[0:99], 0.1)
    # y_hat_2 = lwlr_test(ab_x[100:199], ab_x[0:99], ab_y[0:99], 1)
    # y_hat_3 = lwlr_test(ab_x[100:199], ab_x[0:99], ab_y[0:99], 10)
    # print('k=0.1時,誤差大小為:', rss_error(ab_y[100:199], y_hat_1.T))
    # print('k=1  時,誤差大小為:', rss_error(ab_y[100:199], y_hat_2.T))
    # print('k=10 時,誤差大小為:', rss_error(ab_y[100:199], y_hat_3.T))
    # print('')
    # print('訓(xùn)練集與測試集不同:簡單的線性歸回與k=1時的局部加權(quán)線性回歸對比:')
    # print('k=1時,誤差大小為:', rss_error(ab_y[100:199], y_hat_2.T))
    # ws = stand_regres(ab_x[0:99], ab_y[0:99])
    # y_hat = mat(ab_x[100:199]) * ws
    # print('簡單的線性回歸誤差大小:', rss_error(ab_y[100:199], y_hat.T.A))

    # plot_ridge_regress_mat()

    # plot_stage_wise()

    # lg_x = []
    # lg_y = []
    # set_data_collect(lg_x, lg_y)

    # use_stand_regres()

    # lg_x = []
    # lg_y = []
    # set_data_collect(lg_x, lg_y)
    # cross_validation(lg_x, lg_y)

    # lg_x = []
    # lg_y = []
    # set_data_collect(lg_x, lg_y)
    # print(ridge_test(lg_x, lg_y))

    from sklearn import linear_model
    reg = linear_model.Ridge(alpha=0.5)
    lg_x = []
    lg_y = []
    set_data_collect(lg_x, lg_y)
    reg.fit(lg_x, lg_y)
    print('%f%+f*年份%+f*部件數(shù)量%+f*是否為全新%+f*原價' % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3]))