Kriging模型小结

啰嗦一下

最近看很多文章，在使用回归模型时无外乎以下几种，Kriging、RBFN、高斯过程模型、低阶多项式多重回归模型，其中Kriging被提及最多，以及它有个性质，关于Kriging的不确定性信息似乎有很多故事…所以看来已经到了不得不弄懂的时候了。

Kriging模型是啥

kriging模型是一个插值方法，啥意思呢？就是假设我们已知有一群样本点，它们的分布符合高斯分布，通过这个模型，可以向这些已知点中插入一些预测点，这些点也是和原始的点服从同样的分布的。
为什么它可以作为回归模型？
我们看这样一个式子

y=∑i=0nβixi+b(1)y=\sum _{i=0}^{n} \beta _{i}x_{i}+b\tag{1}

我们并不知道上式中的每一项的系数βi\beta_i和常数项b的具体值，但是我们有一些已知的点(x,y)(\bm {x},y)，我们希望通过这些已知的点来求出这些未知的参数，将式（1）变成一个确定的关系式，这样，我们在得到一个陌生的点xj时\bm {x_j}时，可以根据我们拟合出的确定的式子得到y，就达到了预测的目的。
那一般如何预测这些参数呢？通过最小化预测点到已知点之间的均方误差，也就是我们所说的最小二乘法，得到的参数就可以将（1）确定下来。
所以，如果输入一个未知点，我们通过模型输出的是一个连续的值，那么这样的模型都可以称为回归，自然，kriging也符合回归的定义。

小背景

原来kriging最初是用在地质统计学领域的，当时要挖矿嘛，人们希望通过钻的有限个孔当中有无矿的情况来推测没有钻孔的陌生土地上哪些地方有矿的可能性比较大，进而高效挖矿（物理）。动词是krige，但是更为广泛流传的是其名词形式Kriging，这个k大小写都可以。

Kriging模型建立

说完了简介和定义，大概知道了这个模型是干嘛了，那么如何建立这个模型，其实就是要确定其中的一些参数嘛，那参数怎么确定，不就是要通过优化某个东西使其达到最优嘛。
这里使用的就是最大似然估计。
（ps：这里需要回顾概率统计中的参数估计知识，可以查阅课本或是google一下。贴心地放上一个链接）
最大似然估计介绍
1、多元高斯分布
一维的高斯分布我们很容易知道它的图像长什么样，待确定的参数就是均值μ\mu和方差σ\sigma。但是对于给定的数据集来说，每个x\bm {x}可能是大于一维的，这时候每个样本服从的分布就称为多元高斯分布。这时候需要确定的参数就变成了均值μ\mu和协方差CovCov。
对于给定的数据集X={x1,x2,...,xn}X=\{\bm{x_1},\bm {x_2},...,\bm {x_n}\}，其对应的目标值为Y={y1,y2,...,yn}Y=\{\bm{y_1},\bm {y_2},...,\bm {y_n}\}，有以下定义

Cov=(cor(Y(x1),Y(x1)),…,cor(Y(x1),Y(xn))⋮,⋱,⋮cor(Y(xn),Y(x1)),…,cor(Y(xn),Y(xn)))Cov=\begin{pmatrix} cor(Y(\bm{x_1}),Y(\bm{x_1})),\dots,cor(Y(\bm{x_1}),Y(\bm{x_n})) \\ \vdots, \ddots,\vdots \\ cor(Y(\bm{x_n}),Y(\bm{x_1})),\dots,cor(Y(\bm{x_n}),Y(\bm{x_n})) \end{pmatrix}

注意，这个协方差矩阵主对角线上的值是每个样本的方差

cor[Y(xi),Y(xl)]=exp(−∑j=1kθj∣xij−xlj∣2)cor[Y(\bm{x}_i),Y(\bm{x}_l)]=exp(-\sum_{j=1}^{k}{\theta_j|x^j_i-x^j_l|^2})

其中j表示每个变量在第j维的值

2、最大似然估计
既然已经确定了这个多元高斯分布的参数定义，我们需要确定的先验参数有μ，σ和θ\mu，\sigma和\theta，服从正态分布的似然函数的定义，

L(Y1,Y2,…,Yn∣μ,σ)=1(2πσ2)n/2exp(−Σ(Yi−μ)22σ2)L(\bm{Y}_1,\bm{Y}_2,\dots,\bm{Y}_n|\mu,\sigma)=\frac{1}{(2\pi\sigma^2)^{n/2}}exp(-\frac{\Sigma(\bm{Y_i}-\bm{\mu})^2}{2\sigma^2})

结合上述的参数表示我们可以得到如下的表示

L=1(2πσ2)n/2∣Cov∣1/2exp[−(y−1μ)TCov−1(y−1μ)2σ2]L=\frac{1}{(2\pi\sigma^2)^{n/2}|Cov|^{1/2}}exp[-\frac{(\bm{y}-\bm{1}\mu)^T{Cov}^{-1}(\bm{y}-\bm{1}\mu)}{2\sigma^2}]

一般我们会取个对数，

ln(L)=−n2ln(2π)−n2ln(σ2)−12ln∣Cov∣−(y−1μ)TC−1(y−1μ)2σ2(2)ln(L)=-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^2)-\frac{1}{2}ln|Cov|-\frac{(\bm{y}-\bm{1}\mu)^TC^{-1}(\bm{y}-\bm{1}\mu)}{2\sigma^2} \tag{2}

好啦，现在就是要求使得式（2）最大的μ\mu和σ\sigma的值，就是求偏导啦。
具体过程涉及矩阵的计算，详细这里就不写啦，可以参考这篇详细推导过程
经历了上述计算之后，可以得到

lnL≈−n2ln(σ^2)−12ln(∣Cov∣)(3)lnL\approx -\frac{n}{2}ln(\hat{\sigma}^2)-\frac{1}{2}ln(|Cov|)\tag{3}

其中，

μ^=1TCov−1y1TCov−11\hat{\mu}=\frac{\bm{1}^T{Cov}^{-1}y}{\bm{1}^T{Cov}^{-1}\bm{1}}

σ^2=(y−1μ^)TCov−1(y−1μ^)n\hat{\sigma}^2=\frac{(y-\bm{1}\hat{\mu})^T{Cov}^{-1}(y-\bm{1}\hat{\mu})}{n}

至此，我们只需要求出可以使得式(3)最大的θ\theta，就完全确定了参数啦，就可以进行对未知点的预测了。

3、模型预测
对于一个未知的点xx，通过kriging模型，我们就可以得到

{y^1(x)=μ^+rTCov−1(y−1μ^)s12(x)=σ^[1−rTCov−1r+(1−1TCov−1r)21TCov−11]\left\{ \begin{array}{l} \hat{y}_{1}(x)=\hat{\mu}+r^T{Cov}^{-1}(y-\bm{1}\hat{\mu}) \\[3mm]s_{1}^2(x)=\hat{\sigma}[1-r^T{Cov}^{-1}r+\frac{(1-\bm{1}^T{Cov}^{-1}r)^2}{\bm{1}^T{Cov}^{-1}\bm{1}}] \end{array} \right.

{μ^=1TCov−1y1TCov−11σ^2=(y−1μ^)TCov−1(y−1μ^)n\left\{ \begin{array}{l} \hat{\mu}=\frac{\bm{1}^T{Cov}^{-1}y}{\bm{1}^T{Cov}^{-1}\bm{1}} \\[3mm]\hat{\sigma}^2=\frac{(y-\bm{1}\hat{\mu})^T{Cov}^{-1}(y-\bm{1}\hat{\mu})}{n} \end{array} \right.

这里得到的y^1(x)\hat{y}_{1}(x)和s12(x)s_{1}^2(x)就是对未知点xx预测得到的正态分布的均值和方差，所以说这个模型会有不确定性。
r是数据点和待预测点之间的协方差矩阵。