Interpretacja testu zanurzeniowego Hartigansa

Chciałbym znaleźć sposób na oszacowanie intensywności bimodalności niektórych rozkładów, które uzyskałem empirycznie. Z tego, co przeczytałem, wciąż trwa debata na temat sposobu kwantyfikacji bimodalności. Zdecydowałem się na test zanurzeniowy Hartigansa, który wydaje się być jedynym dostępnym na R (oryginalny artykuł: http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). Test zanurzeniowy Hartigansa definiuje się jako: „Test zanurzeniowy mierzy multimodalność w próbce przez maksymalną różnicę we wszystkich punktach próbki między funkcją rozkładu empirycznego a funkcją rozkładu nieimodalnego, która minimalizuje tę maksymalną różnicę” .

Chciałbym całkowicie zrozumieć, jak powinienem interpretować te statystyki przed ich użyciem. Spodziewałem się, że test zanurzenia wzrośnie, jeśli rozkład jest multimodalny (ponieważ jest on definiowany jako „maksymalna różnica od rozkładu unimodalnego”). Ale : można przeczytać na stronie wikipedii o dystrybucji multimodalnej, że „Wartości mniejsze niż 0,05 wskazują na znaczącą bimodalność, a wartości większe niż 0,05, ale mniejsze niż 0,10 sugerują bimodalność o marginalnym znaczeniu”. . Takie stwierdzenie pochodzi z tego artykułu (ryc. 2). Zgodnie z tym artykułem wskaźnik testu zanurzenia jest bliski 0, gdy rozkład jest bimodalny. To mnie dezorientuje.

Aby poprawnie zinterpretować test zanurzeniowy Hartigansa, skonstruowałem niektóre rozkłady (oryginalny kod stąd ) i zwiększyłem wartość exp (mu2) (zwaną odtąd „Intensywnością bimodularności” - edycja : powinienem nazwać to „Intensywnością” bimodalności ” ), aby uzyskać bimodalność. Na pierwszym wykresie można zobaczyć przykładowy rozkład. Następnie oszacowałem wskaźnik diptestu (drugi wykres) i wartość p (trzeci wykres) powiązaną ( diptest pakietu ) z tymi różnymi symulowanymi rozkładami. Użyty kod R znajduje się na końcu mojego postu.

Pokazuję tutaj, że wskaźnik testu zanurzenia jest wysoki, a wartość Pvalue jest niska, gdy rozkłady są bimodalne. Co jest sprzeczne z tym, co można przeczytać w Internecie.

Nie jestem ekspertem w dziedzinie statystyki, więc ledwo zrozumiałem artykuł Hartigansa. Chciałbym uzyskać komentarze na temat właściwego sposobu interpretacji testu zanurzeniowego Hartigansa. Czy gdzieś się mylę?

Dziękuję wam wszystkim. Pozdrowienia,

TA

Przykład symulowanej dystrybucji:

Indeks testu zanurzeniowego Hartigana związany:

Test zanurzeniowy Hartigana p. Wartość związana:

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

Pan Freeman (autor artykułu , o którym ci mówiłem) powiedział mi, że tak naprawdę patrzy tylko na wartość testu zanurzenia. To zamieszanie wynika z jego zdania:
„Wartości HDS wynoszą od 0 do 1, przy wartościach mniejszych niż 0,05 wskazujących na znaczącą bimodalność, a wartości większych niż .05, ale mniejszych niż .10 sugerujących bimodalność o marginalnym znaczeniu” . Wartości HDS odpowiadają Pvalue, a nie statystykom testu zanurzeniowego. W gazecie nie było jasne.

Moja analiza jest dobra: statystyki testu zanurzenia rosną, gdy rozkład odbiega od rozkładu unimodalnego.

Test bimodalności i test Silvermana można również łatwo obliczyć w R i dobrze wykonać zadanie.