# data from Table 2.6.2 granite <- c(33.63, 39.86, 69.32, 42.13, 58.36, 74.11) basalt <- c(26.15, 18.56, 17.55, 9.84, 28.29, 34.15) # this function conducts the permutation test for the one-sided, upper tail test # based on van der Waerden scores # Note that the function name is vWperm, not vwperm (W is capitalized) vWperm <- function(x,y,iter) { x <- x[!is.na(x)] y <- y[!is.na(y)] vwobs <- sum(qnorm(rank(c(x,y))/(length(c(x,y))+1))[1:length(x)]) print(c("The observed vW value = ",vwobs)) vw <- numeric(iter) for (i in 1:iter) { both <- c(x,y) bothsamp <- sample(both) vwxnew <- sum(qnorm(rank(bothsamp)/(length(bothsamp)+1))[1:length(x)]) vw[i] <- vwxnew # print(bothsamp) ; print(vw[i]) } # print(vw) qvwperm <- quantile(vw,c(.01,.025,.05,.1,.25,.5,.75,.9,.95,.975,.99)) print("Quantiles for the vW statistic ") ; print(qvwperm) permp <- sum((vw >= vwobs)/iter) permp } # these values should be close to those on page 50 vWperm(granite,basalt,1000) rm(granite, basalt) # remove objects when finished