In [1]:

    

boottest2.med <- function(x, y, B) {
    n1 <- length(x)
    n2 <- length(y)
    
    v <- median(y) - median(x)
    z <- c(x,y)
    
    counter <- 0
    teststatall <- rep(0,B)
    
    for(i in 1:B) {
        xstar <- sample(z, n1, replace=TRUE)
        ystar <- sample(z, n2, replace=TRUE)
        vstar <- median(ystar) - median(xstar)
        if (vstar >= v) { counter <- counter + 1 }
        teststatall[i] <- vstar
    }
    pvalue <- counter/B
    list(origtest=v, pvalue=pvalue, teststatall=teststatall)
}


    

In [2]:

    

rcn <- function(n, eps, sigmac) {
  ind <- rbinom(n, 1, eps)
  x <- rnorm(n)
  rcn <- x*(1-ind)+sigmac*x*ind
  rcn
}


    

In [3]:

    

w <- rcn(30, 0.20, 4)
x <- 10 * w[1:15] + 100
y <- 10 * w[16:30] + 120


    

In [4]:

    

boottest2.med(x, y, 3000)$pvalue


    

    
Out[4]:




0.013

In [ ]: