I am doing a project on the Gap Statistic from Tibshirani etc http://www.web.stanford.edu/%7Ehastie/Papers/gap.pdf
On page 4 of the pdf in section 4 on "The computational implementation of the gap statistic" two possible reference distributions to select Monte-Carlo samples from:
a) generate each reference feature uniformly overthe range of observed values for that feature
b) generate reference features from a uniform distribution over a box aligned with the principal components of the data. (The details involve assuming the columns have mean zero and computing the singular value decomposition, transforming the data, picking uniform random samples over the ranges of the columns in the transformed matrix, then back transforming the sample.)
The comments say the latter takes account of the shape of the distribution, and makes the procedure rotationally invariant, as long as the clustering method is invariant.
I can't tell from the help which method is used in the R function clusGap in the cluster library. I have tried looking at the code for the function (see below) and I think it doesn't do the singular value decomposition b) method - but just does the a) method. (I was looking at the part of the code following for (b in 1:B)). My question is am I right?
> clusGap
function (x, FUNcluster, K.max, B = 100, verbose = interactive(),
...)
{
stopifnot(is.function(FUNcluster), length(dim(x)) == 2, K.max >=
2, (n <- nrow(x)) >= 1, ncol(x) >= 1)
if (B != (B. <- as.integer(B)) || (B <- B.) <= 0)
stop("'B' has to be a positive integer")
if (is.data.frame(x))
x <- as.matrix(x)
ii <- seq_len(n)
W.k <- function(X, kk) {
clus <- if (kk > 1)
FUNcluster(X, kk, ...)$cluster
else rep.int(1L, nrow(X))
0.5 * sum(vapply(split(ii, clus), function(I) {
xs <- X[I, , drop = FALSE]
sum(dist(xs)/nrow(xs))
}, 0))
}
logW <- E.logW <- SE.sim <- numeric(K.max)
if (verbose)
cat("Clustering k = 1,2,..., K.max (= ", K.max, "): .. ",
sep = "")
for (k in 1:K.max) logW[k] <- log(W.k(x, k))
if (verbose)
cat("done\n")
xs <- scale(x, center = TRUE, scale = FALSE)
m.x <- rep(attr(xs, "scaled:center"), each = n)
V.sx <- svd(xs)$v
rng.x1 <- apply(xs %*% V.sx, 2, range)
logWks <- matrix(0, B, K.max)
if (verbose)
cat("Bootstrapping, b = 1,2,..., B (= ", B, ") [one \".\" per sample]:\n",
sep = "")
for (b in 1:B) {
z1 <- apply(rng.x1, 2, function(M, nn) runif(nn, min = M[1],
max = M[2]), nn = n)
z <- tcrossprod(z1, V.sx) + m.x
for (k in 1:K.max) {
logWks[b, k] <- log(W.k(z, k))
}
if (verbose)
cat(".", if (b%%50 == 0)
paste(b, "\n"))
}
if (verbose && (B%%50 != 0))
cat("", B, "\n")
E.logW <- colMeans(logWks)
SE.sim <- sqrt((1 + 1/B) * apply(logWks, 2, var))
structure(class = "clusGap", list(Tab = cbind(logW, E.logW,
gap = E.logW - logW, SE.sim), n = n, B = B, FUNcluster = FUNcluster))
}
<bytecode: 0x0000000010248578>
<environment: namespace:cluster>
Warning messages:
1: In .HTMLsearch(query) : Unrecognized search field: title
2: In .HTMLsearch(query) : Unrecognized search field: keyword
3: In .HTMLsearch(query) : Unrecognized search field: alias
