最小二乘法（2）——多項式函數可以擬合非線性問題原理

　　一個複雜的多項式能夠「過擬合」任意數據，言外之意是多項式函數能夠接近於任何函數，這是什麼道理呢？

泰勒公式

　　欲理解多項式函數的過擬合，必先理解泰勒公式。

　　泰勒公式是一種計算近似值的方法，它是一個用函數某點的信息描述在該點附近取值的公式。已知函數在某一點的各階導數值的狀況之下，泰勒公式能夠用這些導數值作係數構建一個多項式來逼近函數在這一點的鄰域中的值。

　　若是f(x)在x₀處具備任意階導數，那麼泰勒公式是這樣的：

　　上式中的冪級數稱爲f(x)在x₀點的泰勒級數。（0的階乘是1）

  　 更多泰勒公式的介紹可參考  單變量微積分筆記31——冪級數和泰勒級數

 

泰勒公式的應用

　　來看一個泰勒公式的應用。假設一個小偷盜取了一輛汽車，他在高速公路上沿着一個方向行駛，車輛的位移s是關於時間t的函數。警方接到報案後立刻調取監控，得知在零點（t=0時刻）小偷距車輛丟失地點的位移是s₀。如今的時間是0:30，警方想要在前方設卡，從而能在凌晨1點攔住小偷，應該在哪裏設卡呢？

　　咱們知道車輛在0點時的位移是s₀，如今想要知道凌晨1點時車輛的位置：

　　能夠直接使用泰勒公式：

　　泰勒公式能夠無限展開，展開得越多，越逼近真實值，而且越到後面的項，對結果的影響越小，咱們認爲0和1很是接近，因此只展開到2階導數：

　　這就是最終結果，在此處設卡最有可能在第一時間攔住小偷。

在0點處的泰勒展開

　　在使用泰勒公式時，常常取x₀=0。

　　f(x)=e^x是一個能夠用泰勒公式展開的例子，下面是e^x在x₀=0處的泰勒展開：

　　當x=1時，還附帶獲得了e的解釋：

　　咱們使用一個很難處理的積分解釋泰勒展開的意義，對正態分佈進行積分：

　　常規的方法很難處理。如今，因爲被積函數與e^x類似，咱們又已經知道e^x的展開式，因此能夠進行下面的變換：

　　將exp(-x²)左右兩側同時積分：

　　很容易計算右側的每一項積分。

　　這個例子展現了冪級數展開的意義——把質的困難轉化成量的複雜。展開前求解函數的值很困難，展開後是冪級數，雖然有不少不少項，可是每一項都是冪函數，都很容易求解，因而，只要對展開後的函數求和，就能獲得展開前的函數的值。

爲何在0點處展開

　　當x₀=0時，能夠極大地簡化泰勒展開式。以前說泰勒公式是一個用函數某點的信息描述在該點附近取值的公式，一個函數中的某點若是距離0很遠怎麼辦呢？實際上泰勒公式也可以逼近函數在距離0很遠處的取值，只不過此時只展開到2階導數是不夠的，須要展開不少項，展開的越多，越能逼近該點。以標準正態分佈函數f(x)=exp(-x²)爲例，雖然它在二階展開使與原函數相差較大，可是當展開到40階時就已經很是接近原函數了。

多項式函數

　　理解了泰勒公式後，再回到問題的原點，看看多項式函數爲何能夠接近於任何函數。

　　仍然以標準正態分佈爲例，它在x₀ = 0點處的10階泰勒展開是：

　　若是將每一項中的xⁱ都看做一個維度，那麼這個多項式函數能夠寫成多元線性迴歸的形式：

　　這就將一個一元的非線性問題轉換成了多元的線性問題，從而利用最小二乘法求得模型參數。

　　下面的代碼以ln(2x) + 2爲原函數，生成40個在-1~1之間隨機震盪的數據點，並使用線性迴歸和多項式迴歸擬合數據點：

 

  1 import numpy as np
  2 import matplotlib.pyplot as plt
  3 
  4 def create_datas():
  5     '''
  6     生成10個待擬合的點
  7     :return: xs, ys
  8     '''
  9     xs = np.arange(0.1, 4, 0.1)
 10     # y = ln(2x) + noize,  -1 <= noize <= 1
 11     ys = np.array([np.log(x * 2) + 2 + np.random.uniform(-1, 1) for x in xs])
 12 
 13     return xs, ys
 14 
 15 class Regression():
 16     ''' 迴歸類 '''
 17     def __init__(self, xs, ys):
 18         '''
 19         :param xs: 輸入數據的特徵集合
 20         :param ys: 輸入數據的標籤集合
 21         '''
 22         self.xs, self.ys = xs, ys
 23         self.theta = None # 模型參數
 24 
 25     def train_datas(self, xs=None):
 26         '''
 27         從新構造訓練樣本的特徵和標籤
 28         :param xs: 輸入數據的特徵集合
 29         :return: 矩陣形式的訓練樣本特徵和標籤
 30         '''
 31         xs = self.xs if xs is None else xs
 32         X = self.train_datas_x(xs)
 33         Y = np.c_[ys] # 將ys轉換爲m行1列的矩陣
 34         return X, Y
 35 
 36     def train_datas_x(self, xs):
 37         '''
 38         從新構造訓練樣本的特徵
 39         :param xs: 輸入數據的特徵集合
 40         :return: 矩陣形式的訓練樣本特徵
 41         '''
 42         m = len(xs)
 43         # 在第一列添加x0,x0=1,並將二維列表轉換爲矩陣
 44         X = np.mat(np.c_[np.ones(m), xs])
 45         return X
 46 
 47     def fit(self):
 48         ''' 數據擬合 '''
 49         X, Y = self.train_datas()
 50         self.theta = (X.T * X).I * X.T * Y
 51 
 52     def predict(self, xs):
 53         '''
 54         根據模型預測結果
 55         :param xs: 輸入數據的特徵集合
 56         :return: 預測結果
 57         '''
 58         X = self.train_datas(xs=xs)[0]
 59         return self.theta.T * X.T
 60 
 61     def show(self):
 62         ''' 繪製擬合結果 '''
 63         plt.figure()
 64         plt.scatter(self.xs, self.ys, color='r', marker='.', s=10)  # 繪製數據點
 65         self.show_curve(plt) # 繪製函數曲線
 66         plt.xlabel('x')
 67         plt.ylabel('y')
 68         plt.axis('equal')
 69         plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤
 70         plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的座標軸負號顯示問題
 71         plt.legend(['擬合曲線', '樣本點'])
 72         plt.show()
 73 
 74     def show_curve(self, plt):
 75         ''' 繪製函數曲線 '''
 76         pass
 77 
 78     def global_fun(self):
 79         ''' 返回目標函數 '''
 80         gf = ['(' + str(t[0, 0]) + str(i) + ')x^' + str(i) for i, t in enumerate(self.theta)]
 81         return ' + '.join(gf)
 82 
 83 class Linear(Regression):
 84     ''' 線性模型'''
 85     def show_curve(self, plt):
 86         '''
 87         繪製擬合結果
 88         :param plt: 輸入數據的特徵集合
 89         '''
 90         xx = [self.xs[0], self.xs[-1]]
 91         yy = self.predict(xx)
 92         plt.plot(xx, np.array(yy)[0])
 93 
 94 class Multinomial(Regression):
 95     ''' 多項式迴歸模型 '''
 96     def __init__(self, xs, ys, n=3):
 97         '''
 98         :param xs: 輸入數據的特徵集合
 99         :param ys: 輸入數據的標籤集合
100         :param n: 多項式的項數
101         '''
102         super().__init__(xs, ys)
103         self.n = n
104 
105     def train_datas_x(self, xs):
106         '''
107         從新構造訓練樣本的特徵
108         :param xs: 輸入數據的特徵集合
109         :return: 矩陣形式的訓練樣本特徵
110         '''
111         X = super().train_datas_x(xs)
112         for i in range(2, self.n + 1):
113             X = np.column_stack((X, np.c_[xs ** i])) # 構造樣本的其餘特徵
114         return X
115 
116     def show_curve(self, plt):
117         ''' 繪製函數曲線 '''
118         xx = np.linspace(self.xs[0], self.xs[-1], len(self.xs) * 20)
119         yy = self.predict(xx)
120         plt.plot(xx, np.array(yy)[0], '-')
121 
122 if __name__ == '__main__':
123     xs, ys = create_datas()
124     regressions = [Linear(xs, ys), Multinomial(xs, ys), Multinomial(xs, ys, n=5), Multinomial(xs, ys, n=10)]
125     for r in regressions:
126         r.fit()
127         r.show()
128         print(r.global_fun())

(1.702537204930)x^0 + (0.75431357262260011)x^1

 

 (0.23422131704216660)x^0 + (3.8713793437217621)x^1 + (-1.51749485964066682)x^2 + (0.206815637166500283)x^3

 (0.0023811193415048670)x^0 + (4.707160334405161)x^1 + (-2.03334533257402762)x^2 + (0.095635349482284143)x^3 + (0.130611330518000564)x^4 + (-0.021122013844903465)x^5

(-4.7285135624557920)x^0 + (77.637488456533421)x^1 + (-377.238590224254552)x^2 + (932.32693158635363)x^3

+ (-1305.30725871564164)x^4 + (1112.9257341435945)x^5 + (-598.57958115210336)x^6

+ (203.91275172701427)x^7 + (-42.641981259587898)x^8 + (4.9915417588645349)x^9 + (-0.250300601937088710)x^10

 　　看來第2、三條曲線的擬合效果比較好，第一幅圖欠擬合，四過擬合。

---

　　做者：我是8位的

　　出處：http://www.cnblogs.com/bigmonkey

　　本文以學習、研究和分享爲主，如需轉載，請聯繫本人，標明做者和出處，非商業用途！ 

　　掃描二維碼關注公做者衆號「我是8位的」