When computing the BLUP we must specify a basis for the space of mean-value functions. In ordinary kriging the random process (Z, X) is assumed to be weakly stationary, and hence to have constant mean. In universal kriging the random process (Z, X) is assumed to be intrinsically stationary with polynomial mean.
Be aware that a basis of size comparable to the size of X will cause the BLUP to over-fit the data.
The basis must be specified using the keyword argument basis
when creating an instance of the Blup
. It cannot be changed afterwards. The basis must be a tuple of functions func(x) -> array_like
with signature (n, d) -> (n,)
, where d
is the dimension of the index set T (The two flavours of linear
prediction <the_two_flavours_of_linear_prediction>
).
The module basis
provides a suite of bases, which may be used when specifying the BLUP using the class Blup
. The available bases are described in the submodule's API documentation (pymimic.basis module <basis_module>
). Currently basis
contains only one basis, consisting of the constantly one function. More will be added in future releases.
Because Numpy and Scipy functions are vectorized they naturally have the required signature. It is therefore convenient to use Numpy and Scipy functions directly, or to construct basis functions using them.
Consider the case of one-dimensional index set T = R, and the monomial basis for the space of second-degree polynomials in one variable, namely (1, x, x2), as follows. We may construct the basis as follows.
python
>>> import numpy as np >>> basis = (np.polynomial.Polynomial([1.]), np.polynomial.Polynomial([0., 1.]), np.polynomial.Polynomial([0., 0., 1.]))
Once we have initiated a Blup
class the coefficients of the best linear unbiased estimator of the mean are stored as the attribute Beta
. Let us return to the example of the Forrester function (Quick-start tutorial <quickstart_tutorial>
). Generate a sample of this function and initialize a Blup
class.
python
>>> ttrain = np.linspace(0., 1., 10) >>> xtrain = mim.testfunc.forrester(ttrain) + 0.5*np.random.randn(10) >>> blup = mim.Blup(ttrain, xtrain, 0.5**2., basis=basis) >>> blup.opt() direc: array([[-0.01099084, 0.00122937], [ 0.00069455, -0.00117173]]) fun: 29.366313678967565 message: 'Optimization terminated successfully.' nfev: 418 nit: 11 status: 0 success: True x: array([1.49828726, 0.738667 ])
As well as computing the BLUP, using xtest
, we may compute the best linear unbiased estimators of the mean.
python
>>> ttest = np.linspace(0., 1.) >>> blup.Beta@[fun(ttrain) for fun in basis] array([ 4.41578258, 0.73007516, ... , 12.72544472])
- C86
Cressie, N. 1986. 'Kriging nonstationary data' in Journal of the American Statistical Association, 81 (395), 625--34. Available at https://doi.org/10.2307/2288990.