直线、抛物线、猜猜看！

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

• 直线抛物线猜猜看
  • 一切从凑数开始
  • 现代人猜是直线
  • 分明是抛物线
  • 猜猜看

一切从凑数开始

机器学习的目的通常有两类：

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-8
import numpy as np
import scipy.optimize as opt 
from pylab import *
#x和y数据
x=np.array([-, -, - , -, -,-, -, -, -, -,-, -, -, -, -,-, -, -, -, -,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ])
y=np.array([ ,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ,,  ,   ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ])
#猜测的模型
def f(k,b,x):
    return k*x+b;
#优化目标函数
def e(p):
    k=p[]
    b=p[]
    return sum((f(k,b,x)-y)**)/;
k=
b=
p=[k,b]
p=opt.fmin(e,p) #求解k和b
k=p[]
b=p[]
figure()
plot(x,y,"r+")
#绘制拟合的直线
nx=np.linspace(-,,)
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-8
import numpy as np
import scipy.optimize as opt 
from pylab import *
#x和y数据
x=np.array([-, -, - , -, -,-, -, -, -, -,-, -, -, -, -,-, -, -, -, -,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,          ])
y=np.array([ ,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ,,  ,   ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ])
#猜测的模型
def f(k,b,x):
    return k*x+b;
#优化目标函数
def e(p):
    k=p[]
    b=p[]
    return sum((f(k,b,x)-y)**)/;
#目标函数梯度
def grade(p):
    k=p[]
    b=p[]
    gk=sum((f(k,b,x)-y)*x)
    gb=sum(f(k,b,x)-y)
    return np.asarray((gk,gb))
k=
b=
p=[k,b]
p=opt.fmin_cg(e,p,fprime=grade) #Newton-CG法 
k=p[]
b=p[]
figure()
plot(x,y,"r+")
nx=np.linspace(-,,)
ny=f(k,b,nx)
plot(nx,ny)
show()
           

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

以上为fmin的迭代次数。

以上为fmin_cg的迭代次数。

分明是抛物线！

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

f(x)=ax2+bx+c

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

# coding=utf-8
import numpy as np
import scipy.optimize as opt 
from pylab import *
#x和y数据
x=np.array([-, -, - , -, -,-, -, -, -, -,-, -, -, -, -,-, -, -, -, -,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ])
y=np.array([ ,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ,,  ,   ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ])
#猜测的模型
def f(a,b,c,x):
    return a*x**+b*x+c
#优化目标函数，
def e(p):
    a=p[]
    b=p[]
    c=p[]
    return sum((f(a,b,c,x)-y)**)/;
a=
b=
c=
p=[a,b,c]
p=opt.fmin(e,p) #求解k和b
a=p[]
b=p[]
c=p[]
figure()
plot(x,y,"r+")
#绘制拟合的抛物线
nx=np.linspace(-,,)
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-8
from keras.models import Sequential
from keras.layers import Dense
import numpy as np
from pylab import *
#x和y数据
x=np.array([-, -, - , -, -,-, -, -, -, -,-, -, -, -, -,-, -, -, -, -,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ])
y=np.array([ ,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,   ,  ,  ,,  ,   ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ,,  ,  ,  ,  ])
trainIdx=arange(,,)
testIdx=arange(,,)
model = Sequential()
model.add(Dense(, activation='tanh',input_shape=(,)))
model.add(Dense(, activation='tanh'))
model.add(Dense(, activation='tanh'))
model.add(Dense(, activation='tanh'))
model.compile(loss='mean_squared_error', optimizer="adam")
model.fit(x[trainIdx], y[trainIdx], batch_size=,epochs=,shuffle=True,verbose=,validation_data=(x[testIdx],y[testIdx]))
print(model.predict(x))
figure()
plot(x,y,"r+")
#绘制拟合的抛物线
nx=np.linspace(-,,)
ny=model.predict(nx)
plot(nx,ny)
show()
           

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

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

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