# rbwheel: Multivariate Barrow Wheel Distribution Random Vectors In robustX: 'eXtra' / 'eXperimental' Functionality for Robust Statistics

## Description

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

## Usage

 ```1 2 3 4 5 6``` ```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”, (= σ `(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,…,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 x p matrix, some intermediate quantities are returned as well.

....

## Value

By default (when `fullResult` is `FALSE`), an n x 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 ~= (1-f)*n, f = `frac`.

If `fullResult` is `TRUE`, a list with components

 `X` the n x 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 x 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

Stahel, W.~A. and M<c3><a4>chler, M. (2009). Comment on “invariant co-ordinate selection”, Journal of the Royal Statistical Society B 71, 584–586.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39``` ```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")) { # 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) ## -------------- ----------------- ```

robustX documentation built on May 29, 2017, 9:37 p.m.