天天看点

直线、抛物线、猜猜看!直线、抛物线、猜猜看!

直线、抛物线、猜猜看!

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

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

一切从凑数开始

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

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个节点。其运行结果如下:

直线、抛物线、猜猜看!直线、抛物线、猜猜看!

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

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

继续阅读