R: Robust Spatial Blind Source Separation

robsbss {SpatialBSS}

R Documentation

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. 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. 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.

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)
}

[Package SpatialBSS version 0.14-0 Index]