Bagaimana cara menguji apakah distribusi saya multimodal?

Saat saya memetakan histogram data saya, ini memiliki dua puncak:

Apakah itu berarti distribusi multi-modal yang potensial? Saya menjalankan dip.testdi R ( library(diptest)), dan hasilnya adalah:

D = 0.0275, p-value = 0.7913

Saya dapat menyimpulkan bahwa data saya memiliki distribusi multi-modal?

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

@NickCox telah menyajikan strategi yang menarik (+1). Saya mungkin mempertimbangkan lebih eksplorasi di alam Namun, karena kekhawatiran bahwa @whuber menunjukkan .

Izinkan saya menyarankan strategi lain: Anda dapat menggunakan model campuran hingga Gaussian. Perhatikan bahwa ini membuat asumsi yang sangat kuat bahwa data Anda diambil dari satu atau lebih normals sejati. Seperti yang ditunjukkan oleh @whuber dan @NickCox dalam komentar, tanpa interpretasi substantif dari data ini — didukung oleh teori yang sudah mapan — untuk mendukung asumsi ini, strategi ini juga harus dipertimbangkan sebagai eksplorasi.

Pertama, mari ikuti saran @ Glen_b dan lihat data Anda menggunakan nampan dua kali lebih banyak:

Kami masih melihat dua mode; jika ada, mereka datang dengan lebih jelas di sini. (Perhatikan juga bahwa garis kepadatan kernel harus identik, tetapi nampak lebih tersebar karena jumlah sampah yang lebih besar.)

Sekarang mari kita muat model campuran hingga Gaussian. Di R, Anda dapat menggunakan Mclustpaket untuk melakukan ini:

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 

Dua komponen normal mengoptimalkan BIC. Sebagai perbandingan, kami dapat memaksakan kecocokan satu komponen dan melakukan uji rasio kemungkinan:

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

Ini menunjukkan sangat tidak mungkin Anda akan menemukan data yang jauh dari unimodal seperti milik Anda jika mereka berasal dari satu distribusi normal sejati tunggal.

Beberapa orang merasa tidak nyaman menggunakan tes parametrik di sini (walaupun jika asumsi itu berlaku, saya tidak tahu masalah apa pun). Salah satu teknik yang sangat luas diterapkan adalah dengan menggunakan Metode Parametrik Bootstrap Cross-fitting (saya jelaskan algoritma di sini ). Kami dapat mencoba menerapkannya pada data ini:

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 

Statistik ringkasan, dan plot kerapatan kernel untuk distribusi sampling menunjukkan beberapa fitur menarik. Kemungkinan log untuk model komponen tunggal jarang lebih besar daripada kecocokan dua komponen, bahkan ketika proses menghasilkan data yang benar hanya memiliki satu komponen tunggal, dan ketika lebih besar, jumlahnya sepele. Gagasan membandingkan model yang berbeda dalam kemampuan mereka untuk mencocokkan data adalah salah satu motivasi di balik PBCM. Dua distribusi sampling hampir tidak tumpang tindih sama sekali; hanya 0,35% dari x2.dyang kurang dari maksimumx1.dnilai. Jika Anda memilih model dua komponen jika perbedaan dalam kemungkinan log> 9.7, Anda akan salah memilih model satu komponen .01% dan dua model komponen .02% dari waktu. Ini sangat diskriminatif. Jika, di sisi lain, Anda memilih untuk menggunakan model satu komponen sebagai hipotesis nol, hasil yang Anda amati cukup kecil sehingga tidak muncul dalam distribusi sampel empiris dalam 10.000 iterasi. Kita dapat menggunakan aturan 3 (lihat di sini ) untuk menempatkan batas atas pada nilai-p, yaitu, kami memperkirakan nilai p Anda kurang dari .0003. Artinya, ini sangat signifikan.

$p < .000001$$p < .001$), dan komponen yang mendasarinya, jika ada, juga tidak dijamin normal. Jika menurut Anda masuk akal bahwa data Anda dapat berasal dari distribusi yang condong positif, dan bukan normal, tingkat bimodalitas ini mungkin berada dalam kisaran variasi yang khas, yang saya curigai katakan dalam tes celup.

Menindaklanjuti ide-ide dalam jawaban dan komentar @ Nick, Anda dapat melihat seberapa lebar bandwidth yang diperlukan untuk hanya meratakan mode sekunder:

Ambil estimasi kerapatan kernel ini sebagai nol proksimal — distribusi yang paling dekat dengan data namun masih konsisten dengan hipotesis nol bahwa itu adalah sampel dari populasi unimodal — dan disimulasikan darinya. Dalam sampel yang disimulasikan, mode sekunder tidak sering terlihat begitu berbeda, dan Anda tidak perlu memperluas lebar pita untuk meratakannya.

Meresmikan pendekatan ini mengarah ke tes yang diberikan dalam Silverman (1981), "Menggunakan perkiraan kepadatan kernel untuk menyelidiki modalitas", JRSS B , 43 , 1. Schwaiger & Holzmann ini silvermantestpaket alat tes ini, dan juga prosedur kalibrasi dijelaskan oleh Hall & York ( 2001), "Pada kalibrasi uji Silverman untuk multimodality", Statistica Sinica , 11 , p 515, yang menyesuaikan konservatisme asimptotik. Melakukan tes pada data Anda dengan hipotesis nol hasil unimodality dalam nilai-p 0,08 tanpa kalibrasi dan 0,02 dengan kalibrasi. Saya tidak cukup terbiasa dengan tes celup untuk menebak mengapa itu mungkin berbeda.

Kode 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)

Hal-hal yang perlu dikhawatirkan antara lain:

1. Ukuran dataset. Itu tidak kecil, tidak besar.

2. Ketergantungan apa yang Anda lihat pada asal histogram dan lebar nampan. Dengan hanya satu pilihan yang jelas, Anda (dan kami) tidak memiliki gagasan tentang sensitivitas.

3. Ketergantungan apa yang Anda lihat pada jenis dan lebar kernel dan apa pun pilihan lain yang dibuat untuk Anda dalam estimasi kepadatan. Dengan hanya satu pilihan yang jelas, Anda (dan kami) tidak memiliki gagasan tentang sensitivitas.

Di tempat lain saya menyarankan secara sementara bahwa kredibilitas mode didukung (tetapi tidak ditetapkan) oleh interpretasi substantif dan oleh kemampuan untuk membedakan modalitas yang sama dalam dataset lain dengan ukuran yang sama. (Lebih besar lebih baik juga ....)

Kami tidak dapat mengomentari salah satu dari mereka di sini. Satu pegangan kecil pada pengulangan adalah untuk membandingkan apa yang Anda dapatkan dengan sampel bootstrap dengan ukuran yang sama. Berikut ini adalah hasil percobaan token menggunakan Stata, tetapi apa yang Anda lihat dibatasi secara sewenang-wenang pada standar Stata, yang didokumentasikan sebagai dicabut dari udara . Saya mendapat estimasi kepadatan untuk data asli dan untuk 24 sampel bootstrap dari yang sama.

Indikasinya (tidak lebih, tidak kurang) adalah apa yang saya pikir analis berpengalaman hanya akan menebak dengan cara apa pun dari grafik Anda. Mode tangan kiri sangat berulang dan tangan kanan jelas lebih rapuh.

Perhatikan bahwa ada hal yang tak terhindarkan tentang ini: karena ada lebih sedikit data di dekat mode kanan, itu tidak akan selalu muncul kembali dalam sampel bootstrap. Tetapi ini juga merupakan poin kunci.

Perhatikan bahwa titik 3. di atas tetap tidak tersentuh. Tetapi hasilnya ada di suatu tempat antara unimodal dan bimodal.

Bagi yang berminat, ini kodenya:

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' 

Identifikasi Mode Nonparametrik LP (nama algoritma LPMode , referensi makalah diberikan di bawah ini)

Mode MaxEnt [Segitiga warna merah dalam plot]: 12783.36 dan 24654.28.

Mode L2 [Segitiga warna hijau dalam plot]: 13054.70 dan 24111.61.

Menarik untuk dicatat bentuk modal, terutama yang kedua yang menunjukkan kemiringan yang cukup besar (model Campuran Gaussian Tradisional cenderung gagal di sini).

Mukhopadhyay, S. (2016) Identifikasi Mode Skala Besar dan Ilmu Pengetahuan Berbasis Data. https://arxiv.org/abs/1509.06428