R help - Jul 2003 - Silverman modality test

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)