sk_var {snapKrig}

Generate a covariance matrix or its factorization

Description

Computes the covariance matrix V (or one of its factorizations) for the non-NA points in sk grid g, given the model parameters list pars

Usage

sk_var(
  g,
  pars = NULL,
  scaled = FALSE,
  fac_method = "none",
  X = NULL,
  fac = NULL,
  sep = FALSE
)

Arguments

g

a sk grid object or a list with entries 'gdim', 'gres', 'gval'

pars

list of form returned by sk_pars (with entries 'y', 'x', 'eps', 'psill')

scaled

logical, whether to scale by 1/pars$psill

fac_method

character, the factorization to return, one of 'none', 'chol', 'eigen'

X

numeric matrix, the X in t(X) %*% V %*% X (default is identity, see details)

fac

matrix or list of matrices, the variance factorization (only used with X)

sep

logical, indicating to return correlation components instead of full covariance matrix

Details

By default the output matrix is V. Alternatively, if X is supplied, the function returns the quadratic form X^T V^-1 X.

When fac_method=='eigen' the function instead returns the eigen-decomposition of the output matrix, and when fac_method=='chol' its lower triangular Cholesky factor is returned. Supplying this factorization in argument fac in a subsequent call with X can speed up calculations. fac is ignored when X is not supplied.

scaled=TRUE returns the matrix scaled by the reciprocal of the partial sill, 1/pars$psill, before factorization. This is the form expected by functions sk_var_mult and sk_LL in argument fac.

Numerical precision issues with poorly conditioned covariance matrices can often be resolved by using 'eigen' factorization method (instead 'chol') and making sure that pars$eps > 0.

If all grid points are observed, then the output V becomes separable. Setting sep=TRUE in this case causes the function to return the x and y component correlation matrices (or their factorizations, as requested in fac_method) separately, in a list. scaled has no effect in this output mode. Note also that sep has no effect when X is supplied.

If the function is passed an empty grid g (all points NA) it returns results for the complete case (no NAs). If it is passed a list that is not a sk grid object, it must include entries 'gdim', 'gres', 'gval' and/or 'idx_grid' (as they are specified in sk), all other entries are ignored in this case.

Value

either matrix V, or X^T V^-1 X, or a factorization ('chol' or 'eigen')

See Also

sk

Other variance-related functions: sk_GLS(), sk_LL(), sk_cmean(), sk_nLL(), sk_sim()

Examples

# define example grid with NAs and example predictors matrix
gdim = c(12, 13)
n = prod(gdim)
n_obs = floor(n/3)
idx_obs = sort(sample.int(n, n_obs))
g = g_empty = sk(gdim)
g[idx_obs] = stats::rnorm(n_obs)
plot(g)

# example kernel
psill = 0.3
pars = utils::modifyList(sk_pars(g), list(psill=psill))

# plot the covariance matrix for observed data, its cholesky factor and eigen-decomposition
V_obs = sk_var(g, pars)
V_obs_chol = sk_var(g, pars, fac_method='chol')
V_obs_eigen = sk_var(g, pars, fac_method='eigen')
sk_plot(V_obs)
sk_plot(V_obs_chol)
sk_plot(V_obs_eigen$vectors)

# empty and complete cases are treated the same

# get the full covariance matrix with sep=FALSE (default)
V_full = sk_var(g_empty, pars)

# check that the correct sub-matrix is there
max(abs( V_obs - V_full[idx_obs, idx_obs] ))

# get 1d correlation matrices with sep=TRUE...
corr_components = sk_var(g_empty, pars, sep=TRUE)
str(corr_components)
sk_plot(corr_components[['x']])

# ... these are related to the full covariance matrix through psill and eps
corr_mat = kronecker(corr_components[['x']], corr_components[['y']])
V_full_compare = pars$psill * corr_mat + diag(pars$eps, n)
max(abs(V_full - V_full_compare))

# ... their factorizations can be returned as (nested) lists
str(sk_var(g_empty, pars, fac_method='chol', sep=TRUE))
str(sk_var(g_empty, pars, fac_method='eigen', sep=TRUE))

# compare to the full covariance matrix factorizations (default sep=FALSE)
str(sk_var(g_empty, pars, fac_method='chol'))
str(sk_var(g_empty, pars, fac_method='eigen'))

# test quadratic form with X
nX = 3
X_all = cbind(1, matrix(stats::rnorm(nX * n), ncol=nX))
cprod_all = crossprod(X_all, chol2inv(chol(V_full))) %*% X_all
abs(max(sk_var(g_empty, pars, X=X_all) - cprod_all ))

# test products with inverse of quadratic form with X
mult_test = stats::rnorm(nX + 1)
cprod_all_inv = chol2inv(chol(cprod_all))
cprod_all_inv_chol = sk_var(g_empty, pars, X=X_all, scaled=TRUE, fac_method='eigen')
sk_var_mult(mult_test, pars, fac=cprod_all_inv_chol) - cprod_all_inv %*% mult_test

# repeat with missing data
X_obs = X_all[idx_obs,]
cprod_obs = crossprod(X_obs, chol2inv(chol(V_obs))) %*% X_obs

abs(max(sk_var(g, pars, X=X_obs) - cprod_obs ))
cprod_obs_inv = chol2inv(chol(cprod_obs))
cprod_obs_inv_chol = sk_var(g, pars, X=X_obs, scaled=TRUE, fac_method='eigen')
sk_var_mult(mult_test, pars, fac=cprod_obs_inv_chol) - cprod_obs_inv %*% mult_test

# `scaled` indicates to divide matrix by psill
print( pars[['eps']]/pars[['psill']] )
diag(sk_var(g, pars, scaled=TRUE)) # diagonal elements equal to 1 + eps/psill
( sk_var(g, pars) - psill * sk_var(g, pars, scaled=TRUE) ) |> abs() |> max()
( sk_var(g, pars, X=X_obs, scaled=TRUE) - ( cprod_obs/psill ) ) |> abs() |> max()

# in Cholesky factor this produces a scaling by square root of psill
max(abs( V_obs_chol - sqrt(psill) * sk_var(g, pars, fac_method='chol', scaled=TRUE) ))

# and in the eigendecomposition, a scaling of the eigenvalues
vals_scaled = sk_var(g, pars, fac_method='eigen', scaled=TRUE)$values
max(abs( sk_var(g, pars, fac_method='eigen')$values - psill*vals_scaled ))