Kriging Weights

1. template<
            class MatrixLibrary,
            class Location,
            class Neighborhood,
            class Covariance,
            class KrigingConstraints,
            class Vector
           >
   int
   kriging_weights(Vector& weights, double& kriging_variance,
                   const Location& center, const Neigborhood& neighbors,
                   Covariance& covar, KrigingConstraints& Kconstraint);
   
   
   
2. template<
            class MatrixLibrary,
            class Location,
            class Neighborhood,
            class Covariance,
            class KrigingConstraints,
            class Vector
           >
   int
   kriging_weights(Vector& weights, 
                   const Location& center, const Neigborhood& neighbors,
                   Covariance& covar, KrigingConstraints& Kconstraint);
   
   
   

The kriging_weights algorithm solves a kriging system :

$$\left\{\vphantom{ \begin{array}{l} Var \Big( \sum_{\alpha=1}^{n}... ... \vdots \\ f_p \Big( \{ \lambda_{\alpha} \} \Big)=0 \end{array} }\right.$$$$\begin{array}{l} Var \Big( \sum_{\alpha=1}^{n} \lambda_{\alpha} ... ...Big)=0 \\ \vdots \\ f_p \Big( \{ \lambda_{\alpha} \} \Big)=0 \end{array}$$

where f_i$\Big(${$\lambda_{\alpha}^{}$}$\Big)$ = 0, i = 1,..., p are linear constraints expressed by the Kriging Constraint Kconstraint. The solution to this system is a set of weights:

($$\lambda_{1}^{}$$,...,$$\lambda_{n}^{}$$,$$\mu_{1}^{}$$,...,$$\mu_{p}^{}$$)

where the weights $\mu_{j}^{}$ are the Lagrange weights used to account for constraints f₁,..., f_p.

The kriging_weights function computes the kriging weights and stores them into Vector weights in the following order: ($\lambda_{1}^{}$,...,$\lambda_{n}^{}$,$\mu_{1}^{}$,...,$\mu_{p}^{}$). Only the weights $\lambda_{1}^{}$,...,$\lambda_{n}^{}$ are useful to compute the kriging estimate. The weights $\mu_{1}^{}$,...,$\mu_{p}^{}$ are used to compute the kriging variance. Vector weights does not need to be of the correct size n + p when passed to the function: the function automatically resizes the vector if necessary. Two versions of the algorithm exist. Version 1 computes the kriging variance and stores it into kriging_variance, while Version 2 does not compute the kriging variance (there are cases where the kriging variance is useless, e.g. in indicator kriging). Both version of the algorithm allow to change the linear algebra library (which defines matrices, matrix inversion routines, ...) used by kriging_weights (see 2.4). The default linear algebra library is the TNTlibrary (slightly modified), a public domain library by Roldan Pozo, Mathematical and Computational Sciences Division,National Institute of Standards and Technology.

The value returned by kriging_weights indicates if any problem were encountered:

• returns 0 if the function executed successfully
• returns 1 if the kriging system could not be solved
• returns 2 if the neighborhood of conditioning data was empty

Where defined

In header file <kriging.h>

Requirements on types

• Vector is a container on which an iterator is defined. It must have three member functions whose prototypes are:
  1. iterator Vector::begin() which returns an iterator to the first element of the Vector
  2. int Vector::size() which returns the size of the Vector.
  3. void Vector::resize(int n) which changes the size of the Vector to n
• Location is a model of Location.
• Neighborhood is a model of Neighborhood. The location type of the geo-values contained in the neighborhood must be location.
• Covariance is a model of Covariance that takes two objects of type location as arguments.
• KrigingConstraints is a model of Kriging Constraints.
• MatrixLibrary specifies the linear algebra library to use. The requirements on matrix_lib are fully defined in 2.4. By default, it is the TNTlibrary (slightly modified).

Remarks

• Vector weights is resized after it is passed to kriging_weights, unless it is of the correct size. The cost of this resize can be decreased by smartly managing the Vector's memory (in the same style as an STL vector): the vector allocates more memory than it needs, hence does not need to re-allocate memory each time its size increases.

  However, each call to kriging_weights implies solving a linear system of equations, hence the cost of the resize would usually be negligible.

• If the same data z($\bf u_{1}^{}$),..., z($\bf u_{d}^{}$) are used for estimating different locations, algorithm global_neigh_kriging_weights can be more suitable than kriging_weights. The kriging matrix does not depend on the location to be estimated, but only on the data locations $\bf u_{1}^{}$,...,$\bf u_{d}^{}$. Hence if the exact same data are used to estimate multiple locations, the kriging matrix remains unchanged: it can be inverted once for all and solving further kriging systems reduces to a mere matrix-vector multiplication. global_neigh_kriging_weights implements such ``global kriging''.

Example

Estimate z(0, 0) from z(2, 3) = 0.21 and z(4, - 7) = 0.09 by ordinary kriging. After the call to kriging_weight, vector weights has size 3; its two first elements are the weights corresponding to z(2, 3) and z(4, - 7) respectively, and its last element is the Lagrange parameter.

typedef std::vector<geo_value2d> neighborhood;

// gaussian covariance of range 4
inline double gauss_covariance(Location2d u1, Location2d u2){
  double square_dist = pow(u1[0]-u2[0], 2) +
                       pow(u1[1]-u2.[1], 2);
  return exp(-3*square_dist/16);                     
}


int main()
{  
  location2d u1(2,3);
  location2d u2(4,-7);

  geo_value2d Z1(u1,0.21);
  geo_value2d Z2(u2,0.09); 
  
  Neighborhood neighbors;
  neighbors.push_back(Z1);
  neighbors.push_back(Z2);

  location2d u(0,0);

  std::vector<double> weights;

  double kvariance;
  kriging_weights(weights, kvariance,
                  u, &neighbors,
                  gauss_covariance, OK_constraint);
}