超定方程组的最小二乘解
超定线性方程组：方程数多于变量数，而且这些方程组往往是不相容的。例如给定一个m * n的方程组 Ax= b，其中m > n。我们要做的就是寻找一个向量x，使得Ax最接近 b 。我们可以构造一个方程：
![][equtation]
[equtation]: http://latex.codecogs.com/svg.latex?r(\bm{x})=\bm{b}-A\bm{x}
b和Ax之间的距离为：
![][equtation2]
[equtation2]: http://latex.codecogs.com/svg.latex?\left|\bm{b}-A\bm{x}\right|=\left|r(\bm{x})\right|
我们希望寻找一个向量x，使得||r(x)||^2最小的，最后得到的x就成为最小二乘解。
根据正交性（证明省略），我们可知得到最小二乘解的矩阵形式是：
![][equtation3]
[equtation3]: http://latex.codecogs.com/svg.latex?\hat{x}=(A{T}A){-1}A^{T}\bm{b}
上式中的\hat{x}就是最优解。通过上式我们可以看到，上述公式中包含(A^T* A)(-1)，所以在回归之前，先要保证AT * A可逆，也就是它的行列式不等于0，（可以理解为0没有倒数）。
TIC
首先构造一个读取数据的函数，这个函数在MLA这本书里很常见。
from numpy import *
def loadDataSet(fileName):#文件读取函数
numFeat = len(open(fileName).readline().split('\t')) - 1#读出每行有n个数据，将n-1作为特征数
dataMat = []; labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr =[]
curLine = line.strip().split('\t')
for i in range(numFeat):
lineArr.append(float(curLine[i]))#将n-1个特征组成一个新list
dataMat.append(lineArr)#形成双层list形成一个矩阵
labelMat.append(float(curLine[-1]))#将最后一个数值作为目标值
return dataMat,labelMat
'''
----------------------------loadDataSet()说明---------------------------------
txt里的数据预览：
每行有n个数据，前n个为特征，最后一个为目标值（默认）
对应回归方程 A x = y
_____________________ ________
n-1个特征(A) | 目标值(y)
_____________________ ________
a1 a2 a3 ... a(n-1) | an
b1 b2 b3 ... b(n-1) | bn
c1 c2 c3 ... c(n-1) | bn
...
g1 g2 g3 ... g(n-1) | gn
_____________________ ________
要强调的是：读取的txt里数据是用tab分割的，用其他符号分隔须将函数加以修改。
'''
拟合函数，结合公式很容易理解：
def standRegres(xArr,yArr):
xMat = mat(xArr); yMat = mat(yArr).T
xTx = xMat.T*xMat
if linalg.det(xTx) == 0.0:
print "This matrix is singular, cannot do inverse"
return
ws = xTx.I * (xMat.T*yMat)
return ws
代码写完了，我们跑一个试下：
不知道哪来的一组数据
if __name__ == '__main__':
xArr , yArr = loadDataSet('ex0.txt')
ws = standRegres(xArr, yArr)
print ws
结果是：
[[ 3.00774324]
[ 1.69532264]]
分析
首先，我们预览一下ex0.txt：
1.000000 0.067732 3.176513
1.000000 0.427810 3.816464
...
1.000000 0.995731 4.550095
我们发现该数据有两个特征，其中第一个都是1，为什么？实际上这个数据是经过预处理的，其实这个数据只有一个特征一个目标值，在每行前面加上1相当于个Ax = y中的A加上了一个常数项使之偏置。我们算出来的回归系数[3.00774324，1.69532264]分别对应y = a + bx 中的a,b。
也就是说，如果每个样本有n个特征，它对应的向量就应该是n+1维。
如果不做这个预处理，那么我们相当于使拟合函数强制过原点。
绘图
if __name__ == '__main__':
xArr , yArr = loadDataSet('ex0.txt')
ws = standRegres(xArr, yArr)
print ws
xMat = mat(xArr)
yMat = mat(yArr)
yHat = xMat*ws
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])
xCopy = xMat.copy()
xCopy.sort(0)
yHat=xCopy*ws
ax.plot(xCopy[:,1],yHat)
plt.show()
得到这样一个图片，忍不住想吐槽真的太丑了：
对于所有的二维数据（一个特征一个目标值），我们都可以这样拟合直线，但实际上我们知道并不是所有数据都符合线性关系，常用相关系数来表明我们的直线和数据是否吻合，Numpy中就自带计算方法。
yHat = xMat*ws
>>> corrcoef(yHat.T, yMat)
>>> array([[1. , 0.98647356],
[0.987356, 1. ]])
对于高维数据，只是矩阵的阶在增加，算法完全相同。
网友评论