機器學習——線性回歸的原理，推導過程，源碼，評價

0.線性回歸

　　做為機器學習入門的經典模型，線性回歸是絕對值得大家深入的推導實踐的，而在眾多的模型中，也是相對的容易。線性回歸模型主要是用於線性建模，假設樣本的特征有n個，我們通常將截距項也添加到特征向量x中，即在x中添加一個全為1的列，這是，我們就能夠將模型表示為如下的形式：

1.殘差的解釋

　　根據上述的模型，我們可以表示出樣本的標簽值與模型預測值之間的表達式，如下所示：

上述式子中，根據殘差的定義：實際值和預測值之間的差值，可知，即為模型的殘差。那么，我們想知道的是，在模型中是怎么樣的分布呢？

　　1，殘差是由模型中的許多誤差的積累的結果，即模型中許多的誤差的累加作用的結果。

　　2，假定這些誤差的分布是相同的。

　　那么，根據中心極限定理，是多個獨立同分布的變量的累加結果，則是服從均值為某值，方差為某值的高斯分布，對於均值，我們總是可以通過改變模型的截距，是的模型得到的平面上下移動，是的得到的殘差的分布的均值為0，假設方差為定值。

 2.中心極限定理

　　在實際的問題中，很多的隨機現象可以看做是眾多因素的獨立影響的綜合反應，往往可以看做是服從正態分布。比如，城市的耗電量：大量用戶的用電量總和。測量誤差：許多觀察不到的，微小誤差的總和。中心極限定理的關鍵是多個隨機變量的和，在有些問題中是乘積誤差，這時則需要鑒別后再使用。

3.極大似然估計

　　我們得到了殘差的分布函數，因而我們可以對殘差進行極大似然估計，就能夠得到似然函數，而似然函數中應該是含有參數，因而我們就能夠通過似然函數求極大值，得到參數的表達式，進而得到模型的解。

　　在上述推導過程中，我們先是對似然函數求極大值，得到的結果中，發現是一個減法計算，因而只需對后面的式子求最小值，則能夠得到線性回歸的代價函數，這也和我們的理解是相符合的，即模型的預測值應該使得它和實際值相差越小越好。

 4.解析解

　　在得到的代價函數中，計算極小值，我們發現這是一個復合函數，二次函數和線性函數的復合函數，根據凸函數的性質：凸函數和仿射函數的復合函數還是凸函數，因而代價函數是凸函數，存在極小值也是最小值。我們通過求極值點就能夠得到最優解。

　　我們可以將M個N維樣本組成矩陣X，m行n列，每一行對應一個樣本，每一列對應所有樣本的一個維度，由於還有截距項，因而有一列全為1。

　　計算梯度

　　可得到參數的解析解：

　　當上述式子不可逆時，可以添加擾動，使其是可逆的

　　看到上述的式子，自然想到了正則化，我們在代價函數里面添加正則項，進行求解剛好得到上述的式子。

L2正則化

 L1正則化

 　　對於L2正則我們可以計算梯度，但是L1正則中我們是無法計算梯度的，一種可行的辦法是使用泰勒公式對后半部分進行近似計算。

 5.梯度下降

　　除了用解析解，在機器學習中，我們更多的時候使用的是學習方法，通過優化的方式得到最優解，下面我們是用梯度下降來進行模型的求解。

　　得到了梯度之后，我們可以通過梯度下降的方式不斷更新參數得到最優解。

6.模型的評價

　　對於m個樣本，通過得到的模型，我們可以計算出m估計值，下面我們定義總平方和TSS：

　　得到了總平方和后，我們可以計算殘差平方和SSE

　　由得到的總平方和與殘差平方和，我們定義R方統計量：

　　R方的值越大，擬合的效果越好，最優值是1。

 7.源碼

def getLinearData(lamda):
    x1 = np.random.random(50) * 20 - 10
    x2 = np.random.random(50) * 7 - 2
    x3 = np.random.random(50) * 36 - 13
    x4 = np.random.random(50) * 9 - 3
    x = np.stack((x1, x2, x3, x4), axis=1)
    b = np.ones((50))
    x = np.c_[x, b]
    y = np.sum(x*lamda, axis=1) + np.random.random(50) * 50 - 25
    return x, y


def predict_y(x, k):
    return np.sum(x*k, axis=1)


def cal_loss(x, y, k):
    y_predict = np.sum(x*k, axis=1)
    loss = 0.5 * np.sum(np.square(y - y_predict))
    return loss


def mini_batch_GD():
    return


def BGD(x, y, k, n, s, loss):
    losslist = []
    for i in range(n):

        y_hat = predict_y(x, k)
        for j, dk in enumerate(k):
            xi = x[:, j]
            dk = dk + s * np.sum((y-y_hat) * xi, axis=0)
            k[j] = dk
        new_loss = cal_loss(x, y, k)
        print(new_loss)
        losslist.append(new_loss)
        if new_loss < loss or (i > 10 and np.abs(new_loss - losslist[-2])< 0.0001):
            return k, new_loss, losslist
    new_loss = cal_loss(x, y, k)
    return k, new_loss, losslist


def SGD(x, y, k, n, s, loss):
    losslist = []
    for i in range(n):

        for j, xi in enumerate(x):
            y_predict = np.sum(xi*k)
            k = k + s * (y[j] - y_predict) * xi
            new_loss = cal_loss(x, y, k)
            losslist.append(new_loss)
        print(new_loss)

        if new_loss < loss :
            return k, new_loss, losslist
    new_loss = cal_loss(x, y_train, k)
    return k, new_loss, losslist


def linear_Regression(x, y):
    size = np.shape(x)[1]
    k = np.ones(size)
    # k, loss, losslist = BGD(x, y, k, 20000000, 0.000005, 10)
    k, loss, losslist = SGD(x, y, k, 20, 0.0005, 10)
    return k, loss, losslist


if __name__ == '__main__':
    lamda = [2, -5, 3, 7, 15]  # a1, a2, a3, a4, b
    x, y = getLinearData(lamda)
    x_train, y_train = x[0:40, :], y[0:40]
    x_test, y_test = x[40:, :], y[40:]
    k, loss, losslist = linear_Regression(x_train, y_train)
    print(k)
    print(loss)
    i = np.arange(np.size(losslist))
    plt.figure()
    # plt.plot(i, losslist, 'r-', label='BGD')
    plt.plot(i, losslist, 'r-', label='SGD')
    plt.legend()
    plt.show()

　　模型運行的結果：

　　在寫梯度下降的代碼中，由於是用所有的樣本計算梯度進行更新的，所以每次得到的梯度值很大，開始時我給出的步長是0.1-0.01，結果出現了跨越極值的情況，最后的到無窮大的結果。以為是代碼的問題，找半天代碼一直沒想通，后面才想到是因為學習率太大的問題。