Kriging模型理論推導

• 1、前言
• 2、條件
• 3、基礎知識
• • • 3.1、方差的理解
    • 3.2、概率密度函數
    • 3.3、多元正態分布
• 4、理論推導
• • • 4.1 模型建立
    • 4.2 模型預測

1、前言

簡介：Kriging模型是一種通過已知試驗點信息來預測未知試驗點上響應的無偏估計模型，其最早是由南非礦業工程師Ｄ.Ｇ.Krige于1951年提出。20世紀70年代，法國的數學家Ｇ.Matheron對Ｄ.Ｇ.Krige的研宄成果進行了進一步的系統化、理論化，并將其命名為Kriging模型。1989年Sacks等將Kriging模型推廣至試驗設計領域，形成了基于計算機仿真和Kriging模型的計算機試驗設計與分析方法。
本文將從原理部分，解析Kriging模型的推導過程。本次克里金模型的推導的參考文獻為：
?A Taxonomy of Global Optimization Methods Based on Response Surfaces。

2、條件

克里金模型在應用時有如下假設條件：
(1)、克里金法假設所有數據之間都服從n維的正態分布。
(2)、無偏。

3、基礎知識

在推導克里金模型之前，先來回顧一些統計學的基礎知識，各位功底深厚的看客老爺可以直接跳過。

3.1、方差的理解

概率論中方差用來度量隨機變量和其數學期望（即均值）之間的偏離程度。機器學習中方差又可以理解為不確定性的一種，即方差越大，不確定性越大。

3.2、概率密度函數

在數學中，連續型隨機變量的概率密度函數是描述這個隨機變量的輸出值，在某個確定的取值點附近的可能性的函數。而隨機變量的取值落在某個區域之內的概率則為概率密度函數在這個區域上的積分。當概率密度函數存在的時候，累積分布函數是概率密度函數的積分。

3.3、多元正態分布

平常我們見的最多的正態分布大多是是一維的，其的概率密度函數(Probability density function,PDF)如下：

其中，μ為均值，σ²為方差。也就是說，在均值和方差確定的條件下，上式f(x)也就確定了，這樣我們就可以知道在該分布下，隨機變量x的可能性大小。
同樣，當拓展到二維正態分布時，相當于添加了一個維度，這是均值仍然為每個維度上的均值組合到一起，而方差則變為了協方差，因為要考慮這兩個維度之間的關系。此時的均值和協方差變為：

此時二元正態分布的概率密度函數(pdf)為：

其中ρ為相關系數，是由這兩個維度上的方差計算得到的，如下圖所示：x是第一維度上的隨機變量，在該維度上，x服從正態分布，同樣的，y是第二維度上的隨機變量，在該維度上，同樣服從正態分布。z就是隨機變量x和y取某一個確定值的可能性大小。

以上為二元正態分布，多元正態分布也是類似，在增加維度即可，不過當維度超過2時，就無法可視化，但并不妨礙我們理解。
多元正態分布的均值向量為：

多元正態分布的協方差矩陣為：

其分布函數為：

也就是說，如果多元正態分布的均值確定了，協方差確定了，那么其分布函數(pdf)就可以確定，我們就可以在這個分布函數上搞點兒事情。比如進一步的進行最大似然估計。

4、理論推導

4.1 模型建立

已知給定了一些標記過的數據集X = { x₁,x₂,…,x_n },其對應的目標函數值為y = { y₁,y₂,…,y_n } ，注意，其中的 x₁ 是一個長度為 n 的向量，y₁ = Y(x₁) 。我們的目標就是想通過這些已知的點，來實現對未知點的預測。
首先，克里金模型假設所有數據服從均值為μ方差為σ²的n元的正態分布，也就是說這個n元正態分布函數的均值可以認為是在[ μ-3σ, μ+3σ ]的范圍內變化(論文原話，實際上刻畫的是不確定性)。現在我們考慮兩個點 x_i 和 x_j ，在我們采樣之前，是不確定這兩個點的目標函數值的，然而，我們假設建模所用的函數是連續的，當距離 || xi-xj || 比較小時，y(xi)和y(xj)也傾向于高度相關。我們可以通過下面的式子來衡量相關性：

上面是論文中的描述，初學者可能會比較蒙，下面我簡單解釋一下：
既然克里金模型假設了所有數據服從n維正態分布，那么對于n維的正態分布，如果想要刻畫其pdf，最重要的就是均值和協方差了了，由于是n維，均值為各個維度的均值的組合，為nx1的矩陣，而協方差矩陣里面，非對角線上的元素就是兩兩隨機變量之間的協方差，對角線上的元素就是各個隨機變量的方差，（如下圖示例，cov(z,x)刻畫的是變量z和變量x之間的相關性）。論文中的式(5),就是一個刻畫隨機變量Y(xi)和變量Y(xj)之間的相關性的函數，屬于協方差矩陣中的一員。我們令i=j，那么corr[Y(xi),Y(xj)]就為1.

緊接著，由于Y是服從n元正態分布，我們將n個已知點放到一起，就變成了：

Y的均值為lμ，其中l是nx1的矩陣。其協方差如下：

注意：文獻中得R乘以了方差，文獻作者應該是想表示協方差矩陣對角線上的值，不過不妨礙我們理解，這里我補充出Cov(Y)的表達式，見下圖：
$$?\begin{array}{ccc}Corr[Y(X1,X1)]& ...& Corr[Y(X1,Xn)]\\ ...& ...& ...\\ Corr[Y(Xn,X1)]& ...& Corr[Y(Xn,Xn)]\end{array}?$$
上式中的對角線上的值就是向量各自的方差。
這里的R的大小為n x n的矩陣，該矩陣中的每個值都是由公式(5)得到的，i和j都是從1取到n。對角線上i=j，所以R為1，那么協方差Cov(Y)的對角線就是方差。
由上公式可知超參數有μ、σ²，θ_l和p_l(l=1,2,3…d),我們用觀測數據y進行最大化似然來估計這些超參數，觀測數據y如下所示：

由于服從多維正態分布，最大似然的式子可以寫為：

為了方便運算，取對數：

下面分別對均值μ、和方差σ²求偏導，即可得到使似然函數最大的均值和方差了，得到結果如下：

最后，將公式(11)和(12)帶入到式(10)中得到log最大似然為：

由式13可知，log最大似然僅和R有關，而R中有參數θ，因此超參數的調節就是選取合適的θ使得log似然最大，可以用遺傳算法或多初始點算法求得。

4.2 模型預測

在4.1中，我們對已知得數據點進行了最大似然估計，得到一些先驗超參數，預測就是利用4.1得到得超參數來對未知數據點進行預測。這里考慮一個點$y$,我們將觀察到的數據和要預測的點放到一起$y$=(y^’,y^*)^T ,則對應的協方差矩陣也發生了改變：

則協方差矩陣變為：

矩陣中得 r^’ 實際是 r^T 的意思。則對應的似然函數為：
**式(10)**在加上下面圖片中的式子：

將$y$和$R$帶入到上式中得：

下面要做的是如何把中間的逆矩陣表示出來，這里作者用了部分求逆的方法，直接上結果如下：

$R$^-1 =

將上式帶入到式(16)中，我們可以得到擴充后的似然函數為：

我們可以看到，式(17)是關于y*的二次函數，對其求導并等于0可得：

從式(18)可以求解出：y*= $y$(x*) = $\mu$ + r^，R^-1(y-l$\mu$)

至此，證明完畢。

本文參考：

(1) https://zhuanlan.zhihu.com/p/90272131
(2) 文獻：A Taxonomy of Global Optimization Methods Based on Response Surfaces

總結

以上是生活随笔為你收集整理的克里金(kriging)模型的推导详解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。