机器学习之线性回归

最近狂看机器学习，对数据挖掘算法甚是着迷。

“协同过滤”，“kNN”，“Bayes”，“k-means”，“决策树”这些都相对比较简单。涉及到距离和概率的计算。

推荐两本书

《数学之美》—— 纲领性书，激起对数学的兴趣。

《写给程序员的数据挖掘实践指南》—— 入门书，通俗易懂。

目前在啃《机器学习实战》，此书一大特色是狂用numpy，代码简洁优美，然而对数学公式解读不强。

看到第五章Logistic回归，一脸懵逼。

细细研究了一个礼拜，发现此章集微积分，线性代数，概率论为一体，相比之前章节，难度升了一个档次。

惭愧的是我还未全部解读通透。数学思维，强求不得，不急于一时，日后慢慢融会贯通。

参考了[置顶] 【机器学习笔记1】Logistic回归总结

地址如下：http://blog.csdn.net/dongtingzhizi/article/details/15962797 

本文举的例子是最简单的回归，线性回归，用的方式是梯度下降。尽可能详细讲解！

给出四个点，用直线拟合，结果如下：

原始数据坐标(4个点)：

x = [1, 3    , 5      , 7    ]
y = [5, 10.5, 16.8 , 23.1]

线性回归的思路：我们假设能用一条直线 y=kx+b 来串起所有的数据点。

那如何确定k和b的值？使得总误差最小！

假设k=3, b=2 .将4个x [1, 3, 5, 7 ]代入，计算得出四个y [5, 11, 17, 23]，目测已经很接近了。

那此时的误差怎么算？我们一般用 [5, 11, 17, 23] 和 原始数据 [5, 10.5, 16.8 , 23.1] 的方差之和来表示：

所以我们的目标就是要使得这个SSE最小，（变量就是k和b！）

插播三个矩阵：第一个θ是我们要求的目标矩阵，后两个X和Y是源数据矩阵。（注意上标和下标）

 

(注：在X的第二列全加1，是为了能够进行矩阵计算，四行二列 乘 二行一列 结果为 四行一列  (4*2)*(2*1)=(4*1)）

回到正题，下面介绍两个专业点的名词，h函数和cost函数。

h函数就是拟合函数：

cost函数就是误差函数，就是前面提到的SSE，记为J(θ)函数，加个1/2系数是为了之后求导能消掉：（m等于几？m=4啦，因为有4个点。(i)就表示是第几个，不是次方。）

我们的目标是什么？让J(θ)函数最小！！！

注意，这里的变量不是xy了，xy都是已知的，这里的变量是k和b，也就是θ！

我们用梯度下降法去求。

好，最深奥的问题出现了，梯度是什么，梯度下降又是什么，梯度怎么求。。。（概念请自行百度，之后会给图）

梯度的算法，就是对θ的每个分量求偏导，合起来的向量就是梯度。

下面对θ的j分量求偏导。（θ一共就2个分量，k和b，所以j=1或2）

如果求导不知道的，那温习一下微积分。。。

最后一步，是因为对θj求导，所以除了θj的系数有用，其他都是常数，常数的导数是0.

以上是求梯度的公式，下面是梯度下降的公式：

表示往θj的方向下降α.

把符号变正，用e来替换y-h，如下：

这只是一个分量上的公式，那把所有分量（这里只有两个分量，j=1或2）展开得：

以下是梯度下降核心代码（附注释）：

weight = ones((2,1)) #初始化θ
alpha = 0.001 #设置α，步长


for k in range(100): #走100步
    h = x*weight #计算h
    err = y-h #计算E
    weight = weight + alpha * x.T * err #迭代计算θ

print weight

（最后会给出所有代码，稍安勿躁。）

公式和理论是很不形象的，接下来上梯度图。

重申一遍：梯度是求SSE最小值时计算的，SSE由向量θ决定，也就是k和b。

那我们看一下，对于不同的k和b，SSE是多少。下面绘制三维曲面图（我只选前20组数据。）

for k in range(100):
    h = x*weight
    err = y-h
    SSE = 0.5*(err.getA()**2).sum()
    print weight[0,0],weight[1,0],SSE
    weight = weight + alpha * x.T * err
print weight

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

a = [1,1.1822,1.3485288,1.50036849,1.638980804,
     1.765517606,1.88103045,1.986479323,2.082740612,2.17061439,
     2.250831056,2.324057406,2.390902171,2.451921072	,2.507621438,
     2.558466417,2.604878829,2.647244671,2.685916329,2.721215502]

b = [1,1.0354,1.0677432,1.097295766,1.124300688,
     1.148979792,1.171535591,1.192152961,1.21100068,1.228232828,
     1.243990066,1.258400809,1.271582287,1.283641524,1.29467622,
     1.304775572,1.314021007,1.322486862,1.330241,1.337345375]

c = [197.95,164.99997168,137.538726568,114.651946034,95.5776181715,
     79.6806666225,66.4318057178,55.3899179207,46.1873668246,38.5177567047,
     32.1257310768,26.7984706073,22.3586072967,18.6583190148,15.5744077648,
     13.0041978085,10.8621170786,9.07684805841,7.5889532661,6.34889628462]

xdata = [1,3,5,7]
ydata = [5,10.5,16.8,23.1]

x1 = []
for xd in xdata:
    x1.append([xd,1])

x = mat(x1)
y = mat(ydata).T

aall = []
ball = []
sall = []
for a1 in a:
    for b1 in b:
        weight = array([[a1],[b1]])
        h = x * weight
        err = y - h
        SSE = SSE = 0.5*(err.getA()**2).sum()
        aall.append(a1)
        ball.append(b1)
        sall.append(SSE)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.view_init(elev=20, azim=70)
ax.scatter(a,b,c)
ax.plot(a,b,c)
ax.plot_trisurf(aall,ball,sall,color='black')
plt.show()

蓝色的点和线就是“机器”根据梯度下降法选择走的路径，从上到下（每次都选择最陡的方向前进）。黑色为曲面。

经过20轮，SSE从200成功降到6.之后下降的趋势越来越缓。

正文完毕

PS1， 本文说的是线性回归，是最简单的一种。但是一元二次回归用梯度也是这么实现。

梯度下降的核心公式便是：

所以，这种回归需要的是事先给出多项式，然后去通过梯度下降求系数。不管多项式是一次而是二次。

另外，重要的是Cost函数，即判断误差大小的函数。平方差和只是一种，梯度都是基于这个函数去求导的。

PS2，稍微说一下书上的Logistic回归。他的作用是根据标签画一条直线把点分为两类。会把值再代入Sigmoid函数，所以理解起来和求导都比较困难。

还有，Logistic回归中的Cost函数和最大似然估计有关。

然而奇妙的是，经过求导，最大似然梯度和最简单的线性回归的梯度公式是一样的。这就匪夷所思了，是选择Cost函数时候故意这样的么？

PS3，代码。

from numpy import *
import matplotlib
import matplotlib.pyplot as plt

xdata = [1, 3    , 5      , 7    ]
ydata = [5, 10.5, 16.8 , 23.1]

x1 = []
for xd in xdata:
    x1.append([xd,1])

x = mat(x1)
y = mat(ydata).T
weight = ones((2,1))
alpha = 0.001

for k in range(100):
    h = x*weight
    err = y-h
    SSE = 0.5*(err.getA()**2).sum()
##    print weight[0,0],weight[1,0],SSE
    weight = weight + alpha * x.T * err
print weight

weightArray = weight.getA()

n = len(xdata)
xcord = []
ycord = []
for i in range(n):
    xcord.append(xdata[i])
    ycord.append(ydata[i])

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord, ycord, s=30, c='red', marker='s')
a = arange(0.0,8.0,0.1)
b = weightArray[0] * a + weightArray[1]
ax.plot(a,b)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()