rbwheel {robustX} | R Documentation |
Multivariate Barrow Wheel Distribution Random Vectors
Description
Generate -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, |
sig1 |
thickness of the “wheel”, ( |
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 |
scaleAfter |
logical indicating if the matrix is re-scaled after
rotation (via |
scaleBefore |
logical indicating if the matrix is re-scaled before
rotation (via |
spherize |
logical indicating if the matrix is to be
“spherized”, i.e., rotated and scaled to have empirical
covariance |
fullResult |
logical indicating if in addition to the |
Details
....
Value
By default (when fullResult
is FALSE
), an
matrix of
sample vectors of the
dimensional barrow wheel distribution, with an attribute,
n1
specifying the exact number of “good” observations,
,
frac
.
If fullResult
is TRUE
, a list with components
X |
the |
X0 |
......... |
A |
the |
n1 |
the number of “good” observations, see above. |
n2 |
the number of “outlying” observations, |
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)
## -------------- -----------------