dip: Compute Hartigans' Dip Test Statistic for Unimodality

Description

Computes Hartigans' dip test statistic for testing unimodality, and additionally the modal interval.

Usage

dip(x, full.result = FALSE, min.is.0 = FALSE, debug = FALSE)

Arguments

numeric; the data.

full.result

logical or string; dip(., full.result=TRUE) returns the full result list; if "all" it additionally uses the mn and mj components to compute the initial GCM and LCM, see below.

min.is.0

logical indicating if the minimal value of the dip statistic $D_n$ can be zero or not. Arguably should be set to TRUE for internal consistency reasons, but is false by default both for continuity and backwards compatibilit

debug

logical; if true, some tracing information is printed (from the C routine).

Value

depending on full.result either a number, the dip statistic, or an object of class "dip" which is a list with components
xthe sorted unname()d data.
nlength(x).
dipthe dip statistic
lo.hiindices into x for lower and higher end of modal interval
xl, xulower and upper end of modal interval
gcm, lcm(last used) indices for greatest convex minorant and the least concave majorant.
mn, mjindex vectors of length n for the GC minorant and the LC majorant respectively.
For full results of class "dip", there are print and plot methods, the latter with its own manual page.

References

P. M. Hartigan (1985) Computation of the Dip Statistic to Test for Unimodality; Applied Statistics (JRSS C) 34, 320--325. Corresponding (buggy!) Fortran code of AS 217 available from Statlib, http://lib.stat.cmu.edu/apstat/217

J. A. Hartigan and P. M. Hartigan (1985) The Dip Test of Unimodality; Annals of Statistics 13, 70--84.

Examples

Run this code

data(statfaculty)
plot(density(statfaculty))
rug(statfaculty, col="midnight blue"); abline(h=0, col="gray")
dip(statfaculty)
(dS <- dip(statfaculty, full = TRUE, debug = TRUE))
plot(dS)
## even more output -- + plot showing "global" GCM/LCM:
(dS2 <- dip(statfaculty, full = "all", debug = 3))
plot(dS2)

data(faithful)
fE <- faithful$eruptions
plot(density(fE))
rug(fE, col="midnight blue"); abline(h=0, col="gray")
dip(fE, debug = 2) ## showing internal work
(dE <- dip(fE, full = TRUE)) ## note the print method
plot(dE, do.points=FALSE)

data(precip)
plot(density(precip))
rug(precip, col="midnight blue"); abline(h=0, col="gray")
str(dip(precip, full = TRUE, debug = TRUE))

##-----------------  The  'min.is.0' option :  ---------------------

##' dip(.) continuity and 'min.is.0' exploration:
dd <- function(x, debug=FALSE) {
   x_ <- x ; x_[1] <- 0.9999999999 * x[1]
   rbind(dip(x , debug=debug),
         dip(x_, debug=debug),
         dip(x , min.is.0=TRUE, debug=debug),
         dip(x_, min.is.0=TRUE, debug=debug), deparse.level=2)
}

dd( rep(1, 8) ) # the 3rd one differs ==> min.is.0=TRUE is *dis*continuous
dd( 1:7 )       # ditto

dd( 1:7, debug=TRUE)
## border-line case ..
dd( 1:2, debug=TRUE)

## Demonstrate that  'min.is.0 = TRUE'  does not change the typical result:
B.sim <- 1000 # or larger
D5  <- {set.seed(1); replicate(B.sim, dip(runif(5)))}
D5. <- {set.seed(1); replicate(B.sim, dip(runif(5), min.is.0=TRUE))}
stopifnot(identical(D5, D5.), all.equal(min(D5), 1/(2*5)))
hist(D5, 64); rug(D5)

D8  <- {set.seed(7); replicate(B.sim, dip(runif(8)))}
D8. <- {set.seed(7); replicate(B.sim, dip(runif(8), min.is.0=TRUE))}
stopifnot(identical(D8, D8.))

Run the code above in your browser using DataCamp Workspace