Chapter 1: Robust Directional Foundations

The foundation functions in HDElliptical use observations in rows and variables in columns. They separate radial magnitude from angular information, which is the central computational advantage of spatial-sign methods under heavy tails.

Simulating an elliptical sample

The stochastic representation is

\[ X = \mu + \xi A U, \qquad AA^\top = \Sigma, \]

where \(U\) is uniform on the unit sphere and \(\xi\geq 0\) is independent of \(U\). The default radius in relliptical() gives a Gaussian sample; supplying a custom radial function changes the tail behavior without changing the shape.

library(HDElliptical)

set.seed(1)
shape <- matrix(c(2, 0.6, 0.6, 1), 2)
x <- relliptical(
  150,
  location = c(1, -1),
  shape = shape,
  radial = function(n) abs(rt(n, df = 3))
)

Location and directional covariance

The spatial median minimizes average Euclidean distance. Its diagnostics are stored as attributes so that the usual return value remains a numeric vector.

center <- spatial_median(x)
center
#> [1]  1.0015968 -0.9671267
#> attr(,"objective")
#> [1] 1.381363
#> attr(,"iterations")
#> [1] 29
#> attr(,"converged")
#> [1] TRUE
#> attr(,"relative_change")
#> [1] 0
#> attr(,"equation_residual")
#> [1] 8.072e-09
attr(center, "converged")
#> [1] TRUE
attr(center, "equation_residual")
#> [1] 8.072e-09

sign_shape <- sscm(x, center = center)
sum(diag(sign_shape))
#> [1] 1

The pairwise spatial Kendall matrix is translation invariant and avoids direct location estimation.

kendall_shape <- spatial_kendall(x)
rank_shape <- spatial_rank_covariance(x)
c(kendall_trace = sum(diag(kendall_shape)),
  rank_trace = sum(diag(rank_shape)))
#> kendall_trace    rank_trace 
#>     1.0000000     0.3882124

Affine-equivariant shape

Tyler’s estimator repeatedly reweights observations by their current Mahalanobis radii and normalizes the result to trace \(p\).

tyler <- tyler_shape(x)
sum(diag(tyler))
#> [1] 2
attr(tyler, "converged")
#> [1] TRUE
attr(tyler, "equation_residual")
#> [1] 4.223832e-09

When both location and affine-equivariant shape are needed, the Hettmansperger-Randles estimator solves the two estimating equations jointly.

hr <- hr_estimator(x)
hr$location
#> [1]  0.9986638 -0.9709963
hr$shape
#>           [,1]      [,2]
#> [1,] 1.1946641 0.3056217
#> [2,] 0.3056217 0.8053359
c(location = hr$location_equation_residual,
  shape = hr$shape_equation_residual)
#>     location        shape 
#> 6.366186e-09 7.873106e-09

Exact Tyler and HR estimators require more observations than variables and general position. The functions stop on zero centered residuals or singular updates rather than returning an invalid matrix.