# data from Table 2.8.1 in the text

before <- c(16.55, 15.36, 15.94, 16.43, 16.01) 


after <- c(16.05, 15.98, 16.10, 15.88, 15.91)


 # a function to calculate rmdhat values 

rmdcalc   <-   function(x,y)
      {
        x <- x[!is.na(x)]
        y <- y[!is.na(y)]
        rmdcalc  <- mean(abs(x -median(x)))/mean(abs(y -median(y)))
        rmdcalc
      }

rmdcalc(before,after)




 # for the one-sided, upper tail test for deviances with unknown location as shown in the text


rmdperm  <-  function(x,y,iter)
  {
     x <- x[!is.na(x)]
     y <- y[!is.na(y)]
     rmdobs <- rmdcalc(x,y) ; print(c("The observed RMDhat value = ",rmdobs))
     rmd  <- numeric(iter)

     for (i in 1:iter)
       {
          devx     <- x - median(x)
          devy     <- y - median(y)
          both     <- c(devx,devy)
          bothsamp <- sample(both)
          devxnew  <- bothsamp[1:length(x)]
          devynew  <- bothsamp[(length(x)+1):(length(x)+length(y))]
          rmd[i]   <- mean(abs(devxnew))/mean(abs(devynew))
         # print(xnew) ; print(ynew) ; print(rmd[i])
       }
     # print(rmd)

     qrmdperm <- quantile(rmd,c(.01,.025,.05,.1,.25,.5,.75,.9,.95,.975,.99))
     print("Quantiles for RMDhat") ; print(qrmdperm)
     permp  <- sum((rmd >= rmdobs)/iter)     
     permp
  }

  # these values should be close to those on page 55

rmdperm(before,after,1000)