Interpretación de la prueba de inmersión de Hartigans

Me gustaría encontrar una manera de cuantificar la intensidad de la bimodalidad de algunas distribuciones que obtuve empíricamente. Por lo que leí, todavía hay cierto debate sobre la forma de cuantificar la bimodalidad. Elegí usar la prueba de inmersión de Hartigans, que parece ser la única disponible en R (documento original: http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). La prueba de inmersión de Hartigans se define como: "La prueba de inmersión mide la multimodalidad en una muestra por la diferencia máxima, sobre todos los puntos de muestra, entre la función de distribución empírica y la función de distribución unimodal que minimiza esa diferencia máxima" .

Me gustaría entender completamente cómo debo interpretar estas estadísticas antes de usarlas. Esperaba que la prueba de inmersión aumentara si la distribución es multimodal (como se define como "la diferencia máxima de la distribución unimodal"). Pero : puede leer en la página de Wikipedia acerca de la distribución multimodal que "Los valores inferiores a 0.05 indican bimodalidad significativa y los valores superiores a 0.05 pero menores a 0.10 sugieren bimodalidad con significación marginal". . Dicha declaración proviene de este documento (Fig. 2). Según este documento, el índice de prueba de inmersión es cercano a 0 cuando la distribución es bimodal. Me confunde.

Para interpretar correctamente la prueba de inmersión de Hartigans, construí algunas distribuciones (el código original es de aquí ) y aumenté el valor de exp (mu2) (llamado 'Intensidad de bimodularidad' de ahora en adelante - Editar: Debería haberlo llamado 'Intensidad de bimodalidad ' ) para obtener bimodalidad. En el primer gráfico, puede ver algunos ejemplos de distribuciones. Luego estimé el índice diptest (segundo gráfico) y el valor p (tercer graphe) asociado (paquete de prueba ) a esas diferentes distribuciones simuladas. El código R utilizado está al final de mi publicación.

Lo que muestro aquí es que el índice de prueba de caída es alto y el Pvalue es bajo cuando las distribuciones son bimodales. Lo cual es contrario a lo que puedes leer en internet.

No soy un experto en estadística, así que apenas entendí el artículo de Hartigans. Me gustaría recibir algunos comentarios sobre la forma correcta en que debemos interpretar la prueba de inmersión de Hartigans. ¿Me equivoco en alguna parte?

Gracias a todos. Saludos,

ejército de reserva

Ejemplo de distribución simulada:

El índice de prueba de inmersión de Hartigan asociado:

Prueba de inmersión de Hartigan con valor de p asociado:

library(diptest)
library(ggplot2)


# CONSTANT PARAMETERS
sig1 <- log(3)
sig2 <- log(3)
cpct <- 0.5
N=1000

#CREATING BIMOD DISTRIBUTION
bimodalDistFunc <- function (n,cpct, mu1, mu2, sig1, sig2) {
  y0 <- rlnorm(n,mean=mu1, sd = sig1)
  y1 <- rlnorm(n,mean=mu2, sd = sig2)

  flag <- rbinom(n,size=1,prob=cpct)
  y <- y0*(1 - flag) + y1*flag 
}

#DIP TEST
DIP_TEST <- function(bimodalData) {
  TEST <- dip.test(bimodalData)
  return(TEST$statistic[[1]])   # return(TEST$p.value[[1]])    to get the p value
}
DIP_TEST(bimodalData)


# SIMULATION
exp_mu1 = 1
max_exp_mu2 = 100
intervStep = 100
repPerInt = 10

# single distibutions
expMu2Value <- c()
bimodalData <- c()
mu1 <- log(exp_mu1)   
mu2 <- log(exp_mu1)
bimodalData <- c(bimodalData,log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)))
expMu2Value <- c(expMu2Value,rep(exp_mu1,length(log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)))))

mu1 <- log(exp_mu1)   
mu2 <- log(max_exp_mu2)
bimodalData <- c(bimodalData,log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)))
expMu2Value <- c(expMu2Value,rep(max_exp_mu2,length(log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)))))

mu1 <- log(exp_mu1)   
mu2 <- log(trunc((max_exp_mu2-exp_mu1)/2+1))
bimodalData <- c(bimodalData,log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)))
expMu2Value <- c(expMu2Value,rep(trunc((max_exp_mu2-exp_mu1)/2+1),length(log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2)))))

tableExamples <- data.frame(expMu2Value,bimodalData)
tableExamples$expMu2Value <- as.factor(tableExamples$expMu2Value)
ExamplePlot <- ggplot(tableExamples)+
  geom_histogram(aes(bimodalData),color='white')+
  ylab("Count")+
  xlab("")+
  facet_wrap(~expMu2Value)+
  ggtitle("Intensity of bimodularity")

# calculation of the dip test index
exp_mu2Int = seq(from=exp_mu1,to=max_exp_mu2,length.out=intervStep)
expmu2Vec = c()
dipStat = c()
testDone = c()
for(exp_mu2 in exp_mu2Int){
  mu1 <- log(exp_mu1)   
  mu2 <- log(exp_mu2)
  for(rep in 1:repPerInt){
    bimodalData <- log(bimodalDistFunc(n=N,cpct,mu1,mu2, sig1,sig2))
    diptestone = DIP_TEST(bimodalData)
    expmu2Vec = c(expmu2Vec,exp_mu2)
    dipStat = c(dipStat,diptestone)
    testDone = c(testDone,"diptest")
  }
}
table = data.frame(expmu2Vec,dipStat,testDone)

IndexPlot <- ggplot(table)+
  geom_point(aes(expmu2Vec,dipStat,color=testDone))+
  ylab("Index")+
  xlab("Intensity of Bimodularity")+
  scale_color_discrete(name="Test")

ExamplePlot
IndexPlot

El Sr. Freeman (autor del artículo que le conté) me dijo que en realidad solo estaba mirando el Pvalue de la prueba de inmersión. Esta confusión proviene de su oración:
"Los valores de HDS varían de 0 a 1 con valores inferiores a .05 que indican una bimodalidad significativa, y valores mayores que .05 pero menores que .10 que sugieren una bimodalidad con significación marginal" . Los valores de HDS corresponden al Pvalue y no a las estadísticas de la prueba de inmersión. No estaba claro en el periódico.

Mi análisis es bueno: las estadísticas de la prueba de inmersión aumentan cuando la distribución se desvía de una distribución unimodal.

La prueba de bimodalidad y la prueba de Silverman también se pueden calcular fácilmente en R y hacer el trabajo bien.