robsbss: Robust Spatial Blind Source Separation

In SpatialBSS: Blind Source Separation for Multivariate Spatial Data

Robust Spatial Blind Source Separation

Description

robsbss is a robust variant of sbss. It estimates the unmixing matrix assuming a spatial blind source separation model by jointly diagonalizing the Hettmansperger-Randles scatter matrix and one/many generalized local sign covariance matrices. These local generalized sign covariance matrices are determined by spatial kernel functions and radial functions. Three types of such kernel functions and three types of radial functions are supported.

Usage

robsbss(x, ...)

## Default S3 method:
robsbss(x, coords, kernel_type = c('ring', 'ball', 'gauss'), 
     kernel_parameters, lcov = c('norm', 'winsor', 'qwinsor'), 
     ordered = TRUE, kernel_list = NULL, ...)
## S3 method for class 'SpatialPointsDataFrame'
robsbss(x, ...)
## S3 method for class 'sf'
robsbss(x, ...)

Arguments

x	either a numeric matrix of dimension c(n, p) where the p columns correspond to the entries of the random field and the n rows are the observations, an object of class SpatialPointsDataFrame or an object of class sf.
coords	a numeric matrix of dimension c(n,2) where each row represents the coordinates of a point in the spatial domain. Only needed if x is a matrix and the argument kernel_list is NULL.
kernel_type	a string indicating which kernel function to use. Either 'ring' (default), 'ball' or 'gauss'.
kernel_parameters	a numeric vector that gives the parameters for the kernel function. At least length of one for 'ball' and 'gauss' or two for 'ring' kernel, see details.
lcov	a string indicating which radial function or type of robust local covariance matrix to use. Either 'norm' (default), 'winsor' or 'qwinsor'. See also local_gss_covariance_matrix for details.
ordered	logical. If TRUE the entries of the latent field are ordered by the sum of squared (pseudo-)eigenvalues of the diagonalized local covariance matrix/matrices. Default is TRUE.
kernel_list	a list of spatial kernel matrices with dimension c(n,n), see details. Usually computed by the function spatial_kernel_matrix.
...	further arguments for the fast real joint diagonalization algorithm that jointly diagonalizes the local covariance matrices. See details and frjd.

Details

robsbss is a robust variant of sbss which uses Hettmansperger-Randles (HR) location and scatter estimates (Hettmansperger & Randles, 2002) for whitening (see white_data for details) and jointly diagonalizes HR scatter matrix and generalized local sign matrices to estimate the unmixing matrix. The generalized local sign matrices are determined by radial functions w(l_i), where l_i = ||x(s_i)-T(x)|| and T(x) is HR location estimator, and kernel functions f(d_{i,j}), where d_{i,j}=||s_i - s_j||. Generalized local sign covariance (gLSCM) matrix is then calculated as

gLSCM(f,w) = 1/(n F^{1/2}_{f,n}) \sum_{i,j} f(d_{i,j}) w(l_i)w(l_j)(x(s_i)-T(x)) (x(s_j)-T(x))'

with

F_{f,n} = 1 / n \sum_{i,j} f^2(d_{i,j}).

Three radial functions (Raymaekers & Rousseeuw, 2019) w(l_i) are implemented, the parameter lcov defines which is used:

• 'norm':

  w(l_i) = 1/l_i

• 'winsor':

  w(l_i) = Q/l_i

• 'qwinsor':

  w(l_i) = Q^2/l_i^2.

The cutoff Q is defined as Q = l_{(h)}, where l_{(h)} is hth order statistic of \{l_1, ..., l_n\} and h = (n + p + 1)/2. In addition, three kernel functions f(d) are implemented, the parameter kernel_type defines which is used:

• 'ring': parameters are inner radius r_{in} and outer radius r_{out}, with r_{in} < r_{out}, and r_{in}, r_{out} \ge 0:

  f(d;r_{in}, r_{out}) = I(r_{in} < d \le r_{out})

• 'ball': parameter is the radius r, with r \ge 0:

  f(d;r) = I(d \le r)

• 'gauss': Gaussian function where 95% of the mass is inside the parameter r, with r \ge 0:

  f(d;r) = exp(-0.5 (\Phi^{-1}(0.95) d/r)^2).

The argument kernel_type determines the used kernel function as presented above, the argument kernel_parameters gives the corresponding parameters for the kernel function. Specifically, if kernel_type equals 'ball' or 'gauss' then kernel_parameters is a numeric vector where each entry corresponds to one parameter. Hence, length(kernel_parameters) local covariance matrices are used. Whereas, if kernel_type equals 'ring', then kernel_parameters must be a numeric vector of even length where subsequently the inner and outer radii must be given (informally: c(r_in1, r_out1, r_in2, r_out2, ...)). In that case length(kernel_parameters) / 2 local covariance matrices are used.

robsbss calls spatial_kernel_matrix internally to compute a list of c(n,n) kernel matrices based on the parameters given, where each entry of those matrices corresponds to f(d_{i,j}). Alternatively, such a list of kernel matrices can be given directly to the function robsbss via the kernel_list argument. This is useful when robsbss is called numerous times with the same coordinates/kernel functions as the computation of the kernel matrices is then done only once prior the actual robsbss calls. For details see also spatial_kernel_matrix.

If more than one generalized local sign covariance matrix is used robsbss jointly diagonalizes these matrices with the function frjd. ... provides arguments for frjd, useful arguments might be:

• eps: tolerance for convergence.

• maxiter: maximum number of iterations.

Value

robsbss returns a list of class 'sbss' with the following entries:

s	object of class(x) containing the estimated source random field.
coords	coordinates of the observations. Is NULL if x was a matrix and the argument kernel_list was not NULL at the robsbss call.
w	estimated unmixing matrix.
weights	numeric vector of length(n) giving the weights for each observation for the robust local covariance estimation.
w_inv	inverse of the estimated unmixing matrix.
pevals	(pseudo-)eigenvalues for each latent field entry.
d	matrix of stacked (jointly) diagonalized local covariance matrices with dimension c(length(kernel_parameters)*p,p) for 'ball' and 'gauss' kernel or c( (length(kernel_parameters) / 2)*p,p) for 'ring' kernel.
diags	matrix of dimension c(length(kernel_parameters),p) where the rows contain the diagonal of the diagonalized local autocovariance matrices.
x_mu	robustly estimated columnmeans of x.
cov_inv_sqrt	square root of the inverse sample covariance matrix of x.

References

Hettmansperger, T. P., & Randles, R. H. (2002). A practical affine equivariant multivariate median. Biometrika, 89 , 851-860. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/biomet/89.4.851")}.

Raymaekers, J., & Rousseeuw, P. (2019). A generalized spatial sign covariance matrix. Journal of Multivariate Analysis, 171 , 94-111. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.jmva.2018.11.010")}.

Sipila, M., Muehlmann, C. Nordhausen, K. & Taskinen, S. (2022). Robust second order stationary spatial blind source separation using generalized sign matrices. Manuscript.

See Also

spatial_kernel_matrix, local_gss_covariance_matrix, sp, sf, frjd

Examples

# simulate coordinates
coords <- runif(1000 * 2) * 20
dim(coords) <- c(1000, 2)
coords_df <- as.data.frame(coords)
names(coords_df) <- c("x", "y")
# simulate random field
if (!requireNamespace('gstat', quietly = TRUE)) {
  message('Please install the package gstat to run the example code.')
} else {
  library(gstat)
  model_1 <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE, beta = 0, 
                   model = vgm(psill = 0.025, range = 1, model = 'Exp'), nmax = 20)
  model_2 <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE, beta = 0, 
                   model = vgm(psill = 0.025, range = 1, kappa = 2, model = 'Mat'), 
                   nmax = 20)
  model_3 <- gstat(formula = z ~ 1, locations = ~ x + y, dummy = TRUE, beta = 0, 
                   model = vgm(psill = 0.025, range = 1, model = 'Gau'), nmax = 20)
  field_1 <- predict(model_1, newdata = coords_df, nsim = 1)$sim1
  field_2 <- predict(model_2, newdata = coords_df, nsim = 1)$sim1
  field_3 <- predict(model_3, newdata = coords_df, nsim = 1)$sim1
  field <- cbind(field_1, field_2, field_3)
  # Generate 5 % local outliers to data
  field_cont <- gen_loc_outl(field, coords, radius = 2, 
                             swap_order = "regular")[,1:3]
  X <- as.matrix(field_cont)
  
  # apply sbss with three ring kernels
  kernel_parameters <- c(0, 1, 1, 2, 2, 3)
  robsbss_result <- 
    robsbss(X, coords, kernel_type = 'ring', kernel_parameters = kernel_parameters)
  
  # print object
  print(robsbss_result)
  
  # plot latent field
  plot(robsbss_result, colorkey = TRUE, as.table = TRUE, cex = 1)
  
  # predict latent fields on grid
  predict(robsbss_result, colorkey = TRUE, as.table = TRUE, cex = 1)
  
  # unmixing matrix
  w_unmix <- coef(robsbss_result)
}

SpatialBSS documentation built on July 26, 2023, 5:37 p.m.