萬字總結複雜而奇妙的高斯過程!

2021-02-20 數據分析1480


高斯過程的理論知識

高斯過程的Python實現

使用Numpy手動實現

使用`Scikit-learn`實現高斯過程

小結


高斯過程GaussianProcess
高斯過程的理論知識非參數方法的基本思想

在監督學習中，我們經常使用包含參數

來解釋數據，並通過 最大似然估計（MLE） 或 最大後驗估計（MPE） 來推斷參數

如果需要，我們還可以推斷出一個完整的後驗分布

來代替點估計

隨著數據複雜性的增加，通常需要使用具有更多參數的模型來合理地解釋數據。在非參數方法中，參數的數量取決於數據集的大小。

例如，在Nadaraya-Watson 核回歸（Kernel regression） 中，將權重

這裡，我們假設了接近 核函數

一種特殊情況是k最近鄰（KNN），其中最接近

非參數方法通常需要處理所有訓練數據以進行預測，因此相比於參數方法，推斷(inference)時速度較慢。另一方面，由於非參數模型僅需要記住訓練數據，在這個層面上來說，訓練速度通常會更快。

注1： 核回歸是一種非參數回歸的方法，這種模型是基於（特徵注2：

高斯過程的基本概念

非參數方法的另一個示例是高斯過程（GPs）。參數方法需要推斷參數的分布，而在非參數方法中，比如高斯過程，它可以直接推斷函數的分布。

高斯過程定義了先驗函數。觀察到某些函數值後，可通過代數運算以將其轉換為後驗函數。

在這種情況下，連續函數值的推斷稱為GP回歸，但GP也可以用於分類。

高斯過程是一種隨機過程，其為任意點

                                                 

在這裡

並且 均值函數，通常使用正定核函數或協方差函數。因此，高斯過程是函數的分布，其形狀（平滑度等等性質）由核函數

這裡有一個假設：

如果通過核函數將點

觀察到一些GP的先驗

         

其中

根據GP的定義，觀測數據

 

給定

其中

其中

其中

 

這是實現高斯過程並將其應用於回歸問題所需的最低知識。下一部分將展示如何從頭開始使用純NumPy實現GP，隨後的部分將展示如何使用scikit-learn中的GP實現。高斯過程的Python實現使用Numpy手動實現定義核函數

在這裡，我們將使用平方指數內核squared exponential kernel，也稱為高斯內核或RBF內核：

 

長度參數 
import numpy as np

def kernel(X1, X2, l=1.0, sigma_f=1.0):
    '''
    各向同性平方指數內核.
    計算點X1與點X2的協方差矩陣.
    
    Args:
        X1: ndArray， m個點 (m x d).
        X2: ndArray， n個點 (n x d).

    返回:
        協方差矩陣 (m x n).
    '''
    sqdist = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
    return sigma_f**2 * np.exp(-0.5 / l**2 * sqdist)
定義先驗
我們首先定義一個均值 0 和一個用內核參數 
為了從該GP提取隨機函數，我們從相應的多元高斯分布抽取隨機樣本。
以下示例繪製了三個隨機樣本，並將其與 0 均值和95％置信區間（根據協方差矩陣的對角線計算）一起繪製。
# 編寫繪圖函數
import matplotlib.pyplot as plt

from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

def plot_gp(mu, cov, X, X_train=None, Y_train=None, samples=[]):
    X = X.ravel()
    mu = mu.ravel()
    uncertainty = 1.96 * np.sqrt(np.diag(cov))
    
    plt.fill_between(X, mu + uncertainty, mu - uncertainty, alpha=0.1)
    
    plt.plot(X, mu, label='均值')
    for i, sample in enumerate(samples):
        plt.plot(X, sample, lw=1, ls='--', label=f'樣本 {i+1}')
    if X_train is not None:
        plt.plot(X_train, Y_train, 'rx')
    plt.legend(bbox_to_anchor=(1.04,0.5), loc="center left")

def plot_gp_2D(gx, gy, mu, X_train, Y_train, title, i):
    ax = plt.gcf().add_subplot(1, 2, i, projection='3d')
    ax.plot_surface(gx, gy, mu.reshape(gx.shape), cmap=cm.coolwarm, linewidth=0, alpha=0.2, antialiased=False)
    ax.scatter(X_train[:,0], X_train[:,1], Y_train, c=Y_train, cmap=cm.coolwarm)
    ax.set_title(title)
%matplotlib inline

# 有限個輸入數據點
X = np.arange(-5, 5, 0.2).reshape(-1, 1)

# 先驗的均值與方差
mu = np.zeros(X.shape)
cov = kernel(X, X)

# 從先驗分布（多元高斯分布）中抽取樣本點
samples = np.random.multivariate_normal(mu.ravel(), cov, 3)

# 畫出GP的均值, 置信區間
plot_gp(mu, cov, X, samples=samples)
基於無噪聲訓練數據進行預測
為了計算充分統計量，即後驗預測分布的均值和協方差矩陣，我們用下面代碼實現公式（4）和（5）
# 倒入計算逆矩陣的函數inv()
from numpy.linalg import inv

def posterior_predictive(X_s, X_train, Y_train, l=1.0, sigma_f=1.0, sigma_y=1e-8):
    '''
    計算後驗分布的充分統計量 
    給定 m 個訓練數據 X_train 與 Y_train 
    給定 n 個新輸入數據 inputs X_s.
    
    Args:
        X_s: 新輸入數據 (n x d).
        X_train: 訓練輸入數據 (m x d).
        Y_train: 訓練輸出數據 (m x 1).
        l: 核函數的長度參數.
        sigma_f: 核函數的縱向波動參數.
        sigma_y: 噪音參數.
    
    返回:
        後驗均值向量 (n x d) 與協方差矩陣 (n x n).
    '''
    K = kernel(X_train, X_train, l, sigma_f) + sigma_y**2 * np.eye(len(X_train))
    K_s = kernel(X_train, X_s, l, sigma_f)
    K_ss = kernel(X_s, X_s, l, sigma_f) + 1e-8 * np.eye(len(X_s))
    K_inv = inv(K)
    
    # 公式 (4)
    mu_s = K_s.T.dot(K_inv).dot(Y_train)

    # 公式 (5)
    cov_s = K_ss - K_s.T.dot(K_inv).dot(K_s)
    
    return mu_s, cov_s
並將它們應用於無噪聲訓練數據X_train和Y_train。以下示例從後驗預測中提取3個樣本，並將它們與均值，置信區間和訓練數據一起繪製。在無噪聲模型中，訓練點的方差為 0
註： 從後驗預測分布提取的所有隨機函數都經過訓練點。
# 無噪音的5個輸入數據
X_train = np.array([-4, -3, -2, -1, 1]).reshape(-1, 1)
# y=sin(x)
Y_train = np.sin(X_train)

# 計算後驗預測分布的均值向量與協方差矩陣
mu_s, cov_s = posterior_predictive(X, X_train, Y_train)

# 從後驗預測分布中抽取3個樣本
samples = np.random.multivariate_normal(mu_s.ravel(), cov_s, 3)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train, samples=samples)
基於帶有噪音的訓練數據進行預測
如果模型中包含一些噪聲，則僅對訓練點進行近似，並且訓練點的方差不為零。
# 定義噪音參數
noise = 0.4

# 帶有噪音的訓練數據
X_train = np.arange(-3, 4, 1).reshape(-1, 1)
Y_train = np.sin(X_train) + noise * np.random.randn(*X_train.shape)

# 可以對比地看一下帶噪音的訓練數據與不帶噪音的訓練數據的區別
plt.figure()
plt.plot(X_train, np.sin(X_train) + noise * np.random.randn(*X_train.shape), lw=1, ls='-.', label='有噪音')
plt.plot(X_train, np.sin(X_train) + 0.0 * np.random.randn(*X_train.shape), lw=2.5, ls='-', label='無噪音')
plt.plot(X_train, np.sin(X_train) + 0.0 * np.random.randn(*X_train.shape), 'rx')
plt.legend(bbox_to_anchor=(1.04,0.5), loc="center left")
plt.show()
# 計算後驗預測分布的均值向量以及方差矩陣
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, sigma_y=noise)

# 從後驗預測分布中抽取3個樣本點
samples = np.random.multivariate_normal(mu_s.ravel(), cov_s, 3)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train, samples=samples)
核函數參數和噪聲參數的影響
以下示例顯示了核函數參數 
（下圖第一行）
（下圖第二行）
（下圖第三行）
params = [
    (0.3, 1.0, 0.2),
    (3.0, 1.0, 0.2),
    (1.0, 0.3, 0.2),
    (1.0, 3.0, 0.2),
    (1.0, 1.0, 0.05),
    (1.0, 1.0, 1.5),
]

plt.figure(figsize=(12, 5))

for i, (l, sigma_f, sigma_y) in enumerate(params):
    mu_s, cov_s = posterior_predictive(X, X_train, Y_train, l=l, 
                                       sigma_f=sigma_f, 
                                       sigma_y=sigma_y)
    plt.subplot(3, 2, i + 1)
    plt.subplots_adjust(top=2)
    plt.title(f'l = {l}, sigma_f = {sigma_f}, sigma_y = {sigma_y}')
    plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)
這些參數的最優值可以通過最大化由[1] [3]給出的邊際對數似然來得到：
 
在下面的代碼中，我們將最小化負邊際對數似然來獲得核函數參數 
from numpy.linalg import cholesky, det, lstsq
from scipy.optimize import minimize

def nll_fn(X_train, Y_train, noise, naive=True):
    '''
    基於給定數據X_train和Y_train以及噪聲水平
    返回一個可以計算負邊際對數似然的函數
    
    Args:
        X_train: 訓練輸入 (m x d).
        Y_train: 訓練輸出 (m x 1).
        noise: Y_train的噪聲水平.
        naive: 如果 True 那麼使用公式(7)來實現
               如果 False 那麼使用數值方法來實現
               
        
    返回:
        最小的目標對象
    '''
    def nll_naive(theta):
        # 使用公式(7)來實現 
        # 與下面的nll_stable的實現相比在數值上不穩定
        K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1]) + \
            noise**2 * np.eye(len(X_train))
        return 0.5 * np.log(det(K)) + \
               0.5 * Y_train.T.dot(inv(K).dot(Y_train)) + \
               0.5 * len(X_train) * np.log(2*np.pi)

    def nll_stable(theta):
        # 數值上更穩定，相比於nll_naive 
        K = kernel(X_train, X_train, l=theta[0], sigma_f=theta[1]) + \
            noise**2 * np.eye(len(X_train))
        L = cholesky(K)
        return np.sum(np.log(np.diagonal(L))) + \
               0.5 * Y_train.T.dot(lstsq(L.T, lstsq(L, Y_train)[0])[0]) + \
               0.5 * len(X_train) * np.log(2*np.pi)
    
    if naive:
        return nll_naive
    else:
        return nll_stable

# 求解可滿足最小化目標函數的參數 l 及 sigma_f.
# 實際上，我們應該使用不同的初始化多次運行最小化，以避免局部最小化，
# 但是為了簡單起見，此處將其跳過
res = minimize(nll_fn(X_train, Y_train, noise), [1, 1], 
               bounds=((1e-5, None), (1e-5, None)),
               method='L-BFGS-B')

# 將優化結果存儲在全局變量中，以便我們以後可以將其與其他實現的結果進行比較
l_opt, sigma_f_opt = res.x

# 使用優化的核函數參數計算後驗預測分布的參數，並繪製結果圖
mu_s, cov_s = posterior_predictive(X, X_train, Y_train, l=l_opt, sigma_f=sigma_f_opt, sigma_y=noise)
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)
其最優化的核函數參數 
print(l_opt, sigma_f_opt)
0.9872536793237083 0.8613778055591963
更高維的高斯過程
以上實現也可以用於更高的輸入數據維度。此處，GP用於擬合在 
下圖顯示了核函數參數優化前後的帶噪聲樣本和後驗預測分布的均值向量。
#噪音參數
noise_2D = 0.1

rx, ry = np.arange(-5, 5, 0.3), np.arange(-5, 5, 0.3)
gx, gy = np.meshgrid(rx, rx)

X_2D = np.c_[gx.ravel(), gy.ravel()]

X_2D_train = np.random.uniform(-4, 4, (100, 2))
Y_2D_train = np.sin(0.5 * np.linalg.norm(X_2D_train, axis=1)) + \
             noise_2D * np.random.randn(len(X_2D_train))

plt.figure(figsize=(14,7))

mu_s, _ = posterior_predictive(X_2D, X_2D_train, Y_2D_train, sigma_y=noise_2D)
plot_gp_2D(gx, gy, mu_s, X_2D_train, Y_2D_train, 
           f'參數優化前: l={1.00} sigma_f={1.00}', 1)

res = minimize(nll_fn(X_2D_train, Y_2D_train, noise_2D), [1, 1], 
               bounds=((1e-5, None), (1e-5, None)),
               method='L-BFGS-B')

mu_s, _ = posterior_predictive(X_2D, X_2D_train, Y_2D_train, *res.x, sigma_y=noise_2D)
plot_gp_2D(gx, gy, mu_s, X_2D_train, Y_2D_train,
           f'參數優化後: l={res.x[0]:.2f} sigma_f={res.x[1]:.2f}', 2)
使用Scikit-learn實現高斯過程
scikit-learn 提供了 GaussianProcessRegressor方法來實現GP回歸模型. 你可以自定義核函數，或者使用它內置的核函數. 核函數也可以疊加. Squared exponential 核函數在scikit-learn中也就是RBF. scikit-learn中的RBF只有長度參數 
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF

rbf = ConstantKernel(1.0) * RBF(length_scale=1.0)
gpr = GaussianProcessRegressor(kernel=rbf, alpha=noise**2)

# 1D 訓練樣本
gpr.fit(X_train, Y_train)

# 計算後驗預測分布的均值向量與協方差矩陣
mu_s, cov_s = gpr.predict(X, return_cov=True)

# 獲得最優核函數參數
l = gpr.kernel_.k2.get_params()['length_scale']
sigma_f = np.sqrt(gpr.kernel_.k1.get_params()['constant_value'])

# 與前面手寫的結果比對
assert(np.isclose(l_opt, l))
assert(np.isclose(sigma_f_opt, sigma_f))

# 繪製結果
plot_gp(mu_s, cov_s, X, X_train=X_train, Y_train=Y_train)
小結
從前面我們可以看出，與常見的機器學習模型不同，用高斯過程做預測的方法是直接生成一個後驗預測分布（依然是高斯分布）。
這也決定了我們可以不僅僅得到一個「光禿禿」的預測值，還可以得到關於這個預測值的不確定性信息，可以利用這些信息繪製error-bar等等！
從統計學的角度上來看，利用高斯過程模型做預測具有很高的價值。

相關焦點

  • 高斯過程
    高斯過程的理論知識高斯過程的Python實現使用Numpy手動實現使用`Scikit-learn`實現高斯過程小結高斯過程GaussianProcess高斯過程的理論知識非參數方法的基本思想非參數方法的另一個示例是高斯過程（GPs）。
  • 高斯過程回歸(GPR)
    我們先從高斯過程開始介紹，理解這兩個視角下的高斯過程回歸，你會發現它其實是貝葉斯線性回歸的非線性形式，且為後面深入理解貝葉斯調參打下基礎。高斯過程（Gaussian Process,GP）是定義在連續域上的無限多個服從高斯分布的隨機變量所組成的隨機過程。
  • 重慶高斯小數
    教師根據課件引導學生將每一個小問題答對，最後再次回到今天要講的問題，做一個強化總結。④檢驗環節講解完知識點之後，增加快問快答環節，通過單步訓練的方式，讓孩子在微緊張的狀態下，對剛吸收的內容在腦海中進行快速的處理。這樣的狀態也容易讓孩子對課堂保持高度集中。
  • NeurIPS 2020 | 近期必讀高斯過程精選論文
    在高斯過程中，連續輸入空間中每個點都是與一個正態分布的隨機變量相關聯。此外，這些隨機變量的每個有限集合都有一個多元正態分布，換句話說他們的任意有限線性組合是一個正態分布。高斯過程的分布是所有那些隨機變量的聯合分布，正因如此，它是連續域上函數的分布。
  • 深度科普:說說高斯過程回歸
    高斯過程是機器學習領域一個基礎的方法，同時又和其他方法有千絲萬縷的聯繫，值得大家研究一下。文中一些細節也歡迎大家和作者一起探討。網上講高斯過程回歸的文章很少，且往往從高斯過程講起，我比較不以為然：高斯過程回歸（GPR）， 終究是個離散的事情，用連續的高斯過程( GP) 來闡述，簡直是殺雞用牛刀。所以我們這次直接從離散的問題搞起，然後把高斯過程逆推出來。這篇博客的主要目的是解釋高斯過程回歸這個主意是怎麼想出來的，模型多了去了，為毛要用它。
  • 一文了解高斯濾波器,附原理及實現過程
    生成的過程，首先根據模板的大小，找到模板的中心位置ksize/2。然後就是遍歷，根據高斯分布的函數，計算模板中每個係數的值。需要注意的是，最後歸一化的過程，使用模板左上角的係數的倒數作為歸一化的係數(左上角的係數值被歸一化為1)，模板中的每個係數都乘以該值(左上角係數的倒數)，然後將得到的值取整，就得到了整數型的高斯濾波器模板。
  • 一種利用推斷網絡對高斯過程模型進行有效推斷的方法
    對於本文使用的場景即在函數空間解決高斯過程推斷問題，隨機變量是函數，因此是無限維的很難得到一個解，作者便使用了推斷網絡，並通過上式計算隨機梯度並更新推斷網絡，這在下一部分進行詳細介紹。4.2 使用minibatch對推斷網絡進行訓練。
  • 圖文詳解高斯過程(一)——含代碼
    為了幫助入門者更好地理解這一簡單易用的方法，近日國外機器學習開發者Alex Bridgland在博客中圖文並茂地解釋了高斯過程，並授權論智將文章分享給中國讀者。註：本文為系列第一篇，雖用可視化形式弱化了數學推導，但仍假設讀者具備一定機器學習基礎。現如今，高斯過程可能稱不上是機器學習領域的炒作核心，但它仍然活躍在研究的最前沿。
  • 高斯模糊算法的全面優化過程分享(二)
    在高斯模糊算法的全面優化過程分享(一)一文中我們已經給出了一種相當高性能的高斯模糊過程，但是優化沒有終點，經過上一個星期的發憤圖強和測試，對該算法的效率提升又有了一個新的高度，這裡把優化過程中的一些心得和收穫用文字的形式記錄下來。
  • 正態分布和高斯分布的作用_高斯分布的定義_誤差服從高斯分布
    打開APP 正態分布和高斯分布的作用_高斯分布的定義_誤差服從高斯分布 發表於 2017-12-04 16:38:44 　　高斯分布的定義 　　高斯分布怎麼來的，很簡單。只要你觀察的系統裡，各種對象之間關聯很弱，那麼他們的總和平均表現，根據中心極限定律，就是高斯或者近高斯的。
  • 《高斯數學》課程
    一、課程簡介高斯數學是江蘇高斯教育以美國著名心理學家吉爾福特的「智力結構理論
  • 數學家高斯
    能夠在頭腦中進行複雜的計算，是上帝賜予他一生的天賦。父親格爾恰爾德·迪德裡赫對高斯要求極為嚴厲，甚至有些過分。高斯尊重他的父親，並且秉承了其父誠實、謹慎的性格。高斯很幸運地有一位鼎力支持他成才的母親。高斯一生下來，就對一切現象和事物十分好奇，而且決心弄個水落石出，這已經超出了一個孩子能被許可的範圍。當丈夫為此訓斥孩子時，她總是支持高斯，堅決反對頑固的丈夫想把兒子變得跟他一樣無知。
  • 從數學到實現,全面回顧高斯過程中的函數最優化
    高斯過程可以被認為是一種機器學習算法，它利用點與點之間同質性的度量作為核函數，以從輸入的訓練數據預測未知點的值。本文從理論推導和實現詳細地介紹了高斯過程，並在後面提供了用它來近似求未知函數最優解的方法。
  • 高斯和小行星的發現
    高斯出生在一個貧苦的家庭裡，祖父是農民，父親做短工，舅舅腓特烈(Friederich) 是一個很有才能的入，自已學會了紡織技術，很快成為一個出色的錦緞織工。他經常教給高斯一些知識，對幼年的高斯影響很大.高斯父親本不打算讓他上學，但高斯很小就顯出有數學才能。
  • 德國數學家高斯的有趣故事
    享有「數學王子」的美譽的高斯，他是近代數學的奠基人。高斯畫像高斯（1777~1855）出生在德國的一個貧困家庭，雖然出生的環境不好，但小高斯在3歲時就喜歡數數，心算能力超出同齡小夥伴們許多，據說在他三歲的時，有一天晚上他在看著父親在反覆數工錢
  • 「數學王子」高斯
    是近代數學奠基者之一，高斯被認為是歷史上最重要的數學家之一，並享有「數學王子」之稱。高斯和阿基米德、牛頓並列為世界三大數學家。一生成就極為豐碩，以他名字「高斯」命名的成果達110個，屬數學家中之最。他對數論、代數、統計、分析、微分幾何、大地測量學、地球物理學、力學、靜電學、天文學、矩陣理論和光學皆有貢獻。 高斯是一對貧窮夫婦的唯一的兒子。
  • 【小新星·推薦】高斯小數,讓孩子愛上數學,主動思考
    教師根據課件引導學生將每一個小問題答對，最後再次回到今天要講的問題，做一個強化總結。講解完知識點之後，增加快問快答環節，通過單步訓練的方式，讓孩子在微緊張的狀態下，對剛吸收的內容在腦海中進行快速的處理。這樣的狀態也容易讓孩子對課堂保持高度集中。升級前的板書是常規的寫字標註形式，升級後通過動態顯示模式，將總結的答案自動顯示到板書上，讓授課更加自動化。
  • 《數學博覽》 數學家高斯
    高理的數論研究 總結 在《算術研究》（1801）中，這本書奠定了近代數論的基礎，它不僅是數論方面的劃時代之作，也是數學史上不可多得的經典著作之一。高斯對代數學的重要貢獻是證明了代數基本定理，他的存在性證明開創了數學研究的新途徑。高斯在1816年左右就得現了著名的柯西積分定理。他還發現橢圓函數的雙周期性，但這些工作在他生前都沒發表出來。
  • 為什麼數據科學家都喜歡高斯分布
    對深度學習和機器學習工程師而言，在世界上所有的概率模型中，高斯分布（Gaussian distribution）模型最為引人注目。即使你從來沒有進行過AI項目，有很大的機率你曾經遇到過高斯模型。高斯分布，又稱為正態分布（Normal distribution），常常可以通過其標誌性的鐘形曲線識別出來。高斯分布如此流行，有三大原因。
  • 高斯數學簡介
    高斯數學是江蘇高斯教育以《3-6歲兒童學習與發展指南》為依據，根據幼兒的心裡、生理特點及認知規律，研發了適合4—7幼兒的系統、科學的學前課程。三、課程理念   高斯數學是幼教專家引領、一線教師二十多年教學經驗的總結和提煉；強大的專家團隊、教師同行的實操驗證；幼兒園的落地使用證明了它的系統性、科學性、可行性。    高斯數學注重「雙重」培養。既注重快速計算能力的訓練，又注重數學思維能力的培養。