1 |
x |
|
y |
|
alpha |
|
q |
|
plotit |
|
pop |
|
fr |
|
rval |
|
xlab |
|
ylab |
|
nboot |
|
SEED |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (x, y, alpha = 0.05, q = 0.25, plotit = FALSE, pop = 0,
fr = 0.8, rval = 15, xlab = "", ylab = "", nboot = 600, SEED = TRUE)
{
library(parallel)
if (SEED)
set.seed(2)
if (q >= 0.5)
stop("q should be less than .5")
if (q <= 0)
stop("q should be greater than 0")
x <- x[!is.na(x)]
y <- y[!is.na(y)]
n1 = length(x)
n2 = length(y)
m <- outer(x, y, FUN = "-")
q2 = 1 - q
est1 = hd(m, q)
est2 = hd(m, q2)
data1 <- matrix(sample(n1, size = n1 * nboot, replace = TRUE),
nrow = nboot)
data2 <- matrix(sample(n2, size = n2 * nboot, replace = TRUE),
nrow = nboot)
data = cbind(data1, data2)
data = listm(t(data))
bvec = NA
bvec <- mclapply(data, cbmhd_subMC, x = x, y = y, q = q,
q2 = q2, n1 = n1, n2 = n2, mc.preschedule = TRUE)
bvec = list2vec(bvec)
p = mean(bvec > 0) + 0.5 * mean(bvec == 0)
p = 2 * min(c(p, 1 - p))
sbv = sort(bvec)
ilow <- round((alpha/2) * nboot)
ihi <- nboot - ilow
ilow <- ilow + 1
ci = sbv[ilow]
ci[2] = sbv[ihi]
if (plotit) {
if (pop == 1 || pop == 0) {
if (length(x) * length(y) > 2500) {
print("Product of sample sizes exceeds 2500.")
print("Execution time might be high when using pop=0 or 1")
print("If this is case, might consider changing the argument pop")
print("pop=2 might be better")
}
}
MM = as.vector(m)
if (pop == 0)
akerd(MM, xlab = xlab, ylab = ylab)
if (pop == 1)
rdplot(MM, fr = fr, xlab = xlab, ylab = ylab)
if (pop == 2)
kdplot(MM, rval = rval, xlab = xlab, ylab = ylab)
if (pop == 3)
boxplot(MM)
if (pop == 4)
stem(MM)
if (pop == 5)
hist(MM, xlab = xlab)
if (pop == 6)
skerd(MM)
}
list(q = q, Est1 = est1, Est2 = est2, sum = est1 + est2,
ci = ci, p.value = p)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.