
R Under development (unstable) (2015-03-03 r67931) -- "Unsuffered Consequences"
Copyright (C) 2015 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #### 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)))
[1] 0.01 0.05 0.10 0.25 0.50 0.75 0.90 0.95 0.99
> 
> 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
Time elapsed:  0.512 0.085 0.597 0 0 
> ## 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)
  1%   5%  10%  25%  50%  75%  90%  95%  99% 
2.29 2.56 2.75 3.08 3.54 4.12 4.70 5.09 5.90 
> 
> 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))
  [1] 191 229 239 246 252 256 261 265 268 272 275 277 280 282 285 287 289 292
 [19] 294 296 298 300 302 304 305 308 310 312 314 315 317 319 321 323 325 327
 [37] 329 331 332 334 336 338 340 341 343 345 347 349 351 352 354 356 358 360
 [55] 362 364 366 368 370 372 374 376 379 381 383 385 387 390 393 395 397 400
 [73] 403 406 409 412 415 418 421 424 427 431 434 438 442 446 450 455 460 464
 [91] 470 476 483 489 499 509 520 539 562 590 773
> 
> ##--- an extreme unimodal case -- i.e. very small dip():
> set.seed(60); x <- rexp(301,1)^3
> hist(x)
> (dt.x <- dip.test(x))

	Hartigans' dip test for unimodality / multimodality

data:  x
D = 0.0072617, p-value = 1
alternative hypothesis: non-unimodal, i.e., at least bimodal

> (dt2 <- dip.test(x, simulate = TRUE))

	Hartigans' dip test for unimodality / multimodality with simulated
	p-value (based on 2000 replicates)

data:  x
D = 0.0072617, p-value = 1
alternative hypothesis: non-unimodal, i.e., at least bimodal

> (dt3 <- dip.test(x, simulate = TRUE, B = 10000))

	Hartigans' dip test for unimodality / multimodality with simulated
	p-value (based on 10000 replicates)

data:  x
D = 0.0072617, p-value = 1
alternative hypothesis: non-unimodal, i.e., at least bimodal

> 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"
Time elapsed:  1.136 0.018 1.158 0 0 
> 
> proc.time()
   user  system elapsed 
  1.765   0.129   1.935 
