机器算法验证 - 解释 Hartigans 的浸渍测试 - 吾爱随笔录

我想找到一种方法来量化我凭经验获得的某些分布的双峰强度。从我读到的内容来看，关于量化双峰的方式仍然存在一些争论。我选择使用 Hartigans 的 dip 测试，这似乎是 R 上唯一可用的测试（原始论文： http: //www.stat.washington.edu/wxs/Stat593-s03/Literature/hartigan85a.pdf）。Hartigans 的倾角测试被定义为：“倾角测试通过经验分布函数与使最大差异最小化的单峰分布函数之间的所有样本点的最大差异来测量样本中的多峰性”。

我想在使用它之前完全理解我应该如何解释这个统计数据。如果分布是多峰的（因为它被定义为“与单峰分布的最大差异”），我期待着倾角测试会增加。但是：您可以在关于多峰分布的维基百科页面中阅读“小于 0.05 的值表示显着的双峰，大于 0.05 但小于 0.10 的值表示具有边际意义的双峰”。. 这种说法来自这篇论文（图2）。根据本文，当分布为双峰时，倾角测试指数接近于0。这让我很困惑。

为了正确解释 Hartigans 的倾角测试，我构建了一些分布（原始代码来自这里）并增加了 exp(mu2) 的值（从现在开始称为“双模强度” -编辑：我应该称之为“强度”的双峰'）获得双峰。在第一张图中，您可以看到一些分布示例。然后我估计了与这些不同的模拟分布相关联的 diptest 指数（第二张图）和 p 值（第三张图）（package diptest ）。使用的 R 代码在我的帖子末尾。

我在这里展示的是，当分布是双峰时，dip test index 很高，Pvalue 很低。这与您在互联网上阅读的内容相反。

我不是统计学专家，所以我几乎不理解 Hartigans 的论文。我想得到一些关于我们应该如何解释 Hartigans 的浸入测试的评论。我在某个地方错了吗？

谢谢你们。问候，

助教

模拟分布示例：

Hartigan 的倾角测试指数相关：在此处输入图像描述

Hartigan 的倾角测试 p.value 相关：在此处输入图像描述

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