最小二乘法（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之間隨機震盪的數據點，並使用線性回歸和多項式回歸擬合數據點：

 

import numpy as np
import matplotlib.pyplot as plt

def create_datas():
    '''
    生成10個待擬合的點
    :return: xs, ys
    '''
    xs = np.arange(0.1, 4, 0.1)
    # y = ln(2x) + noize,  -1 <= noize <= 1
    ys = np.array([np.log(x * 2) + 2 + np.random.uniform(-1, 1) for x in xs])

    return xs, ys

class Regression():
    ''' 回歸類 '''
    def __init__(self, xs, ys):
        '''
        :param xs: 輸入數據的特征集合
        :param ys: 輸入數據的標簽集合
        '''
        self.xs, self.ys = xs, ys
        self.theta = None # 模型參數

    def train_datas(self, xs=None):
        '''
        重新構造訓練樣本的特征和標簽
        :param xs: 輸入數據的特征集合
        :return: 矩陣形式的訓練樣本特征和標簽
        '''
        xs = self.xs if xs is None else xs
        X = self.train_datas_x(xs)
        Y = np.c_[ys] # 將ys轉換為m行1列的矩陣
        return X, Y

    def train_datas_x(self, xs):
        '''
        重新構造訓練樣本的特征
        :param xs: 輸入數據的特征集合
        :return: 矩陣形式的訓練樣本特征
        '''
        m = len(xs)
        # 在第一列添加x0,x0=1,並將二維列表轉換為矩陣
        X = np.mat(np.c_[np.ones(m), xs])
        return X

    def fit(self):
        ''' 數據擬合 '''
        X, Y = self.train_datas()
        self.theta = (X.T * X).I * X.T * Y

    def predict(self, xs):
        '''
        根據模型預測結果
        :param xs: 輸入數據的特征集合
        :return: 預測結果
        '''
        X = self.train_datas(xs=xs)[0]
        return self.theta.T * X.T

    def show(self):
        ''' 繪制擬合結果 '''
        plt.figure()
        plt.scatter(self.xs, self.ys, color='r', marker='.', s=10)  # 繪制數據點
        self.show_curve(plt) # 繪制函數曲線
        plt.xlabel('x')
        plt.ylabel('y')
        plt.axis('equal')
        plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標簽
        plt.rcParams['axes.unicode_minus'] = False  # 解決中文下的坐標軸負號顯示問題
        plt.legend(['擬合曲線', '樣本點'])
        plt.show()

    def show_curve(self, plt):
        ''' 繪制函數曲線 '''
        pass

    def global_fun(self):
        ''' 返回目標函數 '''
        gf = ['(' + str(t[0, 0]) + str(i) + ')x^' + str(i) for i, t in enumerate(self.theta)]
        return ' + '.join(gf)

class Linear(Regression):
    ''' 線性模型'''
    def show_curve(self, plt):
        '''
        繪制擬合結果
        :param plt: 輸入數據的特征集合
        '''
        xx = [self.xs[0], self.xs[-1]]
        yy = self.predict(xx)
        plt.plot(xx, np.array(yy)[0])

class Multinomial(Regression):
    ''' 多項式回歸模型 '''
    def __init__(self, xs, ys, n=3):
        '''
        :param xs: 輸入數據的特征集合
        :param ys: 輸入數據的標簽集合
        :param n: 多項式的項數
        '''
        super().__init__(xs, ys)
        self.n = n

    def train_datas_x(self, xs):
        '''
        重新構造訓練樣本的特征
        :param xs: 輸入數據的特征集合
        :return: 矩陣形式的訓練樣本特征
        '''
        X = super().train_datas_x(xs)
        for i in range(2, self.n + 1):
            X = np.column_stack((X, np.c_[xs ** i])) # 構造樣本的其他特征
        return X

    def show_curve(self, plt):
        ''' 繪制函數曲線 '''
        xx = np.linspace(self.xs[0], self.xs[-1], len(self.xs) * 20)
        yy = self.predict(xx)
        plt.plot(xx, np.array(yy)[0], '-')

if __name__ == '__main__':
    xs, ys = create_datas()
    regressions = [Linear(xs, ys), Multinomial(xs, ys), Multinomial(xs, ys, n=5), Multinomial(xs, ys, n=10)]
    for r in regressions:
        r.fit()
        r.show()
        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

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

---

　　作者：我是8位的

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

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

　　掃描二維碼關注公作者眾號“我是8位的”