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

[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.