Como testar se minha distribuição é multimodal?

Quando plogo um histograma dos meus dados, ele tem dois picos:

Isso significa uma potencial distribuição multimodal? Eu executei o dip.testem R ( library(diptest)) e a saída é:

D = 0.0275, p-value = 0.7913

Posso concluir que meus dados têm uma distribuição multimodal?

DADOS

10346 13698 13894 19854 28066 26620 27066 16658  9221 13578 11483 10390 11126 13487 
15851 16116 24102 30892 25081 14067 10433 15591  8639 10345 10639 15796 14507 21289 
25444 26149 23612 19671 12447 13535 10667 11255  8442 11546 15958 21058 28088 23827 
30707 19653 12791 13463 11465 12326 12277 12769 18341 19140 24590 28277 22694 15489 
11070 11002 11579  9834  9364 15128 15147 18499 25134 32116 24475 21952 10272 15404 
13079 10633 10761 13714 16073 23335 29822 26800 31489 19780 12238 15318  9646 11786 
10906 13056 17599 22524 25057 28809 27880 19912 12319 18240 11934 10290 11304 16092 
15911 24671 31081 27716 25388 22665 10603 14409 10736  9651 12533 17546 16863 23598 
25867 31774 24216 20448 12548 15129 11687 11581

O @NickCox apresentou uma estratégia interessante (+1). No entanto, posso considerá-lo de natureza mais exploratória, devido à preocupação que @whuber aponta .

Deixe-me sugerir outra estratégia: você pode ajustar um modelo de mistura finita gaussiana. Observe que isso pressupõe fortemente que seus dados são extraídos de uma ou mais normais verdadeiras. Como @whuber e @NickCox apontam nos comentários, sem uma interpretação substantiva desses dados - apoiada por uma teoria bem estabelecida - para apoiar essa suposição, essa estratégia também deve ser considerada exploratória.

Primeiro, vamos seguir a sugestão de @ Glen_b e examinar seus dados usando o dobro de posições:

Ainda vemos dois modos; se alguma coisa, eles aparecem mais claramente aqui. (Observe também que a linha de densidade do kernel deve ser idêntica, mas parece mais espalhada devido ao maior número de posições.)

Agora vamos ajustar um modelo de mistura finita gaussiana. Em R, você pode usar o Mclustpacote para fazer isso:

library(mclust)
x.gmm = Mclust(x)
summary(x.gmm)
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm 
# ----------------------------------------------------
#   
# Mclust V (univariate, unequal variance) model with 2 components:
#   
#   log.likelihood   n df       BIC       ICL
#        -1200.874 120  5 -2425.686 -2442.719
# 
# Clustering table:
#  1  2 
# 68 52 

Dois componentes normais otimizam o BIC. Para comparação, podemos forçar um ajuste de um componente e executar um teste de razão de verossimilhança:

x.gmm.1 = Mclust(x, G=1)
logLik(x.gmm.1)
# 'log Lik.' -1226.241 (df=2)
logLik(x.gmm)-logLik(x.gmm.1)
# 'log Lik.' 25.36657 (df=5)
1-pchisq(25.36657, df=3)  # [1] 1.294187e-05

Isso sugere que é extremamente improvável que você encontre dados tão unodais quanto os seus se eles vierem de uma única distribuição normal verdadeira.

Algumas pessoas não se sentem confortáveis ​​usando um teste paramétrico aqui (embora, se as premissas se mantiverem, não conheço nenhum problema). Uma técnica amplamente aplicável é usar o Método de ajuste cruzado de inicialização paramétrica (descrevo o algoritmo aqui ). Podemos tentar aplicá-lo a estes dados:

x.gmm$parameters
# $mean
# 12346.98 23322.06 
# $variance$sigmasq
# [1]  4514863 24582180
x.gmm.1$parameters
# $mean
# [1] 17520.91
# $variance$sigmasq
# [1] 43989870

set.seed(7809)
B = 10000;    x2.d = vector(length=B);    x1.d = vector(length=B)
for(i in 1:B){
  x2      = c(rnorm(68, mean=12346.98, sd=sqrt( 4514863)), 
              rnorm(52, mean=23322.06, sd=sqrt(24582180)) )
  x1      = rnorm( 120, mean=17520.91, sd=sqrt(43989870))
  x2.d[i] = Mclust(x2, G=2)$loglik - Mclust(x2, G=1)$loglik
  x1.d[i] = Mclust(x1, G=2)$loglik - Mclust(x1, G=1)$loglik
}
x2.d = sort(x2.d);  x1.d = sort(x1.d)
summary(x1.d)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -0.29070 -0.02124  0.41460  0.88760  1.36700 14.01000 
summary(x2.d)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  9.006  23.770  27.500  27.760  31.350  53.500 

As estatísticas de resumo e os gráficos de densidade do kernel para as distribuições de amostragem mostram vários recursos interessantes. A probabilidade de log para o modelo de componente único raramente é maior que a do ajuste de dois componentes, mesmo quando o verdadeiro processo de geração de dados tem apenas um componente e, quando é maior, a quantidade é trivial. A idéia de comparar modelos que diferem em sua capacidade de ajustar dados é uma das motivações por trás do PBCM. As duas distribuições de amostragem quase não se sobrepõem; apenas 0,35% x2.dsão inferiores ao valor máximox1.dvalor. Se você selecionasse um modelo de dois componentes se a diferença na probabilidade de log fosse> 9,7, você selecionaria incorretamente o modelo de um componente 0,01% e o modelo de dois componentes 0,02% do tempo. Estes são altamente discrimináveis. Se, por outro lado, você optar por usar o modelo de um componente como hipótese nula, o resultado observado será suficientemente pequeno para não aparecer na distribuição empírica da amostra em 10.000 iterações. Podemos usar a regra 3 (veja aqui ) para colocar um limite superior no valor-p, ou seja, estimamos que seu valor-p seja menor que 0,0003. Ou seja, isso é altamente significativo.

$p < .000001$$p < .001$) e os componentes subjacentes, caso existam, também não são garantidos como perfeitamente normais. Se você acha razoável que seus dados possam vir de uma distribuição positiva inclinada, e não normal, esse nível de bimodalidade pode estar dentro do intervalo típico de variação, que é o que suspeito que o teste de mergulho esteja dizendo.

No seguimento das ideias em @ resposta e os comentários de Nick, você pode ver o quão grande as necessidades de largura de banda para ser a apenas alisar o modo secundário:

Pegue essa estimativa de densidade do kernel como nulo proximal - a distribuição mais próxima dos dados, mas ainda consistente com a hipótese nula de que é uma amostra de uma população unimodal - e simule a partir dela. Nas amostras simuladas, o modo secundário nem sempre parece tão distinto, e você não precisa ampliar tanto a largura de banda para nivelá-la.

A formalização dessa abordagem leva ao teste realizado em Silverman (1981), "Usando estimativas de densidade de núcleo para investigar a modalidade", JRSS B , 43 , 1. O silvermantestpacote de Schwaiger & Holzmann implementa esse teste e também o procedimento de calibração descrito por Hall & York ( 2001), "Na calibração do teste de Silverman para multimodalidade", Statistica Sinica , 11 , p 515, que ajusta o conservadorismo assintótico. A realização do teste em seus dados com uma hipótese nula de unimodalidade resulta em valores de p de 0,08 sem calibração e 0,02 com calibração. Não estou familiarizado o suficiente com o teste de mergulho para adivinhar por que ele pode ser diferente.

Código R:

  # kernel density estimate for x using Sheather-Jones method to estimate b/w:
density(x, kernel="gaussian", bw="SJ") -> dens.SJ
  # tweak b/w until mode just disappears:
density(x, kernel="gaussian", bw=3160) -> prox.null
  # fill matrix with simulated samples from the proximal null:
x.sim <- matrix(NA, nrow=length(x), ncol=10)
for (i in 1:10){
  x.sim[ ,i] <- rnorm(length(x), sample(x, size=length(x), replace=T), prox.null$bw)
}
  # perform Silverman test without Hall-York calibration:
require(silvermantest)
silverman.test(x, k=1, M=10000, adjust=F)
  # perform Silverman test with Hall-York calibration:
silverman.test(x, k=1, M=10000, adjust=T)

As coisas para se preocupar incluem:

1. O tamanho do conjunto de dados. Não é pequeno, nem grande.

2. A dependência do que você vê na origem do histograma e na largura da lixeira. Com apenas uma opção evidente, você (e nós) não tem idéia de sensibilidade.

3. A dependência do que você vê no tipo e largura do kernel e quaisquer outras opções feitas na estimativa de densidade. Com apenas uma opção evidente, você (e nós) não tem idéia de sensibilidade.

Em outros lugares, sugeri provisoriamente que a credibilidade dos modos é suportada (mas não estabelecida) por uma interpretação substantiva e pela capacidade de discernir a mesma modalidade em outros conjuntos de dados do mesmo tamanho. (Maior é melhor também ....)

Não podemos comentar sobre nenhum deles aqui. Uma pequena possibilidade de repetibilidade é comparar o que você obtém com amostras de bootstrap do mesmo tamanho. Aqui estão os resultados de um experimento de token usando o Stata, mas o que você vê é arbitrariamente limitado aos padrões do Stata, que são documentados como arrancados do ar . Eu tenho estimativas de densidade para os dados originais e para 24 amostras de inicialização do mesmo.

A indicação (nem mais, nem menos) é o que acho que analistas experientes adivinhariam de qualquer maneira do seu gráfico. O modo esquerdo é altamente repetível e o direito é claramente mais frágil.

Observe que há uma inevitabilidade nisso: como há menos dados próximos do modo direito, eles nem sempre reaparecem em uma amostra de autoinicialização. Mas este também é o ponto chave.

Observe que o ponto 3. acima permanece intocado. Mas os resultados estão em algum lugar entre unimodal e bimodal.

Para os interessados, este é o código:

clear 
set scheme s1color 
set seed 2803 

mat data = (10346, 13698, 13894, 19854, 28066, 26620, 27066, 16658, 9221, 13578, 11483, 10390, 11126, 13487, 15851, 16116, 24102, 30892, 25081, 14067, 10433, 15591, 8639, 10345, 10639, 15796, 14507, 21289, 25444, 26149, 23612, 19671, 12447, 13535, 10667, 11255, 8442, 11546, 15958, 21058, 28088, 23827, 30707, 19653, 12791, 13463, 11465, 12326, 12277, 12769, 18341, 19140, 24590, 28277, 22694, 15489, 11070, 11002, 11579, 9834, 9364, 15128, 15147, 18499, 25134, 32116, 24475, 21952, 10272, 15404, 13079, 10633, 10761, 13714, 16073, 23335, 29822, 26800, 31489, 19780, 12238, 15318, 9646, 11786, 10906, 13056, 17599, 22524, 25057, 28809, 27880, 19912, 12319, 18240, 11934, 10290, 11304, 16092, 15911, 24671, 31081, 27716, 25388, 22665, 10603, 14409, 10736, 9651, 12533, 17546, 16863, 23598, 25867, 31774, 24216, 20448, 12548, 15129, 11687, 11581)
set obs `=colsof(data)' 
gen data = data[1,_n] 

gen index = . 

quietly forval j = 1/24 { 
    replace index = ceil(120 * runiform()) 
    gen data`j' = data[index]
    kdensity data`j' , nograph at(data) gen(xx`j' d`j') 
} 

kdensity data, nograph at(data) gen(xx d) 

local xstuff xtitle(data/1000) xla(10000 "10" 20000 "20" 30000 "30") sort 
local ystuff ysc(r(0 .0001)) yla(none) `ystuff'   

local i = 1 
local colour "orange" 
foreach v of var d d? d?? { 
    line `v' data, lc(`colour') `xstuff'  `ystuff' name(g`i', replace) 
    local colour "gs8" 
    local G `G' g`i' 
    local ++i 
} 

graph combine `G' 

Identificação não paramétrica do modo LP (nome do algoritmo LPMode , a ref do artigo é fornecida abaixo)

MaxEntModos [triângulos de cor vermelha na plotagem]: 12783.36 e 24654.28.

L2Modos [triângulos de cor verde na plotagem]: 13054.70 e 24111.61.

Interessante notar as formas modais, especialmente a segunda que mostra considerável distorção (o modelo da Mistura Gaussiana Tradicional provavelmente falhará aqui).

Mukhopadhyay, S. (2016) Identificação do modo em larga escala e ciências orientadas a dados. https://arxiv.org/abs/1509.06428