Как проверить, является ли мой дистрибутив мультимодальным?

21

Когда я строю гистограмму моих данных, она имеет два пика:

гистограмма

Означает ли это потенциальное мультимодальное распределение? Я запустил dip.testв R ( library(diptest)), и вывод:

D = 0.0275, p-value = 0.7913

Я могу заключить, что мои данные имеют мультимодальное распределение?

ДАННЫЕ

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
user1260391
источник
3
Используйте больше бинов в своей гистограмме. Я предлагаю вдвое больше
Glen_b
1
В этом ответе упоминаются девять различных тестов , некоторые из которых могут иметь отношение к вашей ситуации.
Glen_b
1
Эта статья может быть полезна для вас, если вы не видели его уже (и это последующие )
Eoin

Ответы:

15

@NickCox представил интересную стратегию (+1). Я мог бы счесть это более исследовательским по своей природе, однако, из-за беспокойства, на которое указывает @whuber .

Позвольте мне предложить другую стратегию: вы можете использовать гауссовскую модель конечной смеси. Обратите внимание, что это делает очень сильное предположение, что ваши данные взяты из одной или нескольких истинных нормалей. Как указывают @whuber и @NickCox в комментариях, без существенной интерпретации этих данных, подкрепленной устоявшейся теорией, в поддержку этого предположения, эту стратегию также следует считать исследовательской.

Во-первых, давайте последуем предложению @ Glen_b и посмотрим на ваши данные, используя в два раза больше бинов:

введите описание изображения здесь

Мы все еще видим два режима; во всяком случае, они проходят здесь более четко. (Обратите также внимание, что линия плотности ядра должна быть одинаковой, но выглядит более разбросанной из-за большего количества бинов.)

Теперь давайте подгоним гауссову модель конечной смеси. В R, вы можете использовать Mclustпакет для этого:

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 

Два обычных компонента оптимизируют БИК. Для сравнения мы можем принудительно выбрать однокомпонентный компонент и выполнить тест отношения правдоподобия:

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

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

Некоторые люди не чувствуют себя комфортно, используя здесь параметрический тест (хотя, если предположения верны, я не знаю ни одной проблемы). Одним из широко применяемых методов является использование метода параметрической начальной загрузки (я опишу алгоритм здесь ). Мы можем попробовать применить его к этим данным:

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 

введите описание изображения здесь

Сводная статистика и графики плотности ядра для распределений выборки показывают несколько интересных особенностей. Вероятность записи в журнал для однокомпонентной модели редко превышает вероятность соответствия двух компонентов, даже если в процессе создания истинных данных имеется только один компонент, а когда он больше, это количество тривиально. Идея сравнения моделей, которые различаются по своей способности подгонять данные, является одной из причин, лежащих в основе PBCM. Два распределения выборки практически не перекрываются; только 0,35% x2.dменьше максимальногоx1.dзначение. Если вы выбрали двухкомпонентную модель, если разница в вероятности записи была> 9,7, вы бы неправильно выбрали однокомпонентную модель 0,01% и двухкомпонентную модель 0,02% времени. Это очень разборчиво. Если, с другой стороны, вы решили использовать однокомпонентную модель в качестве нулевой гипотезы, ваш наблюдаемый результат достаточно мал, чтобы не отображаться в эмпирическом распределении выборки за 10 000 итераций. Мы можем использовать правило 3 (см. Здесь ), чтобы установить верхнюю границу для значения p, а именно, мы оцениваем, что ваше значение p меньше, чем .0003. То есть это очень важно.

п<0,000001п<+0,001) и нижележащие компоненты, если они существуют, также не гарантируют, что они будут совершенно нормальными. Если вы считаете разумным, чтобы ваши данные могли поступать из положительно искаженного распределения, а не из нормального, этот уровень бимодальности вполне может находиться в пределах типичного диапазона отклонения, что, как я подозреваю, говорит тест на провал.

Gung - Восстановить Монику
источник
1
Проблема этого подхода в том, что альтернатива, с которой вы сравниваете гауссову смесь, не очень разумна. Более разумным было бы то, что дистрибутив является каким-то искаженным, например, гамма. Почти очевидно, что смесь будет соответствовать практически любому искаженному набору данных «значительно» лучше, чем один гауссиан .
whuber
Ты прав, @whuber. Я попытался сделать это явно. Я не уверен, как сделать гамма FMM, но это было бы лучше.
gung - Восстановить Монику
1
Поскольку это исследование, одной из мыслей было бы попытаться преобразовать исходное распределение в симметрию (возможно, со смещенным преобразованием Бокса-Кокса, надежно оцененным по нескольким квантилям данных) и повторить попытку. Конечно, вы не будете говорить о «значимости» как таковой, но анализ вероятности все еще может показывать.
whuber
@whuber, я сделал это, но упомянул об этом только мимоходом. (Оптимальным преобразованием Бокса-Кокса является обратный квадратный корень.) Вы получаете тот же результат, но значения p (все еще высокие, но) менее значимы.
gung - Восстановить Монику
3
Мне очень нравится идея, что вы должны моделировать то, что вы считаете процессом генерации. Моя проблема в том, что даже когда гауссовые смеси работают хорошо, я чувствую, что должна быть содержательная интерпретация. Если ФП расскажет нам больше, даже о том, что это за данные, возможно, будут более точные предположения.
Ник Кокс,
10

Следуя идеям в ответе и комментариях @ Nick, вы можете увидеть, насколько широкой должна быть просто сгладить вторичный режим:

введите описание изображения здесь

Возьмите эту оценку плотности ядра как проксимальный нуль - распределение, наиболее близкое к данным, но все еще согласующееся с нулевой гипотезой, что это образец из унимодальной популяции, - и смоделируйте его. В смоделированных выборках вторичный режим не часто выглядит настолько отчетливым, и вам не нужно сильно увеличивать полосу пропускания, чтобы сгладить ее.

<code> введите описание изображения здесь </ code>

Формализация этого подхода приводит к тесту, приведенному в Silverman (1981), «Использование оценок плотности ядра для исследования модальности», JRSS B , 43 , 1. silvermantestПакет Schwaiger & Holzmann реализует этот тест, а также процедуру калибровки, описанную Hall & York ( 2001), «О калибровке критерия Сильвермана для мультимодальности», Statistica Sinica , 11 , p 515, в котором учитывается асимптотический консерватизм. Выполнение теста на ваших данных с нулевой гипотезой унимодальности приводит к значениям р 0,08 без калибровки и 0,02 с калибровкой. Я не достаточно знаком с тестом на падение, чтобы догадаться, почему он может отличаться.

Код 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)
Scortchi - Восстановить Монику
источник
+1. Интересный! Какое ядро ​​используется здесь? Как я смутно помню, есть тонкие причины, по которым ядра Гаусса следует использовать для формальных вариантов этого подхода.
Ник Кокс
@Nick: Гауссово ядро, но я не могу вспомнить, есть ли веская причина для этого. Каждый моделируемый образец масштабируется, и есть поправка для консервативного смещения, которое было в оригинальном тесте - разработанном кем-то по имени Стори, я думаю.
Scortchi - Восстановить Монику
@NickCox: Извините, но не Стори.
Scortchi - Восстановить Монику
@ Scortchi, я немного подправил твой текст и код. Надеюсь, ты не против. +1. Кроме того, вы используете страшный оператор назначения стрелка вправо ?! О человечество ...
бандитский - Восстановить Монику
2
На самом деле это не лучше и не хуже, но в программировании принято указывать ваши переменные слева и назначать их справа. Многие люди взволнованы ->; Я просто ошеломлен.
gung - Восстановить Монику
7

Вещи для беспокойства включают в себя:

  1. Размер набора данных. Это не крошечный, не большой.

  2. Зависимость того, что вы видите от происхождения гистограммы и ширины бина. Поскольку очевиден только один выбор, вы (и мы) понятия не имеем о чувствительности.

  3. Зависимость того, что вы видите, от типа и ширины ядра и любых других вариантов, сделанных для вас при оценке плотности. Поскольку очевиден только один выбор, вы (и мы) понятия не имеем о чувствительности.

В другом месте я предположительно предположил, что достоверность режимов поддерживается (но не установлена) существенной интерпретацией и способностью различать ту же модальность в других наборах данных того же размера. (Чем больше, тем лучше ....)

Мы не можем комментировать ни одного из них здесь. Один небольшой критерий повторяемости заключается в сравнении того, что вы получаете с образцами начальной загрузки того же размера. Вот результаты токен-эксперимента с использованием Stata, но то , что вы видите, произвольно ограничено значениями по умолчанию Stata, которые сами задокументированы как вырванные из воздуха . Я получил оценки плотности для исходных данных и для 24 образцов начальной загрузки из того же.

Показатель (не больше, не меньше) - это то, что, как мне кажется, опытные аналитики могли бы просто угадать из вашего графика. Режим левой руки хорошо повторяется, а правая рука заметно более хрупкая.

Обратите внимание, что в этом есть неизбежность: поскольку в правом режиме меньше данных, они не всегда будут появляться в образце начальной загрузки. Но это также ключевой момент.

введите описание изображения здесь

Обратите внимание, что пункт 3. выше остается нетронутым. Но результаты находятся где-то между унимодальным и бимодальным.

Для тех, кто заинтересован, это код:

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' 
Ник Кокс
источник
+1 Мне нравится ваш подход начальной загрузки: массив графиков помогает всем лучше понять данные. Мне интересно, могут ли эти графики быть чувствительными к тому, как Stata оценивает пропускную способность. Я подозреваю, что это может привести к недостаточной мощности теста, потому что его оценка, вероятно, основана на одномодальном предположении, что приводит к относительно широкой полосе пропускания. Даже немного более узкая оценка пропускной способности может сделать второй режим более заметным во всех выборках начальной загрузки.
whuber
2
@ whuber Спасибо! Как обычно, вы безошибочно сосредотачиваетесь на слабостях, о которых нам нужно беспокоиться, и я согласен. По мере увеличения пропускной способности ядра появление унимодальности имеет тенденцию к неизбежности. И наоборот, небольшие полосы пропускания часто просто указывают на паразитные режимы, которые являются неповторимыми и / или тривиальными. Компромисс действительно деликатный. Я думаю, что главной заслугой этого подхода является риторика "Это воспроизводимо, если мы покачиваемся?" Меня часто беспокоит готовность пользователей программного обеспечения копировать результаты по умолчанию, не задумываясь.
Ник Кокс
2
Существуют систематические подходы к этой проблеме, основанные на постепенном изменении полосы пропускания и отслеживании появления и исчезновения режимов при изменении полосы пропускания. По сути, надежный режим сохраняется, а менее надежный - нет. Это симпатичный подход, но иногда он запускает конструктор туннеля, когда подойдет лопата. Например, если вы изменяете выбор гистограммы и вторичный режим слишком быстро исчезает (или перемещается), не верьте этому.
Ник Кокс
2

LP Непараметрическая идентификация режима

LP Непараметрическая идентификация режима (название алгоритма LPMode , ссылка на статью приведена ниже)

MaxEnt Modes [Красные цветные треугольники на графике]: 12783.36 и 24654.28.

Режимы L2 [Зеленые цветные треугольники на графике]: 13054.70 и 24111.61.

Интересно отметить модальные формы, особенно вторую, которая демонстрирует значительную асимметрию (модель традиционной гауссовой смеси, скорее всего, потерпит неудачу).

Мухопадхяй, С. (2016). Идентификация в крупном масштабе и управляемые данными науки. https://arxiv.org/abs/1509.06428

Глубокое Мукерджи
источник
1
Можете ли вы разработать и предоставить некоторый контекст, чтобы представить и объяснить эти методы? Приятно иметь ссылку на статью, но мы предпочитаем, чтобы наши ответы здесь были автономными, особенно если ссылка не работает.
gung - Восстановить Монику
Контекст является исходным вопросом: существует ли мультимодальность? если так позиции. и актуальность нового метода проистекает из того факта, что охота на неровности непараметрическим способом является сложной задачей моделирования.
Глубокий Мукерджи,
@gung просит вас расширить свой ответ. Например, результатом является метод, описанный в статье, которая не имеет общедоступной версии.
Ник Кокс
2
Нет, я имею в виду, что такое "Идентификация непараметрического режима LP"? Что такое "MaxEnt"? И т.д. В нескольких предложениях, как это работает? Почему / когда это может быть предпочтительнее других методов? И т. Д. Я знаю, что вы ссылаетесь на статью, которая объясняет их, но было бы неплохо иметь пару предложений, чтобы представить их здесь, особенно если ссылка не работает, но даже если не дать будущим читателям понять, будут ли они хочу продолжить этот метод.
gung - Восстановить Монику
2
@DeepMukherjee, вам, конечно, не нужно переписывать всю статью в своем посте. Просто добавьте несколько предложений о том, что это такое и как это работает.
gung - Восстановить Монику