Bases for the space of mean-value functions¶
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.
PyMimic and bases¶
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).
Built-in basis functions¶
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). Currently basis
contains only one basis, consisting of the constantly one
function. More will be added in future releases.
User-defined basis functions¶
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 = \mathbf{R}\), and the monomial basis for the space of second-degree polynomials in one variable, namely \((1, x, x^{2})\), as follows. We may construct the basis as follows.
>>> 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). Generate a sample
of this function and initialize a Blup
class.
>>> 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.
>>> ttest = np.linspace(0., 1.)
>>> blup.Beta@[fun(ttrain) for fun in basis]
array([ 4.41578258, 0.73007516, ... , 12.72544472])
References¶
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.