机器学习-预测-线性系统的预测（最小二乘法、正规方程式实现）

机器学习-预测-线性系统的预测

现在预测学的核心概念：回归。从数学的角度，为事物（系统）的预测提供现代的技术方法。

 

回归与现代预测学

统计学上最初回归的含义由高尔顿（达尔文的表弟）通过研究父母身高与孩子身高得出。

矮个父母所生的儿子往往会比其父母更高，高个父母所生儿子的身高却回降到多数人的平均身高。也就是说，当父母身高走向极端时，子女的身高不会进一步极端化，他们的身高要比父母更接近平均身高，即有”回归“到平均数的趋势，这就是统计学上最初回归的含义。

高尔顿把这一现象叫作”向平均数方向的回归“(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)