# R/rbwheel-fun.R In robustX: 'eXtra' / 'eXperimental' Functionality for Robust Statistics

#### Documented in Qrotrbwheel

###
##
##  Qrot(), originally in  /u/maechler/R/MM/STATISTICS/robust/MV-WSt-ex.R
##  ----
##
## Idea: (alternative:  that rotates  e_i {i-th unit vector} on to
##	e_jk where (e_jk)_i = 1_[i \in {j,k}]

### Construct the Rotation matrix, that rotates (1,0,....0) onto ### (1,1,1,...1)/sqrt(p), or more generally to  u / ||u||  (u := unit.image)
Qrot <- function(p, transpose = FALSE, unit.image = rep(1,p))
{
## Purpose: Construct a  p x p  rotation matrix which rotates
##	(1,0,....0) on (1,1,,...1)/sqrt(p), i.e. the x1-axis onto the diagonal
## ----------------------------------------------------------------------
## Arguments: p : dimension
##          unit.image: vector onto which (1, 0,...,0) is rotated
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 3 Dec 2005; 28 Nov 2008
D <- diag(p)
## for some reason, the " - " below also keeps the sign correct:
Q <- qr.qy(qr(cbind(unit.image, D[,-1])), - D)
if(transpose) Q else t(Q)
}

rbwheel <- function(n,		# observations
p,		# variables
frac = 1/p, # proportion of outliers
sig1 = .05, # thickness of the 'wheel' = sigma(good[,1])
sig2 = 1/10,# diameter of the 'axle' (compared to 1)
rGood = rnorm,# generator for "good" observations
rOut = function(n)sqrt(rchisq(n,p-1))*sign(runif(n,-1,1)),
U1 = rep(1, p), ## Vector to which (1,0,...,0) is rotated
scaleAfter = TRUE, scaleBefore = FALSE, spherize = FALSE,
fullResult = FALSE)
{
## Purpose: simulate data according to the 'wheel' distribution
## ----------------------------------------------------------------------
## Arguments:
## ----------------------------------------------------------------------
## Author: Werner Stahel, Martin Maechler, Date: 28 Nov 2008, 15:27
stopifnot(is.numeric(frac), 0 <= frac, frac < 1,
n >= 1, p >= 2)
## compatibility-warning -- at most once per session :
## if(missing(scaleAfter) &&
##    (is.null(w <- getOption("rbwheel.warn.scaleA")) || isTRUE(w))) {
##     if(is.null(w)) options( rbwheel.warn.scaleA = FALSE)
##     warning("Note: rbwheel() now uses scaleAfter = TRUE  by default")
## }

## a simplified version of scale.default :
scale.simply <- function(x) {
x <- sweep(x, 2, colMeans(x, na.rm = TRUE), check.margin = FALSE)
sdev <- apply(x, 2,
function(v) {
v <- v[!is.na(v)]
sqrt(sum(v^2)/max(1, length(v) - 1L))
})
sweep(x, 2, sdev, "/", check.margin = FALSE)
}

n1 <- pmax(0, pmin(n, round((1-frac)*n)))
n2 <- n-n1 ## ~= frac * n

d0 <- matrix(0, n,p)
d0[ , -1] <- rGood(n*(p-1))
d1 <- sig1 * rGood(n1)
if(n1 < n) { ## have *some* outliers
i <- (n1+1):n
d0[i, -1] <- sig2 * d0[i,-1]
d0[-i, 1] <- d1
d0[ i, 1] <- rOut(n2)
} else { ## n1 == n; n2 == 0: no outliers
d0[, 1] <- d1
}

d1 <- {
if(spherize) { # use Chol(), such that X_1 remains unchanged:
d. <- scale.simply(d0)
t(backsolve(chol(cov(d.)), t(d.)))
}
else if(scaleBefore)
scale.simply(d0)
else d0
}

maybeScale <- function(x) if(scaleAfter) scale.simply(x) else x
if(fullResult) { ## for didactical reasons mainly, see example
A <- Qrot(p, unit.image = U1)
list(X = maybeScale(d1 %*% A), X0 = d0, A = A, n1 = n1, n2 = n2,
scale = c(before=scaleBefore, after=scaleAfter, spherize=spherize))
}
else ## by default -- 'n1' as attribute :
structure(maybeScale(d1 %*% Qrot(p, unit.image = U1)), n1 = n1)
}


## Try the robustX package in your browser

Any scripts or data that you put into this service are public.

robustX documentation built on May 2, 2019, 4:06 a.m.