robustX (version 1.2-4)

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

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

Details

....

References

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

Stahel, W.~A. and M<U+00E4>chler, M. (2009). Comment on ``invariant co-ordinate selection'', Journal of the Royal Statistical Society B 71, 584--586.

Examples

# NOT RUN {
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)
##                                --------------         -----------------
# }