Dear R users, I've written the following functions to implement Silverman's modality test ("Using kernel density estimates to investigate multimodality", J.R. Stat. Soc. B, 43, (1981), 97-99.), and I tested them on the Chondrite data set (Good & Gaskins, J. Amer. Stat. Ass., 75, (1980), 42-56). Values for the critical window width seem OK, which is not the case for the significance levels.
If someone could give me a hint about what is wrong... Or perhaps someone has already done a real implementation of this test? Thanks, Jerome Sackur Inserm U562 Orsay, France nbmodes <- function (x, h) # returns how many modes there are in the kernel density estimate of # vector x, with window width h. { modes <- density (x, bw=h) modes <- diff (diff (modes$y) / abs (diff (modes$y))) modes <- rep(1, length(modes))[modes==-2] modes <- sum (modes) return (modes) } hcrit <- function (x, n, e=.0001) # Returns the minimal window width such that the kernel density estimate # of x has n modes. { minw <- min (abs (diff (x))) maxw <- (max (x) - min (x))/2 winN <- maxw winN1 <- minw while (abs(winN-winN1)>e) { modes <- nbmodes (x, winN) winN1 <- winN if (modes > n) { minw <- winN winN <- (winN + maxw)/2 } else { maxw <- winN winN <- (winN - minw)/2 } } return (winN) } silversignif <- function (x, m, nboot=1000) # Silverman's significance for the null that x has at # most modes. { h <- hcrit (x, m) h0 <- h n <- 0 theta <- function (y, h0=h, sigma2=var(y)) { (y + rnorm (1, mean=0, sd=h0^2)) / sqrt ( 1+ (h0^2/sigma2)) } for (i in 1:nboot) { boot <- theta(sample(x, replace=T)) nb <- nbmodes (boot, h0) if (nb > m) { n <- n+1 } } return (n/nboot) } # Chondrite meteor data chond <- c(20.77, 22.56, 22.71, 22.99, 26.39, 27.08, 27.32, 27.33, 27.57, 27.81, 28.69, 29.36, 30.25, 31.89, 32.88, 33.23, 33.28, 33.40, 33.52, 33.83, 33.95, 34.82) ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help