一、回歸介紹

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

HorsePower = 0.0015 * annualSalary - 0.99 * hoursListeningToPublicRadio

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

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

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

全文代碼在最后第八點(diǎn)里。

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

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

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

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

平方誤差和可以寫做：

用矩陣表示還可以寫做：

找到w，使平方誤差和最小。因?yàn)槲覀冋J(rèn)為平方誤差和越小，說明線性回歸擬合效果越好。

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

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

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

值得注意的是，上述公式中包含逆矩陣，也就是說，這個(gè)方程只在逆矩陣存在的時(shí)候使用，也即是這個(gè)矩陣是一個(gè)方陣，并且其行列式不為0。

上述的最佳w求解是統(tǒng)計(jì)學(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方法計(jì)算這兩個(gè)序列的相關(guān)系數(shù)來比較預(yù)測值和真實(shí)值的相關(guān)性。

可以看到，對角線上的數(shù)據(jù)是1.0，因?yàn)閥Mat和自己的匹配是完美的，而YHat和yMat的相關(guān)系數(shù)為0.98。

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

線性回歸的一個(gè)問題是有可能出現(xiàn)欠擬合現(xiàn)象，因?yàn)樗蟮氖蔷哂凶钚【秸`差的無偏估計(jì)。顯而易見，如果模型欠擬合將不能取得最好的預(yù)測效果。所以有些方法允許在估計(jì)中引入一些偏差，從而降低預(yù)測的均方誤差。

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

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

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

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

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

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

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

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

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

如果數(shù)據(jù)的特征比樣本點(diǎn)還多，就不可以使用線性回歸和之前的方法來做預(yù)測了，這是因?yàn)樵谟?jì)算 (X^T X)^ -1的時(shí)候會(huì)出錯(cuò)。如果特征比樣本點(diǎn)還多（n > m），也就是說輸入數(shù)據(jù)的矩陣X不是滿秩矩陣。非滿秩矩陣在求逆時(shí)會(huì)出現(xiàn)問題。

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

4.1嶺回歸

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

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

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

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

4.2前向逐步線性回歸

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

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

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

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

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

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

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

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

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

這里隨機(jī)選取樣本，因?yàn)槠潆S機(jī)性，所以每次運(yùn)行的結(jié)果可能略有不同。不過整體如上圖所示，可以看出，它與常規(guī)的最小二乘法，即普通的線性回歸沒有太大差異。我們本期望找到一個(gè)更易于理解的模型，顯然沒有達(dá)到預(yù)期效果。

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

看運(yùn)行結(jié)果的第一行，可以看到最大的是第4項(xiàng)，第二大的是第2項(xiàng)。

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

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

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

sklearn.linear_model提供了很多線性模型，包括嶺回歸、貝葉斯回歸、Lasso等?，F(xiàn)在我們用sklearn實(shí)現(xiàn)嶺回歸。

我正則化項(xiàng)系數(shù)設(shè)為0.5，其余參數(shù)使用默認(rèn)?？梢钥吹剑@得的結(jié)果與上面的結(jié)果類似。

七、小結(jié)：

與分類一樣，回歸也是預(yù)測目標(biāo)值的過程?；貧w與分類的不同點(diǎn)在于，前者預(yù)測連續(xù)類型變量，而后者預(yù)測離散類型變量。

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

當(dāng)數(shù)據(jù)的樣本數(shù)比特征數(shù)還少時(shí)候，矩陣 X^TX的逆不能直接計(jì)算。即便當(dāng)樣本數(shù)比特征數(shù)多時(shí)， X^TX的逆仍有可能無法直接計(jì)算，這是因?yàn)樘卣饔锌赡芨叨认嚓P(guān)。這時(shí)可以考慮使用嶺回歸，因?yàn)楫?dāng) X^TX的逆不能計(jì)算時(shí)，它仍保證能求得回歸參數(shù)。

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

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

八、貼上全部代碼：

# -*- 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:
        計(jì)算回歸系數(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)線性回歸計(jì)算回歸系數(shù)w
    Parameters:
        test_point -  測試樣本點(diǎn)
        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
    # 計(jì)算回歸系數(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ù)點(diǎn)
    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)實(shí)例的屬性
    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:
        誤差大小評價(jià)函數(shù)
    Parameters:
        y_arr - 真實(shí)數(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ù)減去均值除以方差實(shí)現(xiàn)標(biāo)準(zhǔn)化
    x_mat = (x_mat - x_means) / x_var
    # 30個(gè)不同的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ù)變化，最初是一個(gè)非常小的數(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')
        # 遍歷每個(gè)特征的回歸系數(shù)
        for j in range(n):
            for sign in [-1, 1]:
                ws_test = ws.copy()
                # 微調(diào)回歸系數(shù)
                ws_test[j] += eps * sign
                # 計(jì)算預(yù)測值
                y_test = x_mat * ws_test
                # 計(jì)算平方誤差
                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 - 原價(jià)
    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)進(jìn)行解析
    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:
        使用簡單的線性回歸進(jìn)行建模
    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*原價(jià)' % (ws[0], ws[1], ws[2], ws[3], ws[4]))


def cross_validation(x_arr, y_arr, num_val=10):
    """
    Function:
        交叉驗(yàn)證嶺回歸
    Parameters:
        x_arr - x數(shù)據(jù)集
        y_arr - y數(shù)據(jù)集
        num_val - 交叉驗(yàn)證次數(shù)
    Returns:
         w_mat - 回歸系數(shù)矩陣
    Modify:
        2018-10-24
    """
    # 統(tǒng)計(jì)樣本個(gè)數(shù)
    m = len(y_arr)
    # 生成索引值列表
    index_list = list(range(m))
    # create error mat 30columns num_val rows
    error_mat = zeros((num_val, 30))
    # 交叉驗(yàn)證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個(gè)不同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)計(jì)誤差
            error_mat[i, k] = rss_error(y_est.T.A, array(test_y))
    # 計(jì)算每次交叉驗(yàn)證的平均誤差
    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*原價(jià)' % (
        (-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時(shí),誤差大小為:', rss_error(ab_y[0:99], y_hat_01.T))
    # print('k=1  時(shí),誤差大小為:', rss_error(ab_y[0:99], y_hat_1.T))
    # print('k=10 時(shí),誤差大小為:', 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時(shí),誤差大小為:', rss_error(ab_y[100:199], y_hat_1.T))
    # print('k=1  時(shí),誤差大小為:', rss_error(ab_y[100:199], y_hat_2.T))
    # print('k=10 時(shí),誤差大小為:', rss_error(ab_y[100:199], y_hat_3.T))
    # print('')
    # print('訓(xùn)練集與測試集不同:簡單的線性歸回與k=1時(shí)的局部加權(quán)線性回歸對比:')
    # print('k=1時(shí),誤差大小為:', 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*原價(jià)' % (reg.intercept_, reg.coef_[0], reg.coef_[1], reg.coef_[2], reg.coef_[3]))