Интерпретация провала Хартиганса

Я хотел бы найти способ количественно оценить интенсивность бимодальности некоторых распределений, которые я получил эмпирически. Из того, что я прочитал, до сих пор идут споры о том, как количественно определить бимодальность. Я решил использовать тест Хартиганса, который кажется единственным, доступным на R (оригинал статьи: http://www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf ). Тест на погружение Хартиганса определяется следующим образом: «Тест на погружение измеряет мультимодальность в выборке по максимальной разности по всем точкам выборки между эмпирической функцией распределения и функцией унимодального распределения, которая минимизирует эту максимальную разницу» .

Я хотел бы полностью понять, как я должен интерпретировать эту статистику перед ее использованием. Я ожидал, что тест на провал увеличится, если распределение будет мультимодальным (так как оно определено как «максимальное отличие от унимодального распределения»). Но : вы можете прочитать на странице википедии о мультимодальном распределении, что «значения менее 0,05 указывают на значительную бимодальность, а значения более 0,05, но менее 0,10 указывают на бимодальность с минимальной значимостью». , Такое утверждение исходит из этой статьи (рис. 2). Согласно этой статье, индекс испытания на провал близок к 0, когда распределение является бимодальным. Это смущает меня.

Чтобы правильно интерпретировать тест падения Хартиганса, я построил несколько распределений (исходный код отсюда ) и увеличил значение exp (mu2) (которое теперь называется «интенсивность бимодулярности») - Правка: я должен был назвать его «Интенсивность». бимодальности » ), чтобы получить бимодальность. На первом графике вы можете увидеть некоторые примеры распределений. Затем я оценил индекс diptest (второй график) и значение p (третий граф), связанный (пакет diptest ) с этими различными моделируемыми распределениями. Код R используется в конце моего поста.

Здесь я показываю, что индекс теста на провал высокий, а значение Pvalue низкое, когда распределения бимодальные. Что противоречит тому, что вы можете прочитать в Интернете.

Я не специалист по статистике, так что я едва понял статью Хартиганса. Я хотел бы получить некоторые комментарии о том, как правильно интерпретировать тест на падение Хартиганса. Я где то не прав?

Спасибо вам всем. С Уважением,

Т.А.

Пример распределения смоделирован:

Индекс теста на провал Хартигана связан с:

Падение Хартигана по п.значению связано:

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

Мистер Фримен (автор статьи, о которой я вам рассказывал) сказал мне, что он на самом деле смотрел только на значение теста на падение. Эта путаница исходит из его предложения:
«Значения HDS варьируются от 0 до 1 со значениями менее 0,05, указывающими значительную бимодальность, и значениями более 0,05, но менее 0,10, указывающими бимодальность с предельным значением» . Значения HDS соответствуют Pvalue, а не статистике теста на падение. Это было неясно в газете.

Мой анализ хорош: статистика теста на провал увеличивается, когда распределение отклоняется от одномодального распределения.

Тест на бимодальность и тест Сильвермана также можно легко вычислить в R и хорошо выполнять свою работу.