# Since Splus has difficulty with looping, we use two functions and sapply.
# In the function quan.call, j is just a dummy variable used for sapply.

get.quantiles_function(alpha, h = 3, m=500000)
{
M <- matrix(1, 1, 1)
k <- 2^(h - 1)
for(j in 2:h)
M <- rbind(rep(1, 2^(j - 1)), cbind(M, - M))
A <- eigen(1/k * t(M) %*% M, symmetric = T)$vectors[1:k,
1:h]
y <- sapply(1:m, quan.call, info = list(alpha = alpha, k
= k, A = A))
quantile(y, c(0.9, 0.95))
}

quan.call_function(info, j)
{
x <- NA
# If alpha is small, occasionally rstab will return NA. We throw out those samples.
while(any(is.na(x))) x <- rstab(info$k, info$alpha)
info$k^(1 - 2/info$alpha) * sum((t(x) %*% info$A)^2)
}