多精度代理模型簡(jiǎn)明教程

引言

本文將簡(jiǎn)要介紹kriging代理模型和多層代理模型的原理及使用方法，以及對(duì)其程序進(jìn)行簡(jiǎn)要說(shuō)明。

目前隨著計(jì)算能力的不斷提升，我們已經(jīng)能夠?qū)σ恍┪锢磉^(guò)程進(jìn)行精細(xì)建模和高精度計(jì)算，但是隨著計(jì)算規(guī)模的提升，高精度計(jì)算帶來(lái)的時(shí)間和計(jì)算成本也常常難以讓人忍受。而使用代理模型可以有效減少計(jì)算量，同時(shí)獲取我們感興趣的參數(shù)。代理模型亦稱為響應(yīng)面模型或元模型，使用代理模型進(jìn)行設(shè)計(jì)優(yōu)化已經(jīng)成為重要的一種優(yōu)化方式。

與其他簡(jiǎn)單形式的代理模型，例如多項(xiàng)式模型和徑向基函數(shù)模型等相比，kriging模型具有如下兩個(gè)特點(diǎn)：1）預(yù)測(cè)精度高，2）可以給出預(yù)測(cè)準(zhǔn)確度估計(jì)。因此，在實(shí)際應(yīng)用中kriging模型正變得越來(lái)越受重視。許多文章和書籍都將kriging模型描述為一種“modelling the function as a realization of a stochastic process”，這是一種基于最小化期望平方誤差（minimize the expected squared prediction error）的模型，并且滿足預(yù)測(cè)結(jié)果是無(wú)偏估計(jì)，而且是其他已知結(jié)果$y_i$的線性組合。通過(guò)一系列初始數(shù)據(jù)可以得到kriging代理模型的參數(shù)，從而能夠構(gòu)建代理模型，這一過(guò)程稱為訓(xùn)練，而這些初始數(shù)據(jù)稱為訓(xùn)練集。

代理模型的預(yù)測(cè)精度取決于其訓(xùn)練集的精度，可以想象訓(xùn)練集的精度越高、數(shù)量越多，訓(xùn)練得到的代理模型精度也會(huì)越高。但是在實(shí)際應(yīng)用過(guò)程中，如果要獲取大量高精度訓(xùn)練集需要耗費(fèi)過(guò)多的時(shí)間和計(jì)算資源，因此為了解決這一問(wèn)題，我們可以考慮使用多精度代理模型（multi-fidelity surrogate model）。多精度代理模型的出發(fā)點(diǎn)是希望能夠用大量低精度的數(shù)據(jù)和少量高精度的數(shù)據(jù)構(gòu)建高精度代理模型，這樣就可以大大江都獲取高精度代理模型的難度。

接下來(lái)我們將簡(jiǎn)要介紹kriging代理模型、多精度kriging代理模型的原理和推導(dǎo)過(guò)程，詳細(xì)內(nèi)容請(qǐng)參見如下兩篇文獻(xiàn)：

• kriging代理模型：A Taxonomy of Global Optimization Methods Based on Response Surfaces by Donald R. Jones.
• 多精度kriging代理模型：Multi-fidelity optimization via surrogate modelling by Alexander I. J. Forrester, Andras Sobester and Andy J. Keane

1. kriging代理模型簡(jiǎn)介

1.1 概述

在工程界我們常常會(huì)遇到這樣一些問(wèn)題，希望通過(guò)建立輸入?yún)?shù)-輸出參數(shù)之間的映射關(guān)系，使得給定新的輸入來(lái)得到新的輸出，即預(yù)測(cè)值。根據(jù)輸出值的類型為連續(xù)型和離散型可以分為回歸問(wèn)題和分類問(wèn)題，處理回歸問(wèn)題的一種重要方法是高斯過(guò)程回歸(Gaussian Process Regression)，kriging代理模型是高斯過(guò)程回歸方法的一種應(yīng)用?；貧w問(wèn)題可以有如下數(shù)學(xué)描述：

假設(shè)有訓(xùn)練集$ D= {(\boldsymbol{x}_i,y_i)|i = 1,2,\cdots,n} = (\boldsymbol{X},\boldsymbol{y}). $其中：$\boldsymbol{x}_i \in R^d$為d維輸入矢量，$\boldsymbol{X} = [\boldsymbol{x}_1, \boldsymbol{x}2, \cdots, \boldsymbol{x}n]$ 為$d\times n$維矩陣，$y_i \in R$為相應(yīng)的輸出標(biāo)量，$\boldsymbol{y}$為輸出矢量?；貧w的任務(wù)是希望通過(guò)訓(xùn)練學(xué)習(xí)建立輸入$\boldsymbol{X}$與輸出$\boldsymbol{y}$之間的之間的映射關(guān)系$(f(\cdot) : R^d \mapsto R)$，預(yù)測(cè)出新測(cè)試點(diǎn)$\boldsymbol{x}*$對(duì)應(yīng)的最有可能的輸出值 $y*.$

高斯過(guò)程是任意有限個(gè)隨機(jī)變量具有聯(lián)合高斯分布的集合，其性質(zhì)完全由均值函數(shù)$m(\boldsymbol{x})$和協(xié)方差函數(shù)$K(\boldsymbol{x},\boldsymbol{x'})$決定，定義如下：
$$
m(\boldsymbol{x}) = E[f(\boldsymbol{x})] \
K(\boldsymbol{x},\boldsymbol{x'}) = E[(f(\boldsymbol{x}) - m(\boldsymbol{x})(f(\boldsymbol{x'}) - m(\boldsymbol{x'})]
$$
接下來(lái)我們簡(jiǎn)要推導(dǎo)一下kriging模型的實(shí)現(xiàn)方法，以對(duì)該方法有一些感性的理解。kriging模型為：
$$
Z(s) = \mu + \epsilon(s)
$$
其中，等式左側(cè)為觀測(cè)值，等式右側(cè)分別為平均值和誤差，在某一空間維度中表示如下1：

image

上面這種模型被稱為普通kriging模型，因?yàn)槲覀冊(cè)诠街幸肓顺?shù)平均值。在很多情況下，平均值函數(shù)還可以表示為一次函數(shù)、二次函數(shù)等不同種類，但是簡(jiǎn)單起見我們這里只討論平均值為常數(shù)的情況，下文中的kriging模型均指普通kriging模型。

1.2 具體方程推導(dǎo)

1.2.1 求解超參數(shù)

接下來(lái)我們推導(dǎo)一下kriging模型的具體表達(dá)形式。假設(shè)我們已知一系列觀測(cè)點(diǎn)${(\boldsymbol{x}_i,y_i)|i = 1,2,\cdots,n} $，同時(shí)我們對(duì)這些觀測(cè)點(diǎn)的取值具有一定的不確定度。為了度量這些不確定度，我們假設(shè)單個(gè)觀測(cè)點(diǎn)的分布服從高斯分布，即$y\sim N(\mu,\sigma^2)$，即我們假設(shè)這些觀測(cè)值有很大概率落在區(qū)間$[\mu-3\sigma,\mu+3\sigma]$內(nèi)?，F(xiàn)在考慮這些觀測(cè)點(diǎn)的協(xié)方差函數(shù)為$K(X,X)$，這是一個(gè)$n\times n$的矩陣。這樣我們可以得到這$n$個(gè)觀測(cè)點(diǎn)的聯(lián)合高斯分布為：
$$
\boldsymbol{y} \sim N(\mu, K(X,X))
$$
現(xiàn)在我們考慮協(xié)方差函數(shù)$K(X,X)$的具體表達(dá)形式。由統(tǒng)計(jì)知識(shí)，
$$
K(X,X) = \sigma^2\boldsymbol{R}
$$
$\boldsymbol{R}$為$n\times n $的相關(guān)系數(shù)矩陣，里面的元素為不同觀測(cè)點(diǎn)之間的相關(guān)系數(shù)，接下來(lái)我們求相關(guān)系數(shù)矩陣的具體表達(dá)形式。現(xiàn)在考慮兩個(gè)觀測(cè)點(diǎn)$\boldsymbol{x}_i$和$\boldsymbol{x}j$，如果我們的預(yù)測(cè)函數(shù)是連續(xù)的，那么當(dāng)這兩個(gè)觀測(cè)點(diǎn)之間的距離$||\boldsymbol{x}i-\boldsymbol{x}j||$之間的距離很小，這兩個(gè)觀測(cè)點(diǎn)的觀測(cè)值$y(\boldsymbol{x_i})$和$y(\boldsymbol{x_j})$也會(huì)很接近，即這兩個(gè)觀測(cè)值對(duì)應(yīng)的隨機(jī)變量之間高度相關(guān)。特別的，我們可以用高斯形式的相關(guān)系數(shù)方程來(lái)描述這樣的關(guān)系：
$$
Corr[y(\boldsymbol{x_i}),y(\boldsymbol{x_j})] = \text{exp}\bigg(-\sum{l=1}^d \theta_l|\boldsymbol{x}{il}-\boldsymbol{x}{jl}|^{{pl}}\bigg)
$$
其中，$d$為觀測(cè)點(diǎn)$\boldsymbol{x}$的維數(shù)。如果兩個(gè)觀測(cè)點(diǎn)之間的距離很近，那么它們的相關(guān)系數(shù)就會(huì)趨近于1，反之則趨近于0.兩個(gè)參數(shù)$\theta_l$和$p_l$比較關(guān)鍵，其中$\theta_l$可以理解為相關(guān)系數(shù)對(duì)第$l$個(gè)維度的坐標(biāo)距離的敏感程度，$\theta_l$越大，相關(guān)系數(shù)就對(duì)該維度的距離越敏感，下降得也越快，反之則下降得越慢。$p_l$可以衡量預(yù)測(cè)函數(shù)的光滑程度，$p_l$越接近$2$，預(yù)測(cè)函數(shù)就越光滑，可以處理一些梯度問(wèn)題，如果越接近$1$，預(yù)測(cè)函數(shù)就越粗糙，不利于處理梯度問(wèn)題。

經(jīng)過(guò)以上分析，我們可以通過(guò)$\mu,\sigma^2,\theta_l,p_l(l=1,2,\cdots,d)$這幾個(gè)參數(shù)確定這些觀測(cè)點(diǎn)的聯(lián)合高斯分布。所以現(xiàn)在需要求解這些參數(shù)來(lái)進(jìn)行之后的預(yù)測(cè)工作，我們的思路是通過(guò)求解最大似然函數(shù)來(lái)求解參數(shù)值。當(dāng)我們已知這些觀測(cè)點(diǎn)的觀測(cè)值時(shí)，意味著我們可以求出他們的似然函數(shù)表示如下：
$$
\frac{1}{(2\pi)^{n/2}(\sigma2)^{{n/2}|\boldsymbol{R}|}{1/2}}\text{exp}\bigg[\frac{-(y-\boldsymbol{1}\mu)'\boldsymbol{R}^{{-1}(y-\boldsymbol{1}\mu)}{2\sigma}2}\bigg]
$$
通過(guò)選取這些參數(shù)使得似然函數(shù)最大化意味著我們得到的這些觀測(cè)點(diǎn)和我們的預(yù)測(cè)模型符合得最好。接下來(lái)我們需要對(duì)上式取對(duì)數(shù)，然后分別對(duì)$\mu$和$\sigma^2$進(jìn)行求導(dǎo)，可以得到這兩個(gè)參數(shù)的最佳觀測(cè)值，它們都是關(guān)于相關(guān)系數(shù)矩陣$\boldsymbol{R}$的函數(shù)：
$$
-\frac{n}{2}\text{log}(\widehat{\sigma}^{2)-\frac{1}{2}\text{log}(|\boldsymbol{R}|)-\frac{(\boldsymbol{y-1}\mu)'\boldsymbol{R}}{-1}(\boldsymbol{y-1}\mu)}{2\sigma^2}\
\widehat{\mu} = \frac{\boldsymbol{1'R^{{-1}y}}{\boldsymbol{1'R}{-1}1}}\
\widehat{\sigma}^2 = \frac{(\boldsymbol{y-1}\widehat{\mu})'\boldsymbol{R}^{-1}(\boldsymbol{y-1}\widehat{\mu})}{n}
$$
將這兩個(gè)參數(shù)代入$\text{log}$形式的似然函數(shù)，得到所謂“concentrated log-likelihood function"，忽略常數(shù)項(xiàng)得到后得到：
$$
-\frac{n}{2}\text{log}(\widehat{\sigma}^2)-\frac{1}{2}\text{log}(|\boldsymbol{R}|)
$$
該方程僅依賴與相關(guān)系數(shù)矩陣$\boldsymbol{R}$，而相關(guān)系數(shù)矩陣依賴于參數(shù)$\theta_l$和$p_l$。這樣，我們的求解思路就很明確了。首先通過(guò)最優(yōu)化方法尋找使concentrated log-likelihood function 最大化的參數(shù)$\theta_l$和$p_l$，可以使用遺傳算法或粒子群算法等隨機(jī)類優(yōu)化算法避免陷入局部最優(yōu)。然后，代入求解參數(shù)$\mu$和$\sigma^{2$的估計(jì)值$\widehat{\mu}$和$\widehat{\sigma}}2$。

1.2.2 預(yù)測(cè)值與誤差估計(jì)

到此為止我們可以得到我們想要的這些參數(shù)，接下來(lái)我們分析一下預(yù)測(cè)函數(shù)的實(shí)現(xiàn)方式。簡(jiǎn)單來(lái)說(shuō)，假設(shè)我們現(xiàn)在有一個(gè)新的測(cè)量點(diǎn)$\boldsymbol{x}^*$，$y$是該測(cè)量點(diǎn)的預(yù)測(cè)值，那么根據(jù)這$(n+1)$個(gè)測(cè)量點(diǎn)我們就得到了一個(gè)新的似然方程，我們的目標(biāo)就是要尋找$y^$，使得新的似然方程最大化。求得的$y^*$也就是我們希望得到的kriging預(yù)測(cè)點(diǎn)。

假設(shè)$\widetilde{\boldsymbol{y}}=(\boldsymbol{y}'\quad y^{*)$為加入新預(yù)測(cè)值$y}$的增廣預(yù)測(cè)值向量，$\boldsymbol{r}$代表新預(yù)測(cè)值與已知預(yù)測(cè)值的相關(guān)系數(shù)向量。這樣我們得到新的$(n+1)\times(n+1)$增廣協(xié)方差矩陣$\widetilde{R}$:
$$
\widetilde{R} = \left( \begin{array}{cc}
\boldsymbol{R}&\boldsymbol{r}\
\boldsymbol{r}'&1\
\end{array}
\right)
$$
回憶之前得到的$\text{log}$形式的似然函數(shù)，如果用增廣協(xié)方差矩陣代替原來(lái)的協(xié)方差矩陣，那么該似然函數(shù)依賴于$y^$的部分為：
$$
-\frac{(\widetilde{\boldsymbol{y}}-\boldsymbol{1}\widehat{\mu})'\widetilde{\boldsymbol{R}}^{{-1}(\widetilde{\boldsymbol{y}}-\boldsymbol{1}\widehat{\mu})}{2\widehat{\sigma}}2}
$$
這一部分主要難點(diǎn)是如何求解增廣協(xié)方差矩陣的逆矩陣$\widetilde{\boldsymbol{R}}^{{-1}$，幸運(yùn)的是已經(jīng)有人幫我們求好了這一部分，由于形式較復(fù)雜，這里就不展出了。由于我們希望求解$y}$，使得增廣似然方程中依賴于$y^$的部分最大化，所以我們將上式具體展開，并對(duì)$y^$求導(dǎo)，令導(dǎo)數(shù)等于零，得到如下方程：
$$
\bigg\frac{-1}{\widehat{\sigma}^{2(1-\boldsymbol{r}'\boldsymbol{R}}{-1}\boldsymbol{r})}\bigg+\bigg[\frac{\boldsymbol{r'\boldsymbol{R}^{{-1}}(\boldsymbol{y-1}\widehat{\mu})}{\widehat{\sigma}}2(1-\boldsymbol{r}'\boldsymbol{R}^{-1}\boldsymbol{r})}\bigg]=0
$$
由上式可以求得kriging預(yù)測(cè)值$\widehat{y}(\boldsymbol{x^})$:
$$
\widehat{y}(\boldsymbol{x^{*})=\widehat{\mu}+r'\boldsymbol{R}}{-1}(\boldsymbol{y-1}\widehat{\mu})
$$
進(jìn)一步地，我們可以求出在測(cè)量點(diǎn)$\boldsymbol{x}^$處的誤差為：
$$
s^{2(\boldsymbol{x}}) = \widehat{\sigma}^{2\bigg[1-\boldsymbol{r}'\boldsymbol{R}}{-1}\boldsymbol{r} + \frac{(1-\boldsymbol{r}'\boldsymbol{R}^{{-1}\boldsymbol{r})}2}{\boldsymbol{1}'\boldsymbol{R}^{-1}\boldsymbol{1}}\bigg]
$$
該誤差方程等號(hào)右端最右側(cè)一項(xiàng)可以理解為對(duì)$\mu$預(yù)測(cè)的不確定性引起的預(yù)測(cè)值的不確定性。值得注意的是，在已知測(cè)量點(diǎn)上，該誤差為零，這一結(jié)論可以根據(jù)誤差函數(shù)嚴(yán)格證明，這里我們就不再詳述。

在我們進(jìn)行數(shù)據(jù)輸入時(shí)，為了消除單位變動(dòng)對(duì)距離的影響，我們常常將輸入數(shù)據(jù)先進(jìn)行歸一化，然后在計(jì)算結(jié)束的時(shí)候再?gòu)?fù)原數(shù)據(jù)。

2. 多精度代理模型cokriging簡(jiǎn)介

2.1 求解超參數(shù)

之前已經(jīng)在序言里闡述過(guò)，多精度代理模型的主要思想是希望能夠用計(jì)算量較低的低精度數(shù)據(jù)代替大部分需要消耗大量計(jì)算資源的高精度數(shù)據(jù)，并達(dá)到和高精度代理模型接近的計(jì)算精度，從而減少高精度數(shù)據(jù)的計(jì)算量。假設(shè)我們有兩套數(shù)據(jù)，分別是在觀測(cè)點(diǎn)$\boldsymbol{X_e}$的高精度觀測(cè)值$\boldsymbol{y}_e$和在觀測(cè)點(diǎn)$\boldsymbol{X}_c$的低精度觀測(cè)值$\boldsymbol{y}_c$，這里最好滿足$(\boldsymbol{X}_e\subset\boldsymbol{X}_c)$.我們將這些觀測(cè)點(diǎn)組合起來(lái)得到聯(lián)合觀測(cè)點(diǎn)：
$$
\boldsymbol{X} = \left(\begin{array}{c}
\boldsymbol{X}_c\
\boldsymbol{X}_e\
\end{array}\right)
=(\boldsymbol{x}_c^{1,\cdots,\boldsymbol{x}_c}{n_c},\boldsymbol{x}_e^{1,\cdots,\boldsymbol{x}_e}{n_e})^T
$$
在kriging模型中，在觀測(cè)點(diǎn)$\boldsymbol{X}$處的預(yù)測(cè)值是一個(gè)高斯隨機(jī)變量，因此這里我們有聯(lián)合預(yù)測(cè)值
$$
\boldsymbol{Y} = \left(\begin{array}{c}
\boldsymbol{Y}_c(\boldsymbol{X}_c)\
\boldsymbol{Y}_e(\boldsymbol{X}_e)\
\end{array}\right)
=(Y_c(\boldsymbol{x}_c^{1),\cdots,Y_c(\boldsymbol{x}_c}{n_c}),Y_c(\boldsymbol{x}_e^{1),\cdots,Y_c(\boldsymbol{x}_e}{n_e}))^T
$$
這里我們假設(shè)在觀測(cè)點(diǎn)$\boldsymbol{x}_i$處的預(yù)測(cè)值具有Markov性質(zhì)，即$cov{Y_e(\boldsymbol{x}^{i),Y_c(\boldsymbol{x})|Y_c(\boldsymbol{x}}i)}=0,\forall\boldsymbol{x}\ne\boldsymbol{x}^i$^[1]，這個(gè)性質(zhì)在我們求解cokriging的協(xié)方差矩陣時(shí)有很大幫助。高斯過(guò)程$Z_c(\cdot)$和$Z_e(\cdot)$代表低精度和高精度kriging模型的當(dāng)?shù)靥卣鳎鼈冎g的關(guān)系表示如下：
$$
Z_e(\boldsymbol{x}) = \rho Z_c(\boldsymbol{x}) + Z_d(\boldsymbol{x}).
$$
其中，$\rho$是一個(gè)標(biāo)量，$Z_d(\boldsymbol{x})$代表高斯過(guò)程$Z_e(\boldsymbol{x})$和$\rho Z_c(\boldsymbol{x})$之間的差別。這樣，我們就可以求出cokriging模型的協(xié)方差矩陣$cov{Y(\boldsymbol{x}),Y(\boldsymbol{x})} $：
$$
\begin{array}{l}
cov{\boldsymbol{Y}_c(\boldsymbol{X}_c),\boldsymbol{Y}_c(\boldsymbol{X}_c)} &=
cov{Z_c(\boldsymbol{X}_c),Z_c(\boldsymbol{X}_c)}\
&= \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_c,\boldsymbol{X}_c)\

cov{\boldsymbol{Y}_e(\boldsymbol{X}_e),\boldsymbol{Y}_c(\boldsymbol{X}_c)} &=
cov{\rho Z_c(\boldsymbol{X}_e) + Z_d(\boldsymbol{X}_e),Z_c(\boldsymbol{X}_c)}\
&=\rho \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_c)\

cov{\boldsymbol{Y}_e(\boldsymbol{X}_e),\boldsymbol{Y}_e(\boldsymbol{X}_e)} &=
cov{\rho Z_c(\boldsymbol{X}_e) + Z_d(\boldsymbol{X}_e),\rho Z_c(\boldsymbol{X}_e) + Z_d(\boldsymbol{X}_e)}\
&= \rho^2\sigma_c2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_e) + \sigma_d^2\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)
\end{array}
$$
這樣，我們就可以得到更為復(fù)雜的協(xié)方差矩陣：
$$
\boldsymbol{C} = \left(\begin{array}{cc}
\sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_c,\boldsymbol{X}_c) &
\rho \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_c,\boldsymbol{X}_e)\
\rho \sigma_c^2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_c)&
\rho^2\sigma_c2\boldsymbol{R}_c(\boldsymbol{X}_e,\boldsymbol{X}_e) + \sigma_d^2\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)\
\end{array}\right)
$$
從協(xié)方差矩陣的構(gòu)成元素中，我們可以清楚地看到這里我們需要求解更多的超參數(shù)，分別是：$\boldsymbol{\theta}_c,\boldsymbol{\theta}_d,\boldsymbol{p}_c,\boldsymbol{p}_d$和標(biāo)量$\rho$。由上一章節(jié)的kriging模型，我們可以根據(jù)低精度數(shù)據(jù)求出$\boldsymbol{\theta}_c$和$\boldsymbol{p}_c$，接下來(lái)我們求解剩下的參數(shù)。首先我們由高斯過(guò)程$Z_c(\cdot)$和$Z_e(\cdot)$之間的關(guān)系，定義一個(gè)新的參數(shù)：
$$
\boldsymbolu0z1t8os = \boldsymbol{y}_e - \rho\boldsymbol{y}_c(\boldsymbol{X}_e)
$$
其中$\boldsymbol{y}_c(\boldsymbol{X}_e)$在觀測(cè)點(diǎn)$\boldsymbol{x}_e$處的低精度觀測(cè)值，一般情況下我們選擇的高精度觀測(cè)點(diǎn)是低精度觀測(cè)點(diǎn)的子集，所以能夠直接得到該值。如果某些高精度觀測(cè)點(diǎn)不屬于低精度觀測(cè)點(diǎn)集，那么我們用構(gòu)造的低精度kriging模型預(yù)測(cè)值來(lái)求出$\boldsymbol{y}_c(\boldsymbol{X}_e)$。為了求出$\boldsymbol{\theta}_d,\boldsymbol{p}_d,\rho$，我們構(gòu)造出參數(shù)$\boldsymbolu0z1t8os$的$\text{log}$形式的似然函數(shù)：
$$
-\frac{n_e}{2}\text{log}(\widehat{\sigma}_d^{2)-\frac{1}{2}\text{log}(|\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)|)-\frac{(\boldsymbol{d-1}\mu_d)'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)}{-1}(\boldsymbol{d-1}\mu_d)}{2\sigma_d^2}
$$
類似地，通過(guò)求解最大似然函數(shù)，分別對(duì)$\mu_d$和$\sigma^2_d$進(jìn)行求導(dǎo)，可以求出：
$$
\begin{array}{l}
\widehat{\mu}_d &= \frac{\boldsymbol{1'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)^{{-1}d}}{\boldsymbol{1'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)}{-1}1}}\
\widehat{\sigma}^2_d &= \frac{(\boldsymbol{d-1}\widehat{\mu}_d)'\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e)^{-1}(\boldsymbol{d-1}\widehat{\mu}_d)}{n_e}
\end{array}
$$
將上式代入似然函數(shù)，可以得到新的似然函數(shù)：
$$
-\frac{n_e}{2}\text{log}(\widehat{\sigma}_d^2)-\frac{1}{2}\text{log}(|\text{det}(\boldsymbol{R}_d(\boldsymbol{X}_e,\boldsymbol{X}_e))|)
$$
通過(guò)最大化求解上面的似然函數(shù)，就可以得到我們需要的參數(shù)$\boldsymbol{\theta}_d,\boldsymbol{p}_d,\rho$，同樣我們可以使用遺傳算法或者粒子群算法等隨機(jī)類算法來(lái)進(jìn)行求解，這樣可以避免陷入局部最優(yōu)點(diǎn)。當(dāng)然如果觀測(cè)點(diǎn)的數(shù)量和維度很大，為了減少計(jì)算量我們也可以假設(shè)參數(shù)$\boldsymbol{\theta}_c,\boldsymbol{\theta}_d$均為標(biāo)量，但是這樣可能會(huì)降低模型的精確度。

2.2 預(yù)測(cè)值與誤差估計(jì)

求出超參數(shù)之后，我們可以給出cokriging模型的預(yù)測(cè)方程：
$$
\widehat{y}_e(\boldsymbol{x}^{n_e+1}) = \widehat{\mu} + \boldsymbol{c}^{T\boldsymbol{C}}{-1}(\boldsymbol{y-1}\widehat{\mu})
$$
其中，
$$
\boldsymbol{c} = \left(\begin{array}{c}
\widehat{\rho}\widehat{\sigma}_c^{2R_c(\boldsymbol{X}_c,\boldsymbol{x}}{n+1})\
\widehat{\rho}\widehat{\sigma}_c^{2R_c(\boldsymbol{X}_e,\boldsymbol{x}}{n+1}) + \widehat{\sigma}_d^{2R_d(\boldsymbol{X}_e,\boldsymbol{x}}{n+1})
\end{array}\right)
$$
該向量是新的預(yù)測(cè)點(diǎn)和已知的低精度和高精度觀測(cè)點(diǎn)的相關(guān)系數(shù)矩陣。平均值常數(shù)$ \widehat{\mu} = \frac{\boldsymbol{1'C^{{-1}y}}{\boldsymbol{1'C}{-1}1}}$ 。

進(jìn)一步地，我們還可以給出在預(yù)測(cè)點(diǎn)的誤差估計(jì)：
$$
s^2(\boldsymbol{x}) = \widehat{\rho}^{2\widehat{\sigma}_c}2 + \widehat{\sigma}_d^2 - \boldsymbol{c^TC{-1}c}
$$
如果新的預(yù)測(cè)點(diǎn)$\boldsymbol{x}^{n_e+1}\in\boldsymbol{X}_e$，那么該點(diǎn)的誤差估計(jì)為$0$，可以根據(jù)上式嚴(yán)格證明，這里我們就不再詳述。

3. 多精度代理模型示例

接下來(lái)，我們希望通過(guò)一個(gè)例子來(lái)闡述一下cokriging模型的具體實(shí)現(xiàn)效果。假設(shè)我們現(xiàn)在分別有高精度模型$f_e$和低精度模型$f_c$，表達(dá)式為：
$$
f_e(x) =(6x-2)^2\sin(12x-4),x\in[0,1]\
f_c(x) = 0.5(6x-2)^2\sin(12x-4)+10(x-0.5)+5,x\in[0,1]
$$
分別取低精度模型和高精度模型的觀測(cè)點(diǎn)$\boldsymbol{X}_c = {0.0,0.1,0.2,0.3,0.4，0.5，0.6，0.7，0.8，0.9，1.0}$和$\boldsymbol{X}_e = {0.0,0.4,0.6,1.0}$，值得注意的是高精度模型的觀測(cè)點(diǎn)包含低精度模型觀測(cè)點(diǎn)的兩個(gè)端點(diǎn)。利用這些觀測(cè)點(diǎn)和觀測(cè)值訓(xùn)練cokriging模型，然后取得剩余區(qū)間內(nèi)的高精度預(yù)測(cè)值曲線，并與真實(shí)的函數(shù)曲線進(jìn)行對(duì)比，考察cokriging模型的準(zhǔn)確度。結(jié)果如下圖所示：

圖片中綠色的點(diǎn)為低精度觀測(cè)值?？梢钥吹剑琧okriging的預(yù)測(cè)曲線與真實(shí)函數(shù)曲線幾乎重合，這說(shuō)明我們的cokriging模型效果很好，僅僅用了四個(gè)高精度觀測(cè)值就取得了很好的預(yù)測(cè)結(jié)果。

誤差曲線如下圖所示?？梢钥吹秸`差在四個(gè)高精度觀測(cè)點(diǎn)為零，最大方差在0.001左右。

這里也給出模型參數(shù)以供參考：
$$
\theta_c = 15.715,p_c =1 .9998\
\theta_d = 0.1,p_d = 1.99985,\rho=1.9836
$$

1. Kennedy M C, O'Hagan A. Predicting the output from a complex computer code when fast approximations are available[J]. Biometrika, 2000, 87(1):1-13. ?