1. The Matérn Covariance Function

The Matérn family is a flexible class of covariance functions for spatial data. The full covariance function is:

\[ C(d \mid \psi) = \sigma^2 \cdot \rho(d \mid \upsilon) \]

where the correlation function \(\rho(d \mid \upsilon)\) is given by:

\[ \rho(d \mid \phi, \nu) = \frac{1}{2^{\nu-1}\Gamma(\nu)} \left( \sqrt{2\nu}\, \phi d \right)^\nu K_\nu\!\left( \sqrt{2\nu}\, \phi d \right) \tag{1} \]

and:

  • \(d \ge 0\): Euclidean distance between two locations
  • \(\phi > 0\): Range parameter (controls decay speed)
  • \(\nu > 0\): Smoothness parameter (controls differentiability)
  • \(\Gamma(\nu)\): Gamma function
  • \(K_\nu(\cdot)\): Modified Bessel function of the second kind, order \(\nu\)
  • \(\sigma^2 > 0\): Marginal variance (sill)

The full parameter vector is \(\psi = (\sigma^2, \phi, \nu)\), with \(\upsilon = (\phi, \nu)\).

2. What This Equation Represents

For two spatial locations \(s_i\) and \(s_j\) separated by distance \(d = \|s_i - s_j\|\), the covariance between the random variables \(Z(s_i)\) and \(Z(s_j)\) is:

\[ \operatorname{Cov}\big(Z(s_i), Z(s_j)\big) = C(d \mid \psi) \]

The equation defines each entry of the covariance matrix, not the matrix itself.

3. Computing the Matérn in R

R provides the Bessel function \(K_\nu\) as besselK(x, nu). Below is an implementation of the Matérn correlation function.

# Matérn correlation function (equation 2.2)
matern_correlation <- function(d, phi, nu) {
  # d: distance (scalar or vector)
  # phi: range parameter (>0)
  # nu: smoothness parameter (>0)
  
  # Handle d = 0 case (correlation = 1)
  if (length(d) == 1 && d == 0) return(1)
  if (any(d == 0)) {
    result <- sapply(d, function(di) {
      if (di == 0) return(1)
      x <- sqrt(2 * nu) * phi * di
      (1 / (2^(nu - 1) * gamma(nu))) * (x^nu) * besselK(x, nu)
    })
    return(result)
  }
  
  x <- sqrt(2 * nu) * phi * d
  rho <- (1 / (2^(nu - 1) * gamma(nu))) * (x^nu) * besselK(x, nu)
  return(rho)
}

# Full covariance function
matern_covariance <- function(d, sigma2, phi, nu) {
  return(sigma2 * matern_correlation(d, phi, nu))
}

4. Example: Fixed Parameters

Let’s use the example: \(\phi\)=0.3, \(\nu\)=0.3, \(\sigma^2\)=2.0

sigma2 <- 2.0
phi <- 0.3
nu <- 1.5

# Create a sequence of distances
distances <- seq(0, 15, length.out = 200)

# Compute correlation and covariance
correlation <- matern_correlation(distances, phi, nu)
covariance <- matern_covariance(distances, sigma2, phi, nu)

# Create data frame for plotting
plot_df <- data.frame(
  distance = distances,
  correlation = correlation,
  covariance = covariance
)

4.1 Correlation vs Distance

library(ggplot2)
ggplot(plot_df, aes(x = distance, y = correlation)) +
  geom_line(color = "steelblue", size = 1.2) +
  labs(
    title = paste0("Matérn Correlation Function (", 
                   "φ = ", phi, ", ν = ", nu, ")"),
    x = "Distance d",
    y = expression(rho(d))
  ) +
  theme_minimal() +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "red") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  annotate("text", x = 12, y = 0.08, label = "Practical range (ρ = 0.05)", 
           color = "red", size = 3)

4.2 Covariance vs Distance

ggplot(plot_df, aes(x = distance, y = covariance)) +
  geom_line(color = "darkgreen", size = 1.2) +
  labs(
    title = paste0("Matérn Covariance Function (",
                   "σ² = ", sigma2, ", φ = ", phi, ", ν = ", nu, ")"),
    x = "Distance d",
    y = expression(C(d))
  ) +
  theme_minimal() +
  geom_hline(yintercept = 0, linetype = "dotted")

5 Effect of varying parameters

5.1 Varying \(\phi\) (range parameter).

Large \(\phi\) means faster decay and shorter range dependence.

phi_values <- c(0.2, 0.5, 1.0)
nu_fixed <- 1.5

phi_df <- data.frame()
for (phi_val in phi_values) {
  corr <- matern_correlation(distances, phi_val, nu_fixed)
  phi_df <- rbind(phi_df, data.frame(
    distance = distances,
    correlation = corr,
    phi = as.factor(phi_val)
  ))
}

ggplot(phi_df, aes(x = distance, y = correlation, color = phi)) +
  geom_line(size = 1.1) +
  labs(
    title = paste0("Effect of φ on Matérn Correlation (ν = ", nu_fixed, ")"),
    x = "Distance d",
    y = "Correlation ρ(d)",
    color = expression(phi)
  ) +
  theme_minimal() +
  scale_color_brewer(palette = "Set1")

5.2 Varying \(\nu\) (Smoothness parameter)

\(\nu\)=0.5: Exponential model (rough, not differentiable)

\(\nu\)=1.5: Smooth, once differentiable

\(\nu\) =2.0: Very smooth, twice differentiable

nu_values <- c(0.5, 1.0, 1.5, 2.0)
phi_fixed <- 0.3

nu_df <- data.frame()
for (nu_val in nu_values) {
  corr <- matern_correlation(distances, phi_fixed, nu_val)
  nu_df <- rbind(nu_df, data.frame(
    distance = distances,
    correlation = corr,
    nu = as.factor(nu_val)
  ))
}

ggplot(nu_df, aes(x = distance, y = correlation, color = nu)) +
  geom_line(size = 1.1) +
  labs(
    title = paste0("Effect of ν on Matérn Correlation (φ = ", phi_fixed, ")"),
    x = "Distance d",
    y = "Correlation ρ(d)",
    color = expression(nu)
  ) +
  theme_minimal() +
  scale_color_brewer(palette = "Set2")

6. Practical Range

The practical range is the distance where correlation drops to 0.05. For the Matérn, it approximates:

Practical Range \(\approx \frac{\sqrt{(8 \nu)}}{\phi}\)

nu_example <- 1.5
phi_example <- 0.3
practical_range <- sqrt(8 * nu_example) / phi_example

cat("For φ =", phi_example, "and ν =", nu_example, "\n")
## For φ = 0.3 and ν = 1.5
cat("Practical range ≈", round(practical_range, 2), "units\n")
## Practical range ≈ 11.55 units
# Verify by finding distance where correlation ≈ 0.05
corr_at_range <- matern_correlation(practical_range, phi_example, nu_example)
cat("Correlation at that distance =", round(corr_at_range, 4), "\n")
## Correlation at that distance = 0.0174

7. Building a Covariance Matrix

For \(n\) spatial locations, the covariance matrix \(\Sigma\) is \(n \times n\) with entries:

\[ \Sigma_{ij} = C(\|s_i - s_j\| \mid \sigma^2, \phi, \nu) \]

# Example: 4 random locations
set.seed(123)
locations <- matrix(runif(8, 0, 10), ncol = 2)  # 4 locations in 2D
colnames(locations) <- c("x", "y")
rownames(locations) <- paste0("s", 1:4)

# Compute distance matrix
dist_matrix <- as.matrix(dist(locations))

# Compute covariance matrix
n <- nrow(locations)
cov_matrix <- matrix(0, n, n)
for (i in 1:n) {
  for (j in 1:n) {
    cov_matrix[i, j] <- matern_covariance(dist_matrix[i, j], sigma2, phi, nu)
  }
}

cat("Locations:\n")
## Locations:
print(round(locations, 2))
##       x    y
## s1 2.88 9.40
## s2 7.88 0.46
## s3 4.09 5.28
## s4 8.83 8.92
cat("\nDistance matrix:\n")
## 
## Distance matrix:
print(round(dist_matrix, 2))
##       s1    s2   s3   s4
## s1  0.00 10.25 4.30 5.97
## s2 10.25  0.00 6.14 8.52
## s3  4.30  6.14 0.00 5.98
## s4  5.97  8.52 5.98 0.00
cat("\nCovariance matrix (σ² =", sigma2, ", φ =", phi, ", ν =", nu, "):\n")
## 
## Covariance matrix (σ² = 2 , φ = 0.3 , ν = 1.5 ):
print(round(cov_matrix, 3))
##       [,1]  [,2]  [,3]  [,4]
## [1,] 2.000 0.061 0.693 0.368
## [2,] 0.061 2.000 0.345 0.130
## [3,] 0.693 0.345 2.000 0.368
## [4,] 0.368 0.130 0.368 2.000

8. Verifying Positive Definiteness

A valid covariance matrix must be positive definite. Let’s check the eigenvalues

eigenvals <- eigen(cov_matrix, symmetric = TRUE, only.values = TRUE)$values
cat("Eigenvalues:\n", round(eigenvals, 4), "\n")
## Eigenvalues:
##  3.0664 1.9648 1.7198 1.2489
cat("All eigenvalues > 0?", all(eigenvals > 1e-10), "\n")
## All eigenvalues > 0? TRUE

9. Summary

Summary of the Matérn covariance function parameters.
Component Symbol Role
Full covariance \(C(d \mid \psi)\) \(\sigma^2 \cdot \rho(d \mid \upsilon)\)
Correlation \(\rho(d \mid \phi, \nu)\) Equation (1) — Bessel function form
Variance \(\sigma^2\) Scales the covariance
Range \(\phi\) Controls decay speed (larger = faster)
Smoothness \(\nu\) Controls differentiability (larger = smoother)

Key insight: Equation (1) defines each entry of the covariance matrix based solely on the distance d between locations-a consequence of stationarity and isotropy.