R Under development (unstable) (2023-11-26 r85638 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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.77 0.04 0.81 NA NA > ## 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.67 0.22 1.91 NA NA > > proc.time() user system elapsed 2.65 0.32 2.98