機器學習基本算法系列之線性迴歸

寫在前面：我將從一個入門者的視角（水平）將機器學習中的經常使用算法娓娓道來。自身水平確實有限，若是其中有什麼錯誤的話但願你們指出，避免誤導你們。對於開篇有比較多的東西想說，因此另取一章介紹一下我我的對機器學習這個領域的理解和一些基本常識介紹。若是你時間比較充裕，能夠看看我前言的文字介紹。

0 前言

介紹一下我的想法

0.1 關於機器學習

毫無疑問，機器學習如今愈來愈熱門了，需求隨着之後人工智能大爆發會持續增加。無論你是否是這個編程領域的人，都有必要了解一下機器學習是怎麼回事，至少有個理性的認識，而不是以爲很神祕。可是涉及到算法數學知識是能夠忽略的，只要知道，咱們經過某個數學手段能夠達到某個目的。好比這裏要介紹的線性迴歸在更新權重時是依據梯度降低算法的，你能夠沒必要知道什麼是梯度降低，知道知道數學上對於求解最優化是有本身的手段的，固然若是能知道更好。因爲機器學習的性質，速成基本不太可能，並且如今處於高速發展的階段，必須保持學習的熱情才能學得下去。另外，這個前期的學習成本是比較高的，在本身寫代碼以前，要想內心有點底，至少要對 微積分、線性代數 和 機率論 的應用要有點本身的理解的。因此這也是一大優點，不太容易被取代。

0.2 關於線性迴歸

對於理性瞭解機器學習的捷徑就是弄懂線性迴歸算法的整個流程，若是你經過這篇文章真的理解了線性迴歸，那麼你應該就理解這句話了。固然，對於大多數閱讀者篇文章的人來講，可能早就對線性迴歸有了本身的理解。可是我這要提的是，若是你沒打算入門機器學習，閱讀一下線性迴歸也是對本身知識面進行拓展的一個好方式，由於基本大多數人都是從這裏開始的。

0.3 關於基礎

首先，你須要一點 基礎。若是沒有也簡單，自學一下就行。而後就是要學習一下如下工具的使用：NumPy，Matplotlib，sklearn，TensorFlow。其中部分能夠參考我之前的文章。

• 從 TensorFlow 入門機器學習
• NumPy 官方快速入門教程(譯)
• Matplotlib 基本操做
• Matplotlib 保存 GIF 動圖

有了這些你還須要一些數學基礎。數學的話我以爲不要求全學會，也不太現實，能夠在寫代碼過程當中碰到什麼算法再去尋求數學證實，前提是以前的知識量最好能達到看懂大多數數學知識的水平。這裏我以爲知乎上有些數學科普其實挺不錯的，畢竟不是數學系，因此理解主要邏輯以後持有一種拿來用的態度我以爲沒問題。

• 機器學習之重拾機率論

0.4 關於參考

這裏很感謝像吳恩達、莫煩等老師的無私奉獻，還有不少在網上寫博客分享知識的人，讓編程這個行業造成一個良好的學習風氣。並且網上還流傳不少優質的學習資源，因此這裏我將我學習過程的參考資料整理進了 GitHub，我想你們奉獻出來就是爲了別人更好的學習吧，若是有侵權之類的話，我會當即刪除的。而後這一系列的文章，我會主要以吳恩達老師的課程爲基礎，以初學者的視角記錄而且實踐。

• 學習參考書籍

0.5 關於我

一個即將畢業的知識水平有限的本科生，因此文章確定有不少不足之處，望你們指正。而後實驗的代碼在這裏。

1 背景

故事就不講了，能夠參考視視頻。下面我講講這個是要解決什麼問題？生活中每每有不少現象，而咱們能夠從現象總結出某些結論，可是，若是咱們看到了一個以前歷來都沒有看到過的現象，咱們就沒法下結論了。但可是一些有經驗的人能夠憑藉豐富的經驗進行預測。就好比，一個處於懵懂期的小孩子，看見一片烏雲不能推測出會下雨，而一個大人就會知道，提醒小孩出門要帶傘。。。好吧，仍是在講故事，那就繼續吧😂😂😂。這個例子中的大人小孩的區別在於年齡不一樣，見識不一樣，知識水平天然不一樣，能作出的判斷也天然不一樣。因此要想達到對現象進行預測，那麼就必須進行學習，等小孩慢慢長大了，看到螞蟻搬家就會帶傘出門了。

好了，咱們總結一下上面的故事：現象能夠比喻成數據，結論形成的行爲也是數據，惟一不一樣的是，是由現象得出的結論，因此能夠將現象的數據看做是輸入，結論的數據當作是輸出。那麼機器學習就是將輸入數據輸入進一個模型，而後將模型的輸出結果和以前的「正確」結論做比對慢慢更改模型，直到這個模型對出入數據有較好的預測，那麼這個模型就能夠拿來用了。抽象出來的話就是找一個輸入 和輸出 之間的映射，使得對於絕大多數 都能經過映射 較好地獲得 。

房價的例子我就不舉了，直接看下面的圖，看看能不能解決

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt

def draw():
    X = np.linspace(0, 10, 50)
    noise = np.random.normal(0, 0.5, X.shape)
    Y = X * 0.5 + 3 + noise

    plt.scatter(X, Y)
    plt.savefig("scattergraph.png")
    plt.show()

if __name__ == "__main__":
    draw()
複製代碼

好了，你可能會發現，要想給出 依賴的是直覺，因此每個人的答案都是不同的，而對於某些標準統一的領域（好比軍事項目），這種依靠直覺判斷的狀況絕對不允許發生。想象在現代信息化戰場上靠直覺怎麼可能做戰成功？

因此咱們須要一個模型，你能夠理解成一個黑盒，把 喂進去，把 取出來，並且不管誰，只要那個黑盒是同一個，那麼每一個人的輸入 相同的話，輸出 就是相同的。 那麼在這裏這個模型大概是什麼形狀的呢？目前這個仍是咱們人爲干預進行選擇的。直觀來講，能夠這個分佈當作一條很粗的線，那麼咱們選擇線性模型來模擬分佈。那麼如何肯定這個線性模型呢？這裏就涉及到今天的主角：線性迴歸

2 一元線性迴歸

先從一元線性函數開始分析，既然假設是一元線性，那麼咱們的擬合函數就能夠定義爲 了。一樣那前一個圖做爲例子，咱們如今的目的是求出，若是求出來了，給定任意一個 就能求出 來了。

2.1 定義評估函數-損失函數

這裏講解如何求出 的整個思考過程。

咱們選取一個點 由於偏差的存在，因此 。同理若是有不少點，那麼就有如下結果，爲了方便，咱們取評估值 。

咱們的目的天然是偏差越小越好，因此採用平方和的形式表達偏差：

到這裏，咱們來看看圖像是怎麼樣的？

PS：圖有點醜，實在是沒辦法了，設置了噪點，使得函數值跨度很大，稍微改一下就會呼臉。只能設置散點圖湊合看了，望知道的小夥伴告知。好了，就不在這浪費過多時間了，直接看代碼。

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel(r'$Weight$')
ax.set_ylabel(r'$bias$')
ax.set_zlabel(r'$loss\ value$')
N = 200

X = np.linspace(0, 10, N)
noise = np.random.normal(0, 0.1, X.shape)
Y = X * 0.5 + 3 + noise

W = np.linspace(-10, 10, N)
b = np.linspace(-2, 8, N)

p = 0.6
plt.xlim(W.min() * p, W.max() * p)
plt.ylim(b.min() * p, b.max() * p)
ax.set_zlim(0, 1000)

h = []
for i in W:
    for j in b:
        h.append(np.sum(np.square(Y - (X * i + j))))

print(np.min(h))
h = np.asarray(h).reshape(N, N)
W, b = np.meshgrid(W, b)

ax.scatter(W, b, h, c='b')
plt.savefig("lossgraph.png")
plt.show()
複製代碼

如今咱們經過可視化的手段知道了大概什麼位置取最小值，其實咱們仍是沒有達到目的，咱們要的是 具體值。

2.2 求解最小值-梯度降低

問題被進一步細化，就是求解下面二元函數在何時取極值：

這時候高數就登場了，找了一張圖

那問題就很簡單了，就是直接求對 的偏導等於 的方程組。可是這是數學方法，不符合計算機思惟，就好像咱們解一個一元一次方程，普通方法是直接拿定義域內的值一個一個試，直到結果符合預期。固然，高級一點點的方法是能夠寫一個解釋器，把人類計算方法程序化。顯然這裏不太現實，並且方程類型一變，解釋器極可能就不適用了。因此咱們仍是採用那種機器承認的「笨」辦法，想起開學時計算機導論老師評價計算機的一句話：fast but stupid。說明在這種狀況下，咱們能夠接受用快來彌補其它劣勢，只要計算的值一步一步靠近最終結果就能知足需求。其實說到這裏，計算機變快的這些特色和如今機器學習、人工智能領域的興起有必定關係。

好了，下面開始想一個辦法如何讓結果慢慢趨近與極值。

咱們能夠這樣想，隨機放一個球在這個平面上，等到球禁止時，它所在的位置就是極值的位置。那咱們如何模仿這一過程呢？

第一步：放球

這個好作，就是直接隨機一個在定義域內的任意位置就行了。

第二步：滾動

滾動的話如何去模擬是一個比較難的點，可是仔細分析就很好理解了。你能夠把本身想象成那個球，開始站在某個地方，這個地方很不穩，那麼你天然反應確定是往平整的地方去，前提是每一個位置的海拔差很少。那這樣就說明那些斜率越大的地方有更大的概率更快跌落谷底。這個例子有個不恰當的地方，能夠思考下二次函數，這裏主要是體會這種思惟，而咱們是能夠變通的。當你有了這種斜率的思惟，那就好辦了，下面給出總結：咱們老是但願在給定步長的狀況下，往水平最低的地方行進。而其中要知道那個是最低，又是一個難題，計算機可沒有直覺，若是數值設置不當，可能致使函數收斂過慢或者直接發散，而這些都是機器學習應該極力避免的。

爲了好分析問題，咱們採用控制變量方法。能夠看到，在損失函數中有兩個變量 ，當 取某一值得時候，是一個二次函數。

# coding: utf-8

import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()

N = 200

X = np.linspace(0, 10, N * 2)
noise = np.random.normal(0, 0.1, X.shape)
Y = X * 0.5 + 3 + noise

W = 1
b = np.linspace(-5, 6, N)

h = []
for i in b:
    h.append(np.sum(np.square(Y - (X * W + i))))

plt.plot(b, h)

plt.savefig("quadraticFunction.png")

plt.show()
複製代碼

咱們來直接看吳恩達老師的課件吧，有點不想作圖了。。。

如今咱們能夠引出更新規則了，在這裏求解最小值，對於一元函數來講就是求導，對於多元函數就是求偏導了。很明顯，當斜率爲正的時候，咱們要往負反向走，反之往正方向走。因此咱們能夠加上一個關於斜率呈反比的函數來跟新值，至於更新力度也就是前面說的一次走多遠就是學習率的設定了。並且離極值點越遠，斜率絕對值越大，步子邁得越大，符合更新邏輯。

這裏可能有點思惟跳躍的是除了一個 ，能夠從數據量的角度思考，也能夠從更新斜率的角度思考，由於在更新權重中，咱們不只僅只是計算一個點的斜率，若是直接求和一定會致使權重過大從而最終結果是發散的，這點你們能夠本身改改代碼試試。好了，既然基本思路都出來了，那就寫代碼實現吧。

# coding: utf-8

import matplotlib.pyplot as plt
import numpy as np

N = 200

X = np.linspace(0, 10, N * 2)
noise = np.random.normal(0, 0.5, X.shape)
Y = X * 0.5 + 3 + noise


def calcLoss(train_X, train_Y, W, b):
    return np.sum(np.square(train_Y - (train_X * W + b)))

def gradientDescent(train_X, train_Y, W, b, learningrate=0.001, trainingtimes=500):
    global loss
    global W_trace
    global b_trace
    size = train_Y.size
    for _ in range(trainingtimes):
        prediction = W * train_X + b
        tempW = W + learningrate * np.sum(train_X * (train_Y - prediction)) / size
        tempb = b + learningrate * np.sum(train_Y - prediction) / size
        W = tempW
        b = tempb
        loss.append(calcLoss(train_X, train_Y, W, b))
        W_trace.append(W)
        b_trace.append(b)


Training_Times = 100
Learning_Rate = 0.002

loss = []
W_trace = [-1]
b_trace = [1]
gradientDescent(X, Y, W_trace[0], b_trace[0], learningrate=Learning_Rate, trainingtimes=Training_Times)
print(W_trace[-1], b_trace[-1])

fig = plt.figure()
plt.title(r'$loss\ function\ change\ tendency$')
plt.xlabel(r'$learning\ times$')
plt.ylabel(r'$loss\ value$')
plt.plot(np.linspace(1, Training_Times, Training_Times), loss)
plt.savefig("gradientDescentLR.png")
plt.show()
複製代碼

到這裏，咱們基本手寫實現了一元變量的線性迴歸，這裏爲了方便訓練次數取得比較少，只有 次，這遠遠不夠，我測試了一下，大概達到 次訓練效果就比較好了。固然也能夠調整學習率，初始值……讀者能夠自行嘗試。

3 多元線性迴歸

這裏咱們先拆分理解多元表明多個自變量，線性表明自變量之間知足齊次性和可加性，迴歸就是一種方法的表述。一元好辦，咱們能夠經過各類工具實現可視化，可是多元的表述就有了一些困難，特別是在高維線性空間咱們想象不出來。因此咱們就必須抽象出來，使用數學符號代替那些複雜的變量。

可能你已經想到了，對，沒錯，咱們要使用線性代數了。若是以前沒有一點基礎的話，能夠稍微花點時間在知乎上了解下線性代數是作什麼的就行了，這裏只涉及線性代數的基礎認識應用。

假如當咱們有三個自變量 和一個因變量 ，且符合線性關係。那麼就會有如下方程成立：

再假設咱們有四組數據，那麼有：

因而咱們能夠抽象如下，寫成矩陣的形式：

再抽象：

其中

咱們的目的是求解 ，也就是把方程 中的 解出來。直接解個方程就能獲得咱們想要的效果，太好了，立刻開始：

若是這裏不理解的話，我簡單說一下，求逆運算必須是形如 的形式，也就是行和列要相等。爲何這樣呢？個人理解是矩陣是對一個線性空間到另外一個線性空間的映射，因此若是一個矩陣的秩小於行數，也就是裏面有線性相關的向量，在對線性空間進行變換時，可能進行降維影響，也就是形成某個維度塌縮，這種影響是不可逆的，因此這種矩陣是沒有逆矩陣進行恢復的。若是要從一個向量映射到另外一個向量，還有從另外一個向量映射回來的話，這個矩陣必須有逆矩陣也就是滿秩的。還有個東西叫奇異矩陣，感興趣的能夠了解一下。

回到主線，直接求解方程吧。

# coding: utf-8

import numpy as np

X1 = np.asarray([2104, 1416, 1534, 852]).reshape(4, 1)
X2 = np.asarray([5, 3, 3, 2]).reshape(4, 1)
X3 = np.asarray([1, 2, 2, 1]).reshape(4, 1)

X = np.mat(np.column_stack((X1, X2, X3, np.ones(shape=(4, 1)))))
noise = np.random.normal(0, 0.1, X1.shape)
Y = np.mat(2.5 * X1 - X2 + 2 * X3 + 4 + noise)
YTwin = np.mat(2.5 * X1 - X2 + 2 * X3 + 4)

W = (X.T * X).I * X.T * Y
WTWin = (X.T * X).I * X.T * YTwin
print(W, "\n", WTWin)

# output:
# [[ 2.50043958]
# [-1.16868808]
# [ 1.79213736]
# [ 4.27637958]] 
# [[ 2.5]
# [-1. ]
# [ 2. ]
# [ 4. ]]
複製代碼

這裏咱們採用吳恩達老師的房子數據，本身生成的數據之間有太大的類似性，算出的結果偏差太大。

咱們基本能夠計算多元線性迴歸，可是一點也體現不出機器學習中學習這個詞，固然咱們能夠本身利用前面給出的例子利用 函數求解。這裏咱們藉助 幫咱們完成。

# coding: utf-8

import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt

N = 1000
train_X1 = np.linspace(1, 10, N).reshape(N, 1)
train_X2 = np.linspace(1, 10, N).reshape(N, 1)
train_X3 = np.linspace(1, 10, N).reshape(N, 1)
train_X4 = np.linspace(1, 10, N).reshape(N, 1)

# train_X = np.column_stack((train_X1, np.ones(shape=(N, 1))))
train_X = np.column_stack((train_X1, train_X2, train_X3, train_X4, np.ones(shape=(N, 1))))

noise = np.random.normal(0, 0.5, train_X1.shape)
# train_Y = 3 * train_X1 + 4
train_Y = train_X1 + train_X2 + train_X3 + train_X4 + 4 + noise

length = len(train_X[0])

X = tf.placeholder(tf.float32, [None, length], name="X")
Y = tf.placeholder(tf.float32, [None, 1], name="Y")

W = tf.Variable(np.random.random(size=length).reshape(length, 1), dtype=tf.float32, name="weight")

activation = tf.matmul(X, W)
learning_rate = 0.006

loss = tf.reduce_mean(tf.reduce_sum(tf.pow(activation - Y, 2), reduction_indices=[1]))
optimizer = tf.train.GradientDescentOptimizer(learning_rate).minimize(loss)

training_epochs = 2000
display_step = 100

loss_trace = []

with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    for epoch in range(training_epochs):
        sess.run(optimizer, feed_dict={X: train_X, Y: train_Y})
        temp_loss = sess.run(loss, feed_dict={X: train_X, Y: train_Y})
        loss_trace.append(temp_loss)
        if 1 == epoch % display_step:
            print('epoch: %4s'%epoch, '\tloss: %s'%temp_loss)
    print("\nOptimization Finished!")
    print("\nloss = ", loss_trace[-1], "\nWeight =\n", sess.run(W, feed_dict={X: train_X, Y: train_Y}))


plt.plot(np.linspace(0, 100, 100), loss_trace[:100])
plt.savefig("tensorflowLR.png")
plt.show()

# output:
# epoch: 1 loss: 118.413925
# epoch: 101 loss: 1.4500043
# epoch: 201 loss: 1.0270562
# epoch: 301 loss: 0.75373846
# epoch: 401 loss: 0.5771168
# epoch: 501 loss: 0.46298113
# epoch: 601 loss: 0.38922414
# epoch: 701 loss: 0.34156123
# epoch: 801 loss: 0.31076077
# epoch: 901 loss: 0.29085675
# epoch: 1001 loss: 0.27799463
# epoch: 1101 loss: 0.26968285
# epoch: 1201 loss: 0.2643118
# epoch: 1301 loss: 0.26084095
# epoch: 1401 loss: 0.2585978
# epoch: 1501 loss: 0.25714833
# epoch: 1601 loss: 0.25621164
# epoch: 1701 loss: 0.2556064
# epoch: 1801 loss: 0.2552152
# epoch: 1901 loss: 0.2549625
# Optimization Finished!
# loss = 0.25480175 
# Weight =
# [[1.0982682 ]
# [0.9760315 ]
# [1.0619627 ]
# [0.87049955]
# [3.9700394 ]]
複製代碼

這裏的擬合效果不太好，不知道是否是數據的問題，由於數據增加的類似性過高，感受可能過擬合了，若是有知道的小夥伴歡迎告知。下圖能夠看到 早早就收斂了。

4 總結

線性迴歸的求解過程很直觀，符合咱們的理解。對於這個問題有一個先入爲主的觀點就是數據必定是擬合成線性的，由於這裏講的是線性迴歸。固然，數據分佈不是線性的這個方法就不適用了，一樣，若是咱們看不到數據的分佈，不能對模型進行預估，那隻能一步一步去試探，因此機器學習開始就是選取適用於當前數據的模型，能夠理解爲找一個輸入輸出的合理映射關係，而後導入數據，構造一個可以正確評估擬合效果的損失函數，接着就是對損失函數進行最優化，這裏有一個問題前面沒有說起的是若是隻是隨機地找，其實咱們不能保證找到的就是全局最優，固然通常狀況下，局部最優獲得的模型已經能很好的預測輸出了。這只是個人我的理解啦，歡迎你們一塊兒討論，而後實驗的代碼在個人 GitHub 上，須要的能夠自取。而後這個系列可能會繼續更新，不過離下一次更新可能要點時間，由於新年快樂^_^，對了參考資料前面又給出，就不一一列舉感謝了。