直線、抛物線、猜猜看！

本文版權屬于重慶大學計算機學院劉骥，禁止轉載

• 直線抛物線猜猜看
  • 一切從湊數開始
  • 現代人猜是直線
  • 分明是抛物線
  • 猜猜看

一切從湊數開始

機器學習的目的通常有兩類：

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個節點。其運作結果如下：

怎麼樣？不算太差吧。但深度神經網絡的使用有一個大前提，那就是資料量足夠多，否則計算機不能推導出模型。初次之外，記得沒次運作程式之前，進行虔誠的祈禱。

如果你是數學白癡，又有大量的資料，那就用深度神經網絡，讓計算機猜猜看吧！