rbwheel function

Multivariate Barrow Wheel Distribution Random Vectors

Multivariate Barrow Wheel Distribution Random Vectors

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

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/p1/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,,0)(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 IpI_p. This means that the principal components are used (before rotation).

  • fullResult: logical indicating if in addition to the nxpn x p matrix, some intermediate quantities are returned as well.

Details

....

Returns

By default (when fullResult is FALSE), an nxpn x p matrix of nn sample vectors of the pp dimensional barrow wheel distribution, with an attribute, n1 specifying the exact number of good observations, n1 =(1f)nn1 ~= (1-f)*n, f=f = frac.

If fullResult is TRUE, a list with components - X: the nxpn 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 pxpp x p rotation matrix, see above.

  • n1: the number of good observations, see above.

  • n2: the number of outlying observations, n2=nn1n2 = n - n1.

References

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

Stahel, W.~A. and , M. (2009). Comment on ``invariant co-ordinate selection'', Journal of the Royal Statistical Society B 71 , 584--586. tools:::Rd_expr_doi("10.1111/j.1467-9868.2009.00706.x")

Author(s)

Werner Stahel and Martin Maechler

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) ## -------------- -----------------
  • Maintainer: Martin Maechler
  • License: GPL (>= 2)
  • Last published: 2023-06-16

Useful links