The Universal Kriging Model

Universal Kriging models are flexible because you can choose from a wide range of correlation functions to build the model; however, Kriging models can be computationally expensive.

The Kriging model has its roots in the field of geostatistics—a hybrid discipline of mining, engineering, geology, mathematics, and statistics (Cressie, 1993)—and is useful in predicting temporally and spatially correlated data. The Kriging model is named after D. G. Krige, a South African mining engineer who, in the 1950s, developed empirical methods for determining true ore grade distributions from distributions based on sample ore grades (Matheron, 1963). Several texts exist that describe the Kriging model and its usefulness for predicting spatially correlated data (Cressie, 1993) and mining (Journel and Huijbregts, 1978). Kriging meta models are extremely flexible because you can choose from a wide range of correlation functions to build the meta model.

In addition, depending on the correlation function that you choose, the meta model can either “honor the data,” (providing an exact interpolation of the data) or “smooth the data” (providing an inexact interpolation) (Cressie, 1993).

Depending on the number of input parameters, the number of design points, and the number of responses (outputs) of the Kriging model, the model building process can be time consuming. As the size of the matrices increases, the CPU power required for manipulating the matrices grows exponentially. Therefore, generating a good Kriging model that uses many design points can take a substantial amount of time even after all data points are analyzed.

The quality of a Kriging model depends on the location of the sample points in the design space. The Kriging model has been observed to perform best with space-filling designs where sample points are placed far apart. When points are clustered together, the matrices used in fitting the Kriging model become ill-conditioned, resulting in a poor fit. To avoid ill-conditioning, you can filter points from the sample based on a smoothing parameter, alpha, that you provide. All points that are closer than the value of alpha are removed from the sample set before fitting. Results Analytics uses other numerical techniques internally to improve the performance and robustness of the approximation.

Once a Kriging model has been built and deemed sufficiently accurate, it is ready to be used in design analysis. Fitting a Kriging model is slower than other interpolation techniques such as Radial Basis Functions; however, the prediction times are comparable.

Kriging postulates a combination of a polynomial model and departures of the following form:

where $y\left(x\right)$ is the unknown function of interest, $f\left(x\right)$ is a known polynomial function of $x$ called the trend, and $Z\left(x\right)$ is the realization of a stochastic process with mean zero, variance ${\sigma }^2$, and nonzero covariance. The $f\left(x\right)$ term in the previous equation is similar to the polynomial model in a response surface, providing a “global” model of the design space. Results Analytics allows a polynomial model in a response surface up to fourth order as the trend function $f\left(x\right)$.

While $f\left(x\right)$ “globally” approximates the design space, $Z\left(x\right)$ creates “localized” deviations so the Kriging model interpolates the $n_s$ sampled data points. The covariance matrix of $Z\left(x\right)$ that dictates the local deviations is

where $\mathbf{R}$ is the correlation matrix, and $R({\mathbf{x}}_i,{\mathbf{x}}_j)$ is the correlation function between any two of the $n_s$ sampled data points ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$. $\mathbf{R}$ is a $[n_s\times n_s]$ symmetric, positive definite matrix with ones along the diagonal.

The following correlation functions are provided by Results Analytics:

Name	Correlation
Gaussian	$${\displaystyle corr\left(X_i,X_j\right)={\displaystyle \prod }{e}^{{-\theta }_k{\left|X_{ik}-X_{jk}\right|}^2}}$$
Exponential	$${\displaystyle corr\left(X_i,X_j\right)={\displaystyle \prod }{e}^{{-\theta }_k\left|X_{ik}-X_{jk}\right|}.}$$
Cubic Spline	$$Corr(X_i,X_j)={\displaystyle \prod g}(x)$$ $$whereg(x)=\{\begin{array}{cc}1-6{({\theta }_k\left|X_{ik}-X_{jk}\right|)}^2+6{({\theta }_k\left|X_{ik}-X_{jk}\right|)}^{3}if{\theta }_k\left|X_{ik}-X_{jk}\right|<0.5& \\ 2{(1-{\theta }_k\left|X_{ik}-X_{jk}\right|)}^{3}if0.5<{\theta }_k\left|X_{ik}-X_{jk}\right|<1.0& \\ 0if{\theta }_k\left|X_{ik}-X_{jk}\right|>=1.0& \end{array}$$
Matern Linear	$${\displaystyle corr\left(X_i,X_j\right)={\displaystyle \prod \left(1+{\theta }_k\left|X_{ik}-X_{jk}\right|\right)}{e}^{{-\theta }_k\left|X_{ik}-X_{jk}\right|}}$$
Matern Cubic	$${\displaystyle corr\left(X_i,X_j\right)={\displaystyle \prod }\left(\frac{1+{\theta }_k\left|X_{ik}-X_{jk}\right|+\frac{1}{2}{\theta }_k^2\left|X_{ik}-X_{jk}\right|}{3}\right)e^{{-\theta }_k\left|X_{ik}-X_{jk}\right|}}$$

In the above table $n_{dv}$ is the number of design variables and ${\theta }_k$ are the unknown correlation parameters used to fit the model.

Once the correlation function has been selected and the best ${\theta }_k$ estimated, the Kriging model can be used to predict the response $y\left(X\right)$ at an untried location $x$ using

where $\widehat{y}$ is the vector of estimated response values at each sample point, $f$ is the vector with values of the trend function evaluated at each sample point, $\widehat{\beta }$ is a constant, and $r^T\left(X\right)$ is the vector of correlation values between the untried location $x$ and the sample data points.

The constant $\widehat{\beta }$ can be estimated using the equation

The estimate of the variance is

The maximum likelihood estimate (i.e., best) for ${\theta }_k$ is obtained by maximizing the likelihood estimate given by