tests/sim1.R

In diptest: Hartigan's Dip Test Statistic for Unimodality - Corrected

#### Very small scale simulation to make the point
#### --> See ../stuff/  for much more
library(diptest)

P.p <- c(1, 5, 10, 25)/100
(P.p <- c(P.p, 1/2, rev(1 - P.p)))

N.sim <- 9999
set.seed(94)
.p0 <- proc.time()
dU100 <- replicate(N.sim, dip(runif(100)))
cat('Time elapsed: ', (p1 <- proc.time()) - .p0,'\n'); .p0 <- p1
## Lynne (2003: P IV, 1.6 GHz): ~7   s
## 2010 (AMD Phenom II X4 925):  1.3 s

100 * round(q100 <- quantile(dU100, p = P.p), 4)

plot(density(sqrt(100) * dU100), lwd = 2, col=2,
     main = expression("Dip  distribution" ~~
         list(sqrt(n)* D[n], ~ n == 100)))
abline(h=0, col="dark gray", lty=3)

round(1e4 * quantile(dU100, p = seq(0,1, by = 0.01), names = FALSE))

##--- an extreme unimodal case -- i.e. very small dip():
set.seed(60); x <- rexp(301,1)^3
hist(x)
(dt.x <- dip.test(x))
(dt2 <- dip.test(x, simulate = TRUE))
(dt3 <- dip.test(x, simulate = TRUE, B = 10000))
stopifnot(dt.x$p.value == 1,## <- gave	NA  earlier
	  dt2$p.value == 1,
	  dt3$p.value == 1)


cat('Time elapsed: ', proc.time() - .p0,'\n') # "stats"