直线、抛物线、猜猜看！

本文版权属于重庆大学计算机学院刘骥，禁止转载

一切从凑数开始

机器学习的目的通常有两类：
1. 模型拟合
2. 分类

这篇文章主要谈谈模型拟合问题。首先来看看什么叫做模型？所谓模型就是：

f(X)

完了。完了？真的完了。我既没有骗你，也没有逗你。模型肯定不是车模，也不是航模或者船模。难道模型不应该是复杂复杂的东西吗？比如什么宇宙模型？给你看一个描述宇宙的基本规律的模型：

F(m1,m2,r)=Gm1m2r2

牛顿的万有引力方程！怎么样高级吧？再来看个更高级的：

E(m)=mc2

爱因斯坦质能方程！怎么样还觉得low吗？物理学家用函数来描述宇宙万物。火箭、卫星、手机、电脑都是依据这些基本模型构造出来的。那么这些模型怎么来的？两种方式：
1. 推导，比如E=mc2就是推导出来的。
2. 拟合，比如通过小车下滑的物理实验，得出v=v0+at。

所谓拟合，就是凑数。这是古代科学家的独门秘籍。托勒密的天体运行模型，就是依据观察数据凑数而来。用托勒密的模型预测行星运行的轨迹竟然能够做到差不太多（但与牛顿力学比差了十万八千里，爱因斯坦又比牛顿提高了头发丝大小的精度）。

今年是2017年，是一个计算机已经被发明了几十年，智能手机极度普及的年份。如果现代人穿越到托勒密的时代，现代人会怎么解决天体运行的问题？

现代人猜是直线

现代人在托勒密时代看到了下面这样的数据：

如果将这个数据画出来，效果就是：

现代人觉得x和y之间的关系是什么？现代人猜是一条直线！像这样的直线：

y=kx+b

是不是有一口鲜血要吐出来的感觉？既然现代人说这是一条直线，我们认了吧。问题是这条直线长什么样？现代人说，两点确定一条直线，任意选两组x和y，就可以确定这条直线。Oh，my god！这里面40个点，任意选择两个那得有多少条直线呢？莫慌，现代人说，可能有些点是噪声，我们要找一条直线，使得这条直线是一条“平均直线”。具体来说，这条直线必然满足如下条件：

argmink,b∑40i=1(yi−kxi−b)240‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√

上述公式称为Root Mean Squared Error(RMSE)，如果不开根号就是 Mean Squared Error(MSE)：

argmink,b∑40i=1(yi−kxi−b)240

更一般的情况下：

RMSE=∑mi=1(yi−f(xi))2m‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√

MSE=∑mi=1(yi−f(xi))2m

换言之，现代人的目标是求RMSE或者MSE关于k和b的最优解。我们发现RMSE和MSE的最优解是相同的，为了计算方便，通常优化的目标简化为：

E=∑i=1m(yi−f(xi))22

这样做不会改变最优解，却可以简化计算（现代人是一个数学家）。现在我们来看看用MSE公式如何求最优解。现代社会是2017年，现代人是穿越回去的，所以他带了一台安装了python程序的电脑。因此问题变得so easy：

# coding=utf-8import numpy as npimport scipy.optimize as opt from pylab import *#x和y数据x=np.array([-1., -0.94871795, -0.8974359 , -0.84615385, -0.79487179,-0.74358974, -0.69230769, -0.64102564, -0.58974359, -0.53846154,-0.48717949, -0.43589744, -0.38461538, -0.33333333, -0.28205128,-0.23076923, -0.17948718, -0.12820513, -0.07692308, -0.02564103,0.02564103,  0.07692308,  0.12820513,  0.17948718,  0.23076923,0.28205128,  0.33333333,  0.38461538,  0.43589744,  0.48717949,0.53846154,  0.58974359,  0.64102564,  0.69230769,  0.74358974,0.79487179,  0.84615385,  0.8974359 ,  0.94871795,  1.])y=np.array([ 0.60653066,  0.63760719,  0.66851557,  0.69908135,  0.72912464,0.75846179,  0.78690719,  0.81427516,  0.84038198,  0.86504789,0.88809911,  0.90936995,  0.92870467,  0.94595947,  0.96100424,0.97372416,  0.98402121,  0.9918154 ,  0.99704579,  0.99967132,0.99967132,  0.99704579,  0.9918154 ,  0.98402121,  0.97372416,0.96100424,  0.94595947,  0.92870467,  0.90936995,  0.88809911,0.86504789,  0.84038198,  0.81427516,  0.78690719,  0.75846179,0.72912464,  0.69908135,  0.66851557,  0.63760719,  0.60653066])#猜测的模型def f(k,b,x):    return k*x+b;#优化目标函数def e(p):    k=p[0]    b=p[1]    return sum((f(k,b,x)-y)**2)/2;k=10b=10p=[k,b]p=opt.fmin(e,p) #求解k和bk=p[0]b=p[1]figure()plot(x,y,"r+")#绘制拟合的直线nx=np.linspace(-1,1,40)ny=f(k,b,nx)plot(nx,ny)show()

最终的执行结果如下：

蓝色的直线穿过所有点的中部，的确能够使MSE最小。fmin采用downhill simplex方法，如果我们采用fmin_cg算法，就可以体现出MSE公式的便利。使用fmin_cg需要对MSE求k和b的偏导数，结果如下：

∂E∂k=∑i=1m2(yi−f(xi))2(−∂f(xi)∂k)=∑i=1m(yi−f(xi))(−xi)=∑i=1m(f(xi)−yi)xi

∂E∂b=∑i=1m2(yi−f(xi))2(−∂f(xi)∂b)=∑i=1m(yi−f(xi))(−1)=∑i=1m(f(xi)−yi)

于是程序可以改写为：

# coding=utf-8import numpy as npimport scipy.optimize as opt from pylab import *#x和y数据x=np.array([-1., -0.94871795, -0.8974359 , -0.84615385, -0.79487179,-0.74358974, -0.69230769, -0.64102564, -0.58974359, -0.53846154,-0.48717949, -0.43589744, -0.38461538, -0.33333333, -0.28205128,-0.23076923, -0.17948718, -0.12820513, -0.07692308, -0.02564103,0.02564103,  0.07692308,  0.12820513,  0.17948718,  0.23076923,0.28205128,  0.33333333,  0.38461538,  0.43589744,  0.48717949,0.53846154,  0.58974359,  0.64102564,  0.69230769,  0.74358974,0.79487179,  0.84615385,  0.8974359 ,  0.94871795,  1.        ])y=np.array([ 0.60653066,  0.63760719,  0.66851557,  0.69908135,  0.72912464,0.75846179,  0.78690719,  0.81427516,  0.84038198,  0.86504789,0.88809911,  0.90936995,  0.92870467,  0.94595947,  0.96100424,0.97372416,  0.98402121,  0.9918154 ,  0.99704579,  0.99967132,0.99967132,  0.99704579,  0.9918154 ,  0.98402121,  0.97372416,0.96100424,  0.94595947,  0.92870467,  0.90936995,  0.88809911,0.86504789,  0.84038198,  0.81427516,  0.78690719,  0.75846179,0.72912464,  0.69908135,  0.66851557,  0.63760719,  0.60653066])#猜测的模型def f(k,b,x):    return k*x+b;#优化目标函数def e(p):    k=p[0]    b=p[1]    return sum((f(k,b,x)-y)**2)/2;#目标函数梯度def grade(p):    k=p[0]    b=p[1]    gk=sum((f(k,b,x)-y)*x)    gb=sum(f(k,b,x)-y)    return np.asarray((gk,gb))k=10b=10p=[k,b]p=opt.fmin_cg(e,p,fprime=grade) #Newton-CG法 k=p[0]b=p[1]figure()plot(x,y,"r+")nx=np.linspace(-1,1,40)ny=f(k,b,nx)plot(nx,ny)show()

两段代码的运行结果相同，但由于求解目标函数的方法不同，因此速度上会有一些差异。

以上为fmin的迭代次数。

以上为fmin_cg的迭代次数。

分明是抛物线！

看起来那个数据的确不能用直线来描述。现代人知道除了直线，抛物线（二次曲线）也是常见的模型，那么试试看吧。抛物线的函数如下：

f(x)=ax2+bx+c

由于只是修改了模型，因此其他部分都不变，程序代码修改如下：

# coding=utf-8import numpy as npimport scipy.optimize as opt from pylab import *#x和y数据x=np.array([-1., -0.94871795, -0.8974359 , -0.84615385, -0.79487179,-0.74358974, -0.69230769, -0.64102564, -0.58974359, -0.53846154,-0.48717949, -0.43589744, -0.38461538, -0.33333333, -0.28205128,-0.23076923, -0.17948718, -0.12820513, -0.07692308, -0.02564103,0.02564103,  0.07692308,  0.12820513,  0.17948718,  0.23076923,0.28205128,  0.33333333,  0.38461538,  0.43589744,  0.48717949,0.53846154,  0.58974359,  0.64102564,  0.69230769,  0.74358974,0.79487179,  0.84615385,  0.8974359 ,  0.94871795,  1.])y=np.array([ 0.60653066,  0.63760719,  0.66851557,  0.69908135,  0.72912464,0.75846179,  0.78690719,  0.81427516,  0.84038198,  0.86504789,0.88809911,  0.90936995,  0.92870467,  0.94595947,  0.96100424,0.97372416,  0.98402121,  0.9918154 ,  0.99704579,  0.99967132,0.99967132,  0.99704579,  0.9918154 ,  0.98402121,  0.97372416,0.96100424,  0.94595947,  0.92870467,  0.90936995,  0.88809911,0.86504789,  0.84038198,  0.81427516,  0.78690719,  0.75846179,0.72912464,  0.69908135,  0.66851557,  0.63760719,  0.60653066])#猜测的模型def f(a,b,c,x):    return a*x**2+b*x+c#优化目标函数，def e(p):    a=p[0]    b=p[1]    c=p[2]    return sum((f(a,b,c,x)-y)**2)/2;a=10b=10c=10p=[a,b,c]p=opt.fmin(e,p) #求解k和ba=p[0]b=p[1]c=p[2]figure()plot(x,y,"r+")#绘制拟合的抛物线nx=np.linspace(-1,1,40)ny=f(a,b,c,x)plot(nx,ny)show()

运行结果如下：

结果让人感动，几乎完美的贴合。上述程序，相对于之前的程序，重新定义了待拟合的模型，其余部分保持不变。那么我们可以得到什么结论呢？

猜猜看

机器学习的目的之一是拟合，而拟合的前提是猜测。对，你没有看错，的确是猜测。拟合的过程都是差不多的，人所需要做的事情就是猜测哪种模型更符合数据的分布。也是古代学者所采用的——凑数。与古代不同，现代有了python这样的高级工具，如果我们能够猜对模型，那么就有能力把它拟合出来。

下面介绍一下高大上的太阳系标准模型的拟合工作。与前面的程序类似，天文学家已经有了大量的观测数据，什么叫做大量呢？几个GB？NO！几个TB？NO！几个PB？NO！每天有几个PB！于是天文学家给出了一个太阳系标准模型，也就是一系列的函数。利用观测数据，我们可以对模型进行拟合。拟合完成之后，利用这个模型进行模拟，再与实际的太阳系的情况进行对比。如果模拟的情况与太阳系的情况吻合，那么这个模型就是正确的。反之，这个模型可能是有问题的。水星进动问题最初被认为是牛顿模型的误差，但爱因斯坦发明相对论之后，就用更好的模型进行了修正。

再来看看，我们前面拟合的结果。二次函数模型并不能完美贴合曲线。所以我们可以进一步修正模型。如果上帝告诉我们，这个数据其实是由函数f(x)=e−x22生成的，那么我们就可以得到一个完美的模型。物理学家一直以来就在做这个工作：妄图猜准宇宙的正确模型。

现在你应该发现：最优化的算法是已经写好的，拟合的过程是标准的，模型是物理学家提出来的，好像没有程序员什么事儿！事实也的确如此。程序员的工作仅仅是写出了python的函数库而已。如果真要程序员发光发热，也只能是如何处理每天几个PB的天文数据。程序员并不擅长提出模型，而这恰好是其他专业的优势所在（比如数学、物理、统计学、经济学）。为什么不擅长呢？老师没教！

有个好消息，程序员发明了深度神经网络！既然程序员不擅长提出模型，那么干脆把提出模型这种事情也交给计算机！

# coding=utf-8from keras.models import Sequentialfrom keras.layers import Denseimport numpy as npfrom pylab import *#x和y数据x=np.array([-1., -0.94871795, -0.8974359 , -0.84615385, -0.79487179,-0.74358974, -0.69230769, -0.64102564, -0.58974359, -0.53846154,-0.48717949, -0.43589744, -0.38461538, -0.33333333, -0.28205128,-0.23076923, -0.17948718, -0.12820513, -0.07692308, -0.02564103,0.02564103,  0.07692308,  0.12820513,  0.17948718,  0.23076923,0.28205128,  0.33333333,  0.38461538,  0.43589744,  0.48717949,0.53846154,  0.58974359,  0.64102564,  0.69230769,  0.74358974,0.79487179,  0.84615385,  0.8974359 ,  0.94871795,  1.])y=np.array([ 0.60653066,  0.63760719,  0.66851557,  0.69908135,  0.72912464,0.75846179,  0.78690719,  0.81427516,  0.84038198,  0.86504789,0.88809911,  0.90936995,  0.92870467,  0.94595947,  0.96100424,0.97372416,  0.98402121,  0.9918154 ,  0.99704579,  0.99967132,0.99967132,  0.99704579,  0.9918154 ,  0.98402121,  0.97372416,0.96100424,  0.94595947,  0.92870467,  0.90936995,  0.88809911,0.86504789,  0.84038198,  0.81427516,  0.78690719,  0.75846179,0.72912464,  0.69908135,  0.66851557,  0.63760719,  0.60653066])trainIdx=arange(0,40,2)testIdx=arange(1,40,2)model = Sequential()model.add(Dense(1, activation='tanh',input_shape=(1,)))model.add(Dense(5, activation='tanh'))model.add(Dense(5, activation='tanh'))model.add(Dense(1, activation='tanh'))model.compile(loss='mean_squared_error', optimizer="adam")model.fit(x[trainIdx], y[trainIdx], batch_size=20,epochs=4000,shuffle=True,verbose=1,validation_data=(x[testIdx],y[testIdx]))print(model.predict(x))figure()plot(x,y,"r+")#绘制拟合的抛物线nx=np.linspace(-1,1,40)ny=model.predict(nx)plot(nx,ny)show()

这是一段用深度神经网络拟合的程序，使用Keras框架。神经网络有两个隐藏层，每个隐藏层有各有5个节点。其运行结果如下：

怎么样？不算太差吧。但深度神经网络的使用有一个大前提，那就是数据量足够多，否则计算机不能推导出模型。初次之外，记得没次运行程序之前，进行虔诚的祈祷。

如果你是数学白痴，又有大量的数据，那就用深度神经网络，让计算机猜猜看吧！