What is \(\mathbf{H}(\phi)\)?

Before solving the problem, let’s clarify the notation.

\(\mathbf{H}(\phi)\) is the spatial correlation matrix of the underlying Gaussian process. It is an \(n \times n\) matrix with entries:

\[ [\mathbf{H}(\phi)]_{ij} = \rho(\phi; d_{ij}) \]

where: - \(\rho(\phi; d_{ij})\) is a correlation function - \(\phi\) is the spatial range parameter (or vector of parameters) - \(d_{ij}\) is the distance between locations \(\mathbf{s}_i\) and \(\mathbf{s}_j\)

Properties: - Diagonal elements: \(\rho(\phi; 0) = 1\) - Symmetric: \(\rho(\phi; d_{ij}) = \rho(\phi; d_{ji})\) - Positive definite (for valid correlation functions) - Usually decreases with distance

Common Correlation Functions:

Name Formula
Exponential \(\rho(\phi; d) = \exp(-d/\phi)\)
Gaussian \(\rho(\phi; d) = \exp(-d^2/\phi^2)\)
Matern \(\rho(\phi; d) = \frac{2^{1-\nu}}{\Gamma(\nu)} (d/\phi)^\nu K_\nu(d/\phi)\)
Spherical \(\rho(\phi; d) = 1 - \frac{3}{2}(d/\phi) + \frac{1}{2}(d/\phi)^3\) for \(d \leq \phi\); 0 for \(d > \phi\)

Role of \(\mathbf{H}(\phi)\) in the Model

In our model:

\[ \boldsymbol{\Sigma} = \sigma^2 \mathbf{H}(\phi) + \tau^2 \mathbf{I} \]

Components:

  1. \(\sigma^2 \mathbf{H}(\phi)\): The spatial covariance component
    • \(\sigma^2\) is the partial sill (spatial variance)
    • \(\mathbf{H}(\phi)\) captures the correlation structure between locations
    • This represents the spatial dependence in the underlying process
  2. \(\tau^2 \mathbf{I}\): The nugget effect (measurement error or microscale variation)
    • \(\tau^2\) is the nugget variance
    • \(\mathbf{I}\) is the identity matrix
    • This represents independent noise at each location

Interpretation:

  • When \(\tau^2 = 0\): No nugget effect, the model interpolates exactly through data points
  • When \(\tau^2 > 0\): There is measurement error or microscale variation, the model smooths rather than interpolates

Example: Constructing \(\mathbf{H}(\phi)\)

Suppose we have 3 locations, and we use an exponential correlation function with \(\phi = 1\) and \(\sigma^2 = 2\):

The correlation function is: \[ \rho(\phi; d) = \exp(-d/\phi) = \exp(-d) \]

The correlation matrix is:

\[ \mathbf{H}(\phi) = \begin{bmatrix} 1 & e^{-d_{12}} & e^{-d_{13}} \\ e^{-d_{21}} & 1 & e^{-d_{23}} \\ e^{-d_{31}} & e^{-d_{32}} & 1 \end{bmatrix} \]

If the distances are: - \(d_{12} = 0.5\) - \(d_{13} = 1.0\) - \(d_{23} = 0.7\)

Then:

\[ \mathbf{H}(\phi) = \begin{bmatrix} 1 & e^{-0.5} & e^{-1.0} \\ e^{-0.5} & 1 & e^{-0.7} \\ e^{-1.0} & e^{-0.7} & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0.6065 & 0.3679 \\ 0.6065 & 1 & 0.4966 \\ 0.3679 & 0.4966 & 1 \end{bmatrix} \]

Then the covariance matrix of the observed data is:

\[ \boldsymbol{\Sigma} = 2\mathbf{H}(\phi) + \tau^2 \mathbf{I} \]

If \(\tau^2 = 0\) (noiseless case):

\[ \boldsymbol{\Sigma} = \begin{bmatrix} 2 & 1.213 & 0.7358 \\ 1.213 & 2 & 0.9932 \\ 0.7358 & 0.9932 & 2 \end{bmatrix} \]

If \(\tau^2 = 0.5\) (with noise):

\[ \boldsymbol{\Sigma} = \begin{bmatrix} 2.5 & 1.213 & 0.7358 \\ 1.213 & 2.5 & 0.9932 \\ 0.7358 & 0.9932 & 2.5 \end{bmatrix} \]

Notice that when \(\tau^2 > 0\), the diagonal elements increase (representing additional variance from noise), but the off-diagonal elements remain the same (representing spatial correlation).

Problem Statement

We are given the joint normal distribution:

\[ \begin{pmatrix} Y_1 \\ Y_2 \end{pmatrix} = \begin{pmatrix} Y(\mathbf{s}_0) \\ \mathbf{y} \end{pmatrix} \sim N\left( \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix}, \begin{pmatrix} \Omega_{11} & \Omega_{12} \\ \Omega_{21} & \Omega_{22} \end{pmatrix} \right) \]

with:

  • \(\mu_1 = E[Y(\mathbf{s}_0)] = \mathbf{x}_0^T \boldsymbol{\beta}\)
  • \(\mu_2 = E[\mathbf{y}] = \mathbf{X} \boldsymbol{\beta}\)
  • \(\Omega_{11} = \text{Var}(Y(\mathbf{s}_0)) = \sigma^2 + \tau^2\)
  • \(\Omega_{12} = \text{Cov}(Y(\mathbf{s}_0), \mathbf{y}) = \boldsymbol{\gamma}^T\)
  • \(\Omega_{21} = \text{Cov}(\mathbf{y}, Y(\mathbf{s}_0)) = \boldsymbol{\gamma}\)
  • \(\Omega_{22} = \text{Var}(\mathbf{y}) = \boldsymbol{\Sigma} = \sigma^2 \mathbf{H}(\phi) + \tau^2 \mathbf{I}\)

We need to verify expressions (2.22), (2.23), and (2.24).

Step 1: Verify Expression (2.22) - Conditional Mean

For a joint normal distribution, the conditional mean is:

\[ E[Y_1 | Y_2 = \mathbf{y}] = \mu_1 + \Omega_{12} \Omega_{22}^{-1} (\mathbf{y} - \mu_2) \]

Substituting the given values:

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \mathbf{x}_0^T \boldsymbol{\beta} + \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \]

Therefore:

\[ \boxed{E[Y(\mathbf{s}_0) | \mathbf{y}] = \mathbf{x}_0^T \boldsymbol{\beta} + \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})} \]

This matches (2.22) exactly. ✓

Step 2: Verify Expression (2.23) - Conditional Variance

For a joint normal distribution, the conditional variance is:

\[ \text{Var}[Y_1 | Y_2 = \mathbf{y}] = \Omega_{11} - \Omega_{12} \Omega_{22}^{-1} \Omega_{21} \]

Substituting the given values:

\[ \text{Var}[Y(\mathbf{s}_0) | \mathbf{y}] = \sigma^2 + \tau^2 - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma} \]

Therefore:

\[ \boxed{\text{Var}[Y(\mathbf{s}_0) | \mathbf{y}] = \sigma^2 + \tau^2 - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma}} \]

This matches (2.23) exactly. ✓

Step 3: Verify Expression (2.24) - Kriging Weights \(\lambda\)

We want to express the conditional mean as a linear combination of \(\mathbf{y}\):

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \lambda^T \mathbf{y} \]

Step 3.1: Use the GLS Estimator

The Generalized Least Squares (GLS) estimator of \(\boldsymbol{\beta}\) is:

\[ \hat{\boldsymbol{\beta}}_{\text{GLS}} = \left( \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{y} \]

Step 3.2: Rewrite the Conditional Mean

Starting from (2.22):

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \mathbf{x}_0^T \boldsymbol{\beta} + \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \]

Replace \(\boldsymbol{\beta}\) with \(\hat{\boldsymbol{\beta}}_{\text{GLS}}\):

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \mathbf{x}_0^T \hat{\boldsymbol{\beta}}_{\text{GLS}} + \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \left( \mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}}_{\text{GLS}} \right) \]

Step 3.3: Expand and Group Terms

Expanding the second term:

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \mathbf{x}_0^T \hat{\boldsymbol{\beta}}_{\text{GLS}} + \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{y} - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}}_{\text{GLS}} \]

Grouping terms involving \(\hat{\boldsymbol{\beta}}_{\text{GLS}}\):

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{y} + \left( \mathbf{x}_0^T - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right) \hat{\boldsymbol{\beta}}_{\text{GLS}} \]

This gives us:

\[ \boxed{E[Y(\mathbf{s}_0) | \mathbf{y}] = \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{y} + \left( \mathbf{x}_0^T - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right) \hat{\boldsymbol{\beta}}_{\text{GLS}}} \]

Note: The key step is factoring out \(\hat{\boldsymbol{\beta}}_{\text{GLS}}\) from \(\mathbf{x}_0^T \hat{\boldsymbol{\beta}}_{\text{GLS}} - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}}_{\text{GLS}}\) to get \((\mathbf{x}_0^T - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X}) \hat{\boldsymbol{\beta}}_{\text{GLS}}\).

Step 3.4: Substitute \(\hat{\boldsymbol{\beta}}_{\text{GLS}}\)

Substitute the GLS estimator into the expression:

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{y} + \left( \mathbf{x}_0^T - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right) \left( \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{y} \]

Step 3.5: Factor Out \(\mathbf{y}\)

Factor \(\mathbf{y}\) from both terms:

\[ E[Y(\mathbf{s}_0) | \mathbf{y}] = \left[ \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} + \left( \mathbf{x}_0^T - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right) \left( \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \right] \mathbf{y} \]

Step 3.6: Identify \(\lambda^T\)

Since \(E[Y(\mathbf{s}_0) | \mathbf{y}] = \lambda^T \mathbf{y}\), we have:

\[ \lambda^T = \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} + \left( \mathbf{x}_0^T - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right) \left( \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right)^{-1} \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \]

Step 3.7: Take the Transpose

Taking the transpose of both sides:

\[ \lambda = \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma} + \boldsymbol{\Sigma}^{-1} \mathbf{X} \left( \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right)^{-1} \left( \mathbf{x}_0 - \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma} \right) \]

Therefore:

\[ \boxed{\lambda = \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma} + \boldsymbol{\Sigma}^{-1} \mathbf{X} \left( \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X} \right)^{-1} \left( \mathbf{x}_0 - \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma} \right)} \]

This matches (2.24) exactly. ✓

Summary of Verification

Expression Result Status
(2.22) \(E[Y(\mathbf{s}_0) \| \mathbf{y}] = \mathbf{x}_0^T \boldsymbol{\beta} + \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})\) ✓ Verified
(2.23) \(\text{Var}[Y(\mathbf{s}_0) \| \mathbf{y}] = \sigma^2 + \tau^2 - \boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma}\) ✓ Verified
(2.24) \(\lambda = \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma} + \boldsymbol{\Sigma}^{-1} \mathbf{X} (\mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} (\mathbf{x}_0 - \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma})\) ✓ Verified

Key Insights

  1. Conditional Mean: The conditional mean of \(Y(\mathbf{s}_0)\) given observed data \(\mathbf{y}\) consists of:

    • A trend component: \(\mathbf{x}_0^T \boldsymbol{\beta}\)
    • A correction term: \(\boldsymbol{\gamma}^T \boldsymbol{\Sigma}^{-1} (\mathbf{y} - \mathbf{X} \boldsymbol{\beta})\)

  2. Conditional Variance: The conditional variance is the unconditional variance minus the variance explained by the observed data.

  3. Kriging Weights: The weight vector \(\lambda\) accounts for both:

    • Spatial correlation: \(\boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma}\)
    • Trend adjustment: \(\boldsymbol{\Sigma}^{-1} \mathbf{X} (\mathbf{X}^T \boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} (\mathbf{x}_0 - \mathbf{X}^T \boldsymbol{\Sigma}^{-1} \boldsymbol{\gamma})\)

These results form the foundation of universal kriging with measurement error.