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
