chol_s2Dp {sim2Dpredictr}

Build and Take the Cholesky Decomposition of a Covariance Matrix

Description

The function first builds a correlation matrix with correlation.builder, converts that matrix to a covariance matrix if necessary, and then takes the Cholesky decomposition of the matrix using either base R or the R package spam. Note that spam is particularly effective when the matrix is sparse.

Usage

chol_s2Dp(
  matrix.type = "cov",
  im.res,
  use.spam = FALSE,
  corr.structure = "ar1",
  rho = NULL,
  phi = NULL,
  tau = 1,
  alpha = 0.75,
  corr.min = NULL,
  neighborhood = "none",
  w = NULL,
  h = NULL,
  r = NULL,
  print.R = FALSE,
  print.S = FALSE,
  print.Q = FALSE,
  sigma = 1,
  triangle = "upper",
  print.all = FALSE,
  round.d = FALSE,
  return.cov = TRUE,
  return.prec = TRUE
)

Arguments

matrix.type

Determines whether to build a covariance matrix, "cov", or a precision matrix, "prec". See correlation_builder{sim2Dpredictr} and precision_builder{sim2Dpredictr} for more details.

im.res

A vector defining the dimension of spatial data. The first entry is the number of rows and the second entry is the number of columns.

use.spam

If use.spam = TRUE then use tools from the R package spam; otherwise, base R functions are employed. For large dimension MVN with sparse correlation structure, spam is recommended; otherwise, base R may be faster. Defaults to FALSE.

corr.structure

One of "ar1", exponential, gaussian, or "CS". Correlations between locations i and j are rho^{d} for corr.structure = "ar1", exp(-phi * d) for corr.structure = "exponential", exp(-phi * d ^ 2) for corr.structure = "gaussian", and rho when corr.structure = "CS". Note that d is the Euclidean distance between locations i and j.

rho

This is the maximum possible correlation between locations i and j. For all i,j rho MUST be between -1 and 1.

phi

A scalar value greater than 0 that determines the decay rate of correlation. This argument is only utilized when corr.structure %in% c("exponential", "gaussian").

tau

A vector containing precision parameters. If of length 1, then all precisions are assumed equal. Otherwise the length of tau should equal the number of variables.

alpha

A scalar value between 0 and 1 that defines the strength of correlations. Note that when alpha = 0 the data are independent and when alpha = 1, the joint distribution is the improper Intrinsic Autoregression (IAR), which cannot be used to generate data. Note also that while alpha does control dependence it is not interpretable as a correlation.

corr.min

Scalar value to specify the minimum non-zero correlation. Any correlations below corr.min are set to 0. Especially for high image resolution using this option can result in a sparser covariance matrix, which may significantly speed up draws when using spam. This option is preferred to using neighborhood and associated arguments when the primary concern is to avoid very small correlations and improve computation efficiency. Default is NULL, which places no restrictions on the correlations.

neighborhood

Defines the neighborhood within which marginal correlations are non-zero. The default is "none", which allows marginal correlations to extend indefinitely. neighborhood = "round" defines a circular neighborhood about locations and neighborhood = "rectangle" defines a rectangular neighborhood about locations. Note that this argument differs from that in precision_builder, in which neighborhood defines conditional non-zero correlations.

w, h

If neighborhood = "rectangle" then w and h are the number of locations to the left/right and above/below a location i that define its neighborhood. Any locations outside this neighborhood have have zero correlation with location i.

r

If neighborhood = "round", then if locations i,j are separated by distance d \ge r, the correlation between them is zero.

print.R, print.S, print.Q

Logical. When TRUE, then print the correlation, covariance, or precision matrix before taking the Cholesky decomposition. If sigma = 1, then S = R.

sigma

Specify the desired standard deviations; the default is 1, in which case the Cholesky decomposition is of a correlation matrix. If sigma != 1, then the Cholesky decomposition is of a covariance Matrix.

• If sigma is a vector then length(sigma) must be equal to the total number of locations, i.e. (n.row * n.col) by (n.row * n.col).

• sigma can take any scalar value when specifying common standard deviation.

triangle

Determine whether to output an upper (triangle = "upper") or lower (triangle = "lower") triangular matrix.

print.all

If print.all = TRUE, then prints each correlation and allows you to check whether the correlations are as you intended. This option is NOT recommended for large point lattices/images.

round.d

If round.d = TRUE, then d is rounded to the nearest whole number.

return.cov, return.prec

Logical. When TRUE, also return the covariance or precision matrix, respectively. This is recommended when using spam to generate draws from the MVN.

Value

Matrix of dimension (n.row x n.col) x (n.row x n.col). If either return.cov or return.prec is TRUE, then returns a list where the first element is the covariance or precision matrix, and the second element is the Cholesky factor.

References

Banerjee S, Carlin BP, Gelfand AE (2015). Hierarchical Modeling and Analysis for Spatial Data, Second edition. Chapman & Hall/CRC, Boca Raton, Florida.

Ripley BD (1987). Stochastic Simulation. John Wiley & Sons. doi:10.1002/9780470316726.

Rue H (2001). “Fast Sampling of Gaussian Markov Random Fields.” Journal of the Royal Statistical Society B, 63, 325-338. doi:10.1111/1467-9868.00288.

Furrer R, Sain SR (2010). “spam: A Sparse Matrix R Package with Emphasis on MCMC Methods for Gaussian Markov Random Fields.” Journal of Statistical Software, 36(10), 1-25. https://www.jstatsoft.org/v36/i10/.

Examples

## Use R package spam for Cholesky decomposition
R_spam <- chol_s2Dp(im.res = c(3, 3), matrix.type = "prec",
                    use.spam = TRUE, neighborhood = "ar1",
                    triangle = "upper")

## Use base R for Cholesky decomposition
R_base <- chol_s2Dp(corr.structure = "ar1", 
                    im.res = c(3, 3), rho = 0.15,
                    neighborhood = "round", 
                    r = 3, use.spam = FALSE)

## Specify standard deviations instead of default of sigma = 1.
R_sd <- chol_s2Dp(corr.structure = "ar1",
                  im.res = c(3, 3), rho = 0.15,
                  neighborhood = "round", r = 3, 
                  sigma = runif(9, 1.1, 4))
## Not run: 
## Print options ON
R_pr_on <- chol_s2Dp(corr.structure = "ar1", 
                     im.res = c(3, 3), rho = 0.15,
                     sigma = 1:9, neighborhood = "round", 
                     r = 3, print.R = TRUE, print.S = TRUE)
 
## End(Not run)