rbwheel {robustX}

Multivariate Barrow Wheel Distribution Random Vectors

Description

Generate p-dimensional random vectors according to Stahel's Barrow Wheel Distribution.

Usage

rbwheel(n, p, frac = 1/p, sig1 = 0.05, sig2 = 1/10,
        rGood = rnorm,
        rOut = function(n) sqrt(rchisq(n, p - 1)) * sign(runif(n, -1, 1)),
        U1 = rep(1, p),
        scaleAfter = TRUE, scaleBefore = FALSE, spherize = FALSE,
        fullResult = FALSE)

Arguments

n

integer, specifying the sample size.

p

integer, specifying the dimension (aka number of variables).

frac

numeric, the proportion of outliers. The default, 1/p, corresponds to the (asymptotic) breakdown point of M-estimators.

sig1

thickness of the “wheel”, (= \sigma (good[,1])), a non-negative numeric.

sig2

thickness of the “axis” (compared to 1).

rGood

function; the generator for “good” observations.

rOut

function, generating the outlier observations.

U1

p-vector to which (1,0,\dots,0) is rotated.

scaleAfter

logical indicating if the matrix is re-scaled after rotation (via scale()).. Default TRUE; note that this used to be false by default in the first public version.

scaleBefore

logical indicating if the matrix is re-scaled before rotation (via scale()).

spherize

logical indicating if the matrix is to be “spherized”, i.e., rotated and scaled to have empirical covariance I_p. This means that the principal components are used (before rotation).

fullResult

logical indicating if in addition to the n \times p matrix, some intermediate quantities are returned as well.

Details

....

Value

By default (when fullResult is FALSE), an n \times p matrix of n sample vectors of the p dimensional barrow wheel distribution, with an attribute, n1 specifying the exact number of “good” observations, n1 \approx (1-f)\cdot n, f = frac.

If fullResult is TRUE, a list with components

X

the n \times p matrix of above, X = X0 %*% A, where A <- Qrot(p, u = U1), and X0 is the corresponding matrix before rotation, see below.

X0

.........

A

the p \times p rotation matrix, see above.

n1

the number of “good” observations, see above.

n2

the number of “outlying” observations, n2 = n - n1.

Author(s)

Werner Stahel and Martin Maechler

References

http://stat.ethz.ch/people/maechler/robustness

Stahel, W.~A. and Mächler, M. (2009). Comment on “invariant co-ordinate selection”, Journal of the Royal Statistical Society B 71, 584–586. doi:10.1111/j.1467-9868.2009.00706.x

Examples

set.seed(17)
rX8 <- rbwheel(1000,8, fullResult = TRUE, scaleAfter=FALSE)
with(rX8, stopifnot(all.equal(X, X0 %*% A,    tol = 1e-15),
                    all.equal(X0, X %*% t(A), tol = 1e-15)))
##--> here, don't need to keep X0 (nor A, since that is Qrot(p))

## for n = 100,  you  don't see "it", but may guess .. :
n <- 100
pairs(r <- rbwheel(n,6))
n1 <- attr(r,"n1") ; pairs(r, col=1+((1:n) > n1))

## for n = 500, you *do* see it :
n <- 500
pairs(r <- rbwheel(n,6))
## show explicitly
n1 <- attr(r,"n1") ; pairs(r, col=1+((1:n) > n1))

## but increasing sig2 does help:
pairs(r <- rbwheel(n,6, sig2 = .2))

## show explicitly
n1 <- attr(r,"n1") ; pairs(r, col=1+((1:n) > n1))

set.seed(12)
pairs(X <- rbwheel(n, 7, spherize=TRUE))
colSums(X) # already centered

if(require("ICS") && require("robustbase")) {
  # ICS: Compare M-estimate [Max.Lik. of t_{df = 2}] with high-breakdown :
  stopifnot(require("MASS"))
  X.paM <- ics(X, S1 = cov, S2 = function(.) cov.trob(., nu=2)$cov, stdKurt = FALSE)
  X.paM.<- ics(X, S1 = cov, S2 = function(.) tM(., df=2)$V, stdKurt = FALSE)
  X.paR <- ics(X, S1 = cov, S2 = function(.) covMcd(.)$cov, stdKurt = FALSE)
  plot(X.paM) # not at all clear
  plot(X.paM.)# ditto
  plot(X.paR)# very clear
}
## Similar such experiments --->  demo(rbwheel_d)  and   demo(rbwheel_ics)
##                                --------------         -----------------