機器學習-預測-線性系統的預測（最小二乘法、正規方程式實現）

機器學習-預測-線性系統的預測

現在預測學的核心概念：回歸。從數學的角度，為事物（系統）的預測提供現代的技術方法。

 

回歸與現代預測學

統計學上最初回歸的含義由高爾頓（達爾文的表弟）通過研究父母身高與孩子身高得出。

矮個父母所生的兒子往往會比其父母更高，高個父母所生兒子的身高卻回降到多數人的平均身高。也就是說，當父母身高走向極端時，子女的身高不會進一步極端化，他們的身高要比父母更接近平均身高，即有”回歸“到平均數的趨勢，這就是統計學上最初回歸的含義。

高爾頓把這一現象叫作”向平均數方向的回歸“(Regression Toward Mediocrity)

回歸反應了系統的隨機運動總是趨向於其整體運動規律的趨勢。在數學上，就是根據系統的總體靜態觀測值，通過算法去除隨機性的噪聲，發現系統整體運動規律的過程。

回歸是現代預測學的基礎，對各種數據集的預測結果，可以理解為對其回歸函數的計算。

最小二乘法

最小二乘法是目前已知最古老的回歸預測方法，由高斯於1809年發表於《天體運動論》中。簡單地說，最小二乘的思想就是使觀測點和估計點的距離平方和最小。

古漢語中，”平方“稱為”二乘“，這里是指用平方值來度量預測點與估計點的遠近（這就消除了距離中的正負號），”最小“指的是參數的估計值要保證各個觀測點與估計點的距離平方和達到最小。

解得Y=aX+b，就稱為該樣本集的最小二乘法，求解這個函數的解法也稱為最小二乘法。

對於某個數據集(xi,yi)(i=0,1,...,n)，我們需要找到一條趨勢線，能夠表達出數據集(xi, yi)這些點所指向的方向。

我們現用一個直線函數表示這條曲線：Y = aX + b

數據集的點一定位於這條趨勢線的上下兩側，或者與趨勢線重合。我們把某個樣本點xi到這條趨勢線的垂直距離定義為殘差，那么過這一點與趨勢線平行的樣本函數為

 

如果這個樣本點位於趨勢線的上側，在(殘差i) > 0,反之則(殘差i) < 0,如果樣本點位於趨勢線上則(殘差i) = 0。

求解這條趨勢線，因為是線性函數，所以也就是求解a、b這兩個值。

 

因為殘差有正負號的問題，所以統一用平方和來計算，即殘差平方和；

 

很明顯這個二次函數式一個凸函數(單峰函數)，接下來對該函數求極值，即它的一階導數等於0.

 

 

 

兩個方程聯立，令

 

解得a、b的值為

 

 

 

xi和yi均為已知的樣本集

解得直線函數Y=aX+b，就稱為該樣本集的最小二乘解，求解這個函數的解法也稱為最小二乘法。

 

正規方程組法

假設所有的樣本點構成一個線性方程組矩陣，其形式化的表示：aX+b = 0，這里未知數是a、b,X是已知的樣本向量矩陣。

為了簡化矩陣的形式，我們從X中抽取最后一列Xn,作為輸出的列向量。在X的第一列前增加一列b,這樣b就被合並到X中，改變一下方程組的形式為：Y = XA，因為b被合並到X中，這里的A就包含原來的a、b兩個向量，我們要求解得變量就從a、b變成A。

下面兩邊都乘以樣本矩陣X的轉置：

 

如果方程組有解，的行列式要大於0，即這個對稱矩陣是非奇異的。

下面我們讓右乘它的逆：

這樣等式右邊就變為EA，而等式左邊變為正規方程組的形式。

 

因為X、Y都是已知的，所以可以直接得到A這個包含a、b的系數矩陣，它也就是方程組的最小二乘解。

 

最小二乘法代碼實現：

# 最小二乘法
from numpy import  *
import matplotlib.pyplot as plt

# 加載數據集
def loadDataSet(fileName):
    X = []; Y = []
    fr = open(fileName, 'r')

    content = fr.read()
    fr.close()
    rowlist = content.splitlines()  # 按行轉換為一維表
    recordlist = array([row.split("\t") for row in rowlist if row.strip()])
    for line in recordlist:
        line = array(line)
        X.append(float(line[0]))
        Y.append(float(line[-1]))
    return X, Y


# 繪制圖形
def plotscatter(Xmat, Ymat, a, b, plt):
    fig = plt.figure()
    ax = fig.add_subplot(111)   # 繪制圖形位置
    ax.scatter(Xmat, Ymat, c='blue', marker='o')    # 繪制散點圖
    Xmat.sort() # 對Xmat各元素進行排序
    yhat = [a * float(xi) + b for xi in Xmat]   # 計算預測值
    plt.plot(Xmat, yhat, 'r')   # 繪制回歸線
    plt.show()
    return yhat


Xmat, Ymat = loadDataSet("/Users/FengZhen/Desktop/accumulate/機器學習/predict/線性回歸測試集.txt")    # 導入數據文件
meanX = mean(Xmat)  # 原始數據集的均值
meanY = mean(Ymat)
dX = Xmat - meanX   # 各元素與均值的差
dY = Ymat - meanY

sumXY = vdot(dX, dY)    # 返回兩個向量的點乘
sqX = sum(power(dX, 2)) # 向量的平方 (X-meanX)^2

# 計算斜率和截距
a = sumXY / sqX
b = meanY - a * meanX
print(a, b) # 1.199346340945075 -0.6940230973871095
plotscatter(Xmat, Ymat, a, b, plt)

 

正規方程組的代碼實現

#正規方程組
from numpy import  *
import matplotlib.pyplot as plt

# 加載數據集
def loadDataSet(fileName):
    X = []; Y = []
    fr = open(fileName, 'r')

    content = fr.read()
    fr.close()
    rowlist = content.splitlines()  # 按行轉換為一維表
    recordlist = array([row.split("\t") for row in rowlist if row.strip()])
    for line in recordlist:
        line = array(line)
        X.append(float(line[0]))
        Y.append(float(line[-1]))
    return X, Y

# 繪制圖形
def plotscatter(Xmat, Ymat, a, b, plt):
    fig = plt.figure()
    ax = fig.add_subplot(111)   # 繪制圖形位置
    ax.scatter(Xmat, Ymat, c='blue', marker='o')    # 繪制散點圖
    Xmat.sort() # 對Xmat各元素進行排序
    yhat = [a * float(xi) + b for xi in Xmat]   # 計算預測值
    plt.plot(Xmat, array(yhat), 'r')   # 繪制回歸線
    plt.show()
    return yhat


xArr, yArr = loadDataSet("/Users/FengZhen/Desktop/accumulate/機器學習/predict/線性回歸測試集.txt")
print(xArr, yArr)
m = len(xArr)
Xmat = mat(ones((m, 2)))  # 生成x坐標
for i in range(m):
    Xmat[i, 1] = xArr[i]
Ymat = mat(yArr).T  # 轉換為Y列
xTx = Xmat.T * Xmat # 矩陣左乘自身的轉置
ws = [] # 直線的斜率和截距
if linalg.det(xTx) != 0.0:  # 行列式不為0
    ws = xTx.I * (Xmat.T * Ymat)    # 矩陣正規方程組公式：inv(X.T * X) * X.T*Y
else:
    print("矩陣為奇異矩陣，無逆矩陣")
    sys.exit(0) # 退出程序
print(ws)
print(ws[1, 0], ws[0,0])    # 截距：斜率
plotscatter(xArr, yArr, ws[1, 0], ws[0,0], plt)