机器学习系列-Locally weighted linear regression(2)

机器学习
这是记录自学的过程，目前的理论基础就是：大学高等数学+线性代数+概率论。编程基础：C/C++,python
在观看机器学习实战这本书，慢慢介入。相信有读过以上三门课的人完全可以开始自学机器学习了，当然我上面这三门课学的一般，所以你只知道有这么一个公式或名词，不懂可以百度之深究之。在写这篇文章的时候作者机器学习还没学完，故文章中的错误还请不吝指出。再次声明，系列文章只是分享学习过程，学习点滴，不能保证文章的技术含量。后续再技术的不断完善中，我会重新再捋一遍这个学习过程，纠正其中错误。

目前的学习方法如下：

理论：斯坦福Andrew NG machine learning，系列课程。

代码实战：引用书籍 机器学习实战（python）。

真枪实战： 研究一个案例，建模选择合适的机器学习算法，分析数据。
演变：基于C/C++ 在linux下实现部分机器学习过程。将在学习中期阶段去实现。

一、Locally weighted linear regression

              局部加权线性回归(LWLR)的提出无疑是为了改进线性回归的一些不可避免的缺点，比如欠拟合。Andrew ng课程通过类比 1）y=k1+k2x   2)y=k1+k2x+k3x^2....的图像效果拟合程度，以此推断当我们加入更多的features的时候，结果会更加接近真实结果。当然features达到一定数量时，就会出现过拟合的现象。
       简单回顾一下线性回归的算法：（讲义截图）



LWLR在此基础上乘上一个权值w



在一个数据集中，如果采用线性回归的算法，数据之间的影响是比较大的，此算法对所有数据一视同仁，并不会很理智的去筛选出有用的数据。二LWLR加入了权值概念，为每一个数据定义一个权值系数。如果认为这个数据有用，那这个数据对应的权值相应就比较大；相反如果认为这个数据偏差太大，其对应的权值系数就会很小。这样一来，让有用的数据影响性增大，偏差太大的数据影响性减小，最后会使整个系统更加准确的预计出y值。

能实现以上描述现象w的数学模型，这里给出一个标准的选择高斯核：



这是一个很像标准正态分布的模型，我们知道正太分布的图像是这样：



τ的作用是控制随样本点与待测值的距离的递增，权值w成指数级衰减的速率，x1，x2是我们要估计的数据。可以看到当 |x(i) − x|越小，说明数据离我们的输入数据越近（简单说就是有利数据的推测，我们需要留下它），此时可以看出w的值逼近1；相反 |x(i) − x|越大，说明数据距我们的输入数据比较远（简单说就是噪声数据，我们需要尽量忽略它），此时可以看出w的值逼近0. 

     上图中的浅蓝色点就是可能的噪声数据，它们对应的权值很小，无限接近0，当这些数据乘以它们的权值后，它们对X1，X2的影响就大大降低了。所谓的局部加权就是上述的整个过程，局部就是指输入数据的一定范围内的数据是对该输入数据有影响的。

理论推导同线性回归差不多再其基础之上乘以一个权值w。

先看看线性回归的理论推导结果：



因为后面会用python来检测这个算法，所以我们统一用矩阵公式，LWLR在其基础上乘以W得出来的的结果是：



这里省略了，LWLR的结论推导过程，感兴趣的同学可以去了解下。虽然我可能理解不够透彻，无法很顺畅的把推导过程写出来，这是件很遗憾的事，我相信后面慢慢会理解这些的。

二、用Python实现的过程

    我们先来看下我们要使用的数据的结构：



........

前面两列代表X1和X2，会被读到一个矩阵里面xArr。第三列代表Y，被读到yArr.第一列是作为常数用的，所以数值都是1。

 先看一下对单点进行回归的主要代码：

def lwlr(testPoint,xArr,yArr,k=1.0):
xMat = mat(xArr); yMat = mat(yArr).T
m = shape(xMat)[0]
weights = mat(eye((m)))
for j in range(m):                      #next 2 lines create weights matrix
diffMat = testPoint - xMat[j,:]     #
weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))
xTx = xMat.T * (weights * xMat)
if linalg.det(xTx) == 0.0:
print "This matrix is singular, cannot do inverse"
return
ws = xTx.I * (xMat.T * (weights * yMat))
return testPoint * ws

深入分析：

这段代码主要作用是，任意给定空间一点，计算出对应的预测值yHat。

xMat = mat(xArr); yMat = mat(yArr).T   把数组变成举证此时xMat是n行两列，yMat是n行一列。

m = shape(xMat)[0]      读取xMat矩阵的行数，如果后面跟着[1]则是列数，shape保存着这样一组数据。

weights = mat(eye((m)))  创造一个m*m的单位矩阵，也即是对角权值矩阵

for j in range(m):                      #next 2 lines create weights matrix

        diffMat = testPoint - xMat[j,:]     #

        weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))

上面是一个for循环，为了创造w这个矩阵，也就是上面所说的高斯核。先求出被测数据与已知数据的差值，因为都是1*2的矩阵，所以可以直接相减。接着diffMat*diffMat.T其实就是求|diffMat|的过程，求出来是一个值。套入高斯公式，就可以循环求出高斯核的每一个数据对应的权值。

xTx = xMat.T * (weights * xMat)    xMat转置后是一个2*n的矩阵；eights*xMat 就是n*n 与n*2相乘，结果是n*2的矩阵。最后两个相乘，变成2*2的矩阵。

if linalg.det(xTx) == 0.0:

        print "This matrix is singular, cannot do inverse"

        return

判断刚才生成的2*2矩阵是否存在逆矩阵，也就是求他的行列值是否为零。

ws = xTx.I * (xMat.T * (weights * yMat))   这就是求回归系数的公式套用，xTx.I 就是求逆矩阵的意思，后面也就是如上面所说。

最后返回预测的的值，testPoint*ws。

以上是对一个值进行预测，如果多个值可以调用下面这个方法：

def lwlrTest(testArr,xArr,yArr,k=1.0): #loops over all the data points and applies lwlr to each one
m = shape(testArr)[0]
yHat = zeros(m)
for i in range(m):
yHat[i] = lwlr(testArr[i],xArr,yArr,k)
return yHat

上述也就是输入变成了一个数组，然后循环调用刚才那个方法，最后得到yHat的矩阵集。
从上面也可以看出，这个算法对每个点进行预测的时候，都必须使用整个数据集，虽然有些数据集会被判定为可忽略数据即权值很低的数据，但还是会遍历计算一遍。

看一下上米哦按在K=1.0时计算后得出来的数据的图：



可以看出，K=1时导致大部分数据的权值都比较大，如同将所有数据视为等同权重，得出的最佳拟合直线和标准的回归一致，为一条直线（蓝色）。所以需要微调一下K的值，比如调成0.003的结果：



可以看出，这个时候纳入了太多的噪声数据，处于过拟合状态。再把k调成0.01时：



这个时候数据呈现出来的曲线比较很好的拟合了大多数数据点.

下面是画图的代码，原本没有提供函数形式的，手输入比较麻烦，整理成函数形式比较方便：

def plotshow(xMat,yHat,yArr):
srtInd = xMat[:,1].argsort(0)
xSort = xMat[srtInd][:,0,:]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xSort[:,1],yHat[srtInd])
ax.scatter(xMat[:,1].flatten().A[0],mat(yArr).T.flatten().A[0],s=2,c='red')
plt.show()