Viacvrcholové rozdelenia sú častým prípadom rozdelenia početností ekonomických javov. Poukazujú na na nehomogénny súbor údajov. Každá vrcholová oblasť a jeho blízke okolie charakterizuje niečo typické.


Poznámky


Stáva sa, že na základe histogramu alebo jadrovej funkcie hustoty sa nevieme rozhodnúť, či sa jedná o jednovrcholové rozdelenie alebo dvojvrcholové, resp. viacvrcholové pravdepodobnostné rozdelenie. Ako k riešeniu problému pristupovať si ukážeme na príklade.

Príklad
Budeme pracovať s údajmi o cene práce (pre odvetvia Priemysel, stavebníctvo a služby (okrem verejnej správy, obrany, povinného sociálneho zabezpečenia)) za rok 2008 a 2020 z Eurostatu.

Tabuľka. Cena práce za krajiny EÚ v rokoch 2008 a 2020 (v Eurách).

Zdroj údajov:https://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=lc_lci_lev&lang=en

Načítanie údajov
V prvom kroku načítame údaje do datovej tabuľky udajeEU. Budeme potrebovať package “xlsx”. Údaje sú v prvom hárku súboru cenaEU.xlsx, ktorý sa nachádza adresári bs na disku c.

library(xlsx)
udajeEU = read.xlsx("c:/bs/cenaEU.xlsx",1) 
attach(udajeEU)


Načítajme názvy štátov

rownames(udajeEU) = c("BE", "BG", "CZ", "DK", "DE", "EE", "IE", "EL", "ES", "FR", "HR", "IT", "CY", "LV", "LT", "LU", "HU","MT", "NL", "AT", "PL", "PT", "RO", "SI", "SK", "FI", "SE")


Vypíšme datovú tabuľku udajeEU. 

udajeEU
##    Rok2008 Rok2020
## BE    32.9    41.1
## BG     2.6     6.5
## CZ     9.2    14.1
## DK    34.6    45.8
## DE    27.9    36.6
## EE     7.9    13.6
## IE    28.9    32.3
## EL    16.8    16.9
## ES    19.4    22.8
## FR    31.2    37.5
## HR     9.2    10.8
## IT    25.2    29.8
## CY    16.7    17.0
## LV     5.9    10.5
## LT     5.9    10.1
## LU    32.3    42.1
## HU     7.8     9.9
## MT    11.4    14.5
## NL    29.8    36.8
## AT    26.4    36.7
## PL     7.6    11.0
## PT    12.2    15.7
## RO     4.2     8.1
## SI    13.9    19.9
## SK     7.0    13.4
## FI    27.1    34.3
## SE    31.6    37.3


Histogram a jadrové funkcie hustoty
Pomocou nasledujúcich príkazov vykreslíme histogram a jadrové funkcie hustoty. 

hist(udajeEU$Rok2008, # histogram   
     col="gray", # sivá farba stĺpcov v histograme  
     border="black",    # rám stĺpcov v histograme čiernou farbou    
     prob = TRUE, # na osi y sa znázornia relatívne početnosti a nie početnosti ( pre tie je hodnota prob=FALSE)    
     xlab = "cena práce (v Eurách)", # popis osi x  
     xlim=c(0,50), # hodnoty na osi x budú od 0 do 50 (v pripade potreby porovnať posun je vhodné mať zhodné hodnoty na osi x)
     main = "Cena práce v krajinách EÚ v roku 2008")

lines(density(udajeEU$Rok2008), # vykreslí funkciu hustoty
     lwd = 2, # hrúbka čiary    
     col = "green") # farba čiary


hist(udajeEU$Rok2020, # histogram   
     col="gray", # sivá farba stĺpcov v histograme  
     border="black",    # rám stĺpcov v histograme čiernou farbou    
     prob = TRUE, # na osi y sa znázornia relatívne početnosti a nie početnosti ( pre tie je hodnota prob=FALSE)    
     xlab = "cena práce (v Eurách)", # popis osi x  
     xlim=c(0,50), # hodnoty na osi x budú od 0 do 50 (v pripade potreby porovnať posun je vhodné mať zhodné hodnoty na osi x)
     main = "Cena práce v krajinách EÚ v roku 2020")


lines(density(udajeEU$Rok2020), # vykreslí funkciu hustoty
     lwd = 2, # hrúbka čiary    
     col = "green") # farba čiary


Overovanie hypotézy o viacvrcholovom rozdelení
Nulová hypotéza pre Hartiganov test (Hartigan a Hartigan, 1985) je Ho: rozdelenie je jednovrcholové a alternatívna hypotéza je, že nie je jednovrcholové; aspoň dvojvrcholové. 

library(diptest)
dip.test(udajeEU$Rok2008)  # testovanie hypotézy, že rozdelenie údajov za rok 2008 je jednovrcholové a alternatívna hypotéza je, že nie je jednovrcholové pomocou Hartiganovho testu
## 
##  Hartigans' dip test for unimodality / multimodality
## 
## data:  udajeEU$Rok2008
## D = 0.08892, p-value = 0.07352
## alternative hypothesis: non-unimodal, i.e., at least bimodal
dip.test(udajeEU$Rok2020)  #testovanie hypotézy, že rozdelenie údajov za rok 2020 je jednovrcholové a alternatívna hypotéza je, že nie je jednovrcholové pomocou Hartiganovho testu
## 
##  Hartigans' dip test for unimodality / multimodality
## 
## data:  udajeEU$Rok2020
## D = 0.084463, p-value = 0.1097
## alternative hypothesis: non-unimodal, i.e., at least bimodal


P hodnota pre cenu prace v roku 2008 je 0,07532. Znamená to, že na hladine významnosti 0,05 nemôžeme nulovú hypotézu o jednovrcholovom rozdelení zamietnuť. Na hladine významnosti 0,10 však môžeme nulovú hypotézu zamietnuť. Rozdelenie nie je jednovrcholové. To znamená, že existujúu dve skupiny štátov z hľadiska výšky ceny práce.

P hodnota pre cenu prace v roku 2020 je 0,1097. Znamená to, že na hladine významnosti 0,05 a ani na hladine významnosti 0,10 nemôžeme nulovú hypotézu o jednovrcholovom rozdelení zamietnuť. Pohľad na funkciu hustoty však vyvoláva pochybnosti. Okrem Hartiganovho testu existuju aj ďalšie testy.

Knižnica multimode obsahuje okrem Hartiganovho testu HH aj ďalšie Silvermanov SI (Silverman, 1981), Hallov a Yarkov HY (Hall and York, 2001), Fisherov a Marronov FM (Fisher and Marron, 2001), Chengov a Hallov CH (Cheng and Hall,1998), Ameijeirasov-Alonsov ACR (Ameijeiras-Alonso et al., 2019). Pozrime si vysledky pre cenu prace za rok 2020 pomocou nich. 

library(multimode)
modetest(udajeEU$Rok2020,mod0=1,method="SI") # Silvermanov test, že rozdelenie má jeden vrchol (mod0=1)
## 
##  Silverman (1981) critical bandwidth test
## 
## data:  udajeEU$Rok2020
## Critical bandwidth = 8.9006, p-value = 0.3
## alternative hypothesis: true number of modes is greater than 1
modetest(udajeEU$Rok2020,mod0=1,method="HY") # Hallov a Yarkov test, že rozdelenie má jeden vrchol (mod0=1)
## Warning in modetest(udajeEU$Rok2020, mod0 = 1, method = "HY"): Arguments
## 'lowsup' and 'uppsup' should be finite unless the density function has a bounded
## support
## 
##  Hall and York (2001) critical bandwidth test
## 
## data:  udajeEU$Rok2020
## Critical bandwidth = 8.9006, p-value = 0.076
## alternative hypothesis: true number of modes is greater than 1
modetest(udajeEU$Rok2020,mod0=1,method="FM") # Fisherov a Marronov test, že rozdelenie má jeden vrchol (mod0=1)
## 
##  Fisher and Marron (2001) Cramer-von Mises test
## 
## data:  udajeEU$Rok2020
## Cramer-von Mises = 0.20988, p-value = 0.006
## alternative hypothesis: true number of modes is greater than 1
modetest(udajeEU$Rok2020,mod0=1,method="CH") # Chengov a Hallov, že rozdelenie má jeden vrchol (mod0=1)
## Warning in excessmassex(data, mod0): A modification of the data was made in
## order to compute the excess mass statistic
## 
##  Cheng and Hall (1998) excess mass test
## 
## data:  udajeEU$Rok2020
## Excess mass = 0.17004, p-value = 0.004
## alternative hypothesis: true number of modes is greater than 1
modetest(udajeEU$Rok2020,mod0=1,method="ACR") # Ameijeirasov-Alonsov test, že rozdelenie má jeden vrchol (mod0=1)
## Warning in excessmassex(data, mod0): A modification of the data was made in
## order to compute the excess mass statistic
## 
##  Ameijeiras-Alonso et al. (2019) excess mass test
## 
## data:  udajeEU$Rok2020
## Excess mass = 0.16759, p-value = 0.046
## alternative hypothesis: true number of modes is greater than 1


Naozaj, možeme urobiť záver, že cena práce v krajiných EÚ má dvojvrcholové rozdelenie. Nájdime tieto vrcholy pomocou bootstrapingu. 

locmodes(udajeEU$Rok2008,mod0=2,display=TRUE)  # určenie vrcholov v pripade, že počet vrcholov je 2 (mod0=2)
## Warning in locmodes(udajeEU$Rok2008, mod0 = 2, display = TRUE): If the density
## function has an unbounded support, artificial modes may have been created in the
## tails

## 
## Estimated location
## Modes: 7.556531  30.55259 
## Antimode:  21.9193 
## 
## Estimated value of the density
## Modes: 0.05373277  0.04102266 
## Antimode: 0.00689931 
## 
## Critical bandwidth:  2.02202
locmodes(udajeEU$Rok2020,mod0=2,display=TRUE)  # určenie vrcholov v pripade, že počet vrcholov je 2 (mod0=2)
## Warning in locmodes(udajeEU$Rok2020, mod0 = 2, display = TRUE): If the density
## function has an unbounded support, artificial modes may have been created in the
## tails

## 
## Estimated location
## Modes: 10.89206  36.80524 
## Antimode: 26.28913 
## 
## Estimated value of the density
## Modes: 0.05098975  0.04494425 
## Antimode: 0.002336202 
## 
## Critical bandwidth: 1.753967


Vrcholy v prípade ceny práce v roku 2008 sú 7,556531 a 30,55259. Hodnota funkcie hustoty v prvom vrchole jw 0,05373277 a v prípade druhého vrcholu 0,04102266.
Vrcholy v prípade ceny práce v roku 2020 sú 10,89206 a 36,80524. Hodnota funkcie hustoty v prvom vrchole jw 0,05098975 a v prípade druhého vrcholu 0,04494425.
Hľadané vrcholy sú 10,89206 a 36,80524. Zodpovedá to členeniu krajín Belgicko, Dánsko, Francúzsko, Nemecko, Írsko, Taliansko, Luxembursko, Holandsko a ostatné krajiny. V prípade roku 2008 sú vrcholy 7,556531 a 30,55259. Pre obe skupiny krajín došlo medzi rokmi 2008 a 2020 k nárastu ceny práce. Hodnota antimodusu (bod s najnižšou relatívnou početnosťou) sa z hodnoty 21,91933 posunula na hodnotu 26,28913.

Zmes normálnych rozdelení
Keby sme predpokladali, že výsledné pravdepodobnostné rozdelenie je zmesou dvoch normálnych rozdelení, môžeme použiť knižnicu mixtools. 

library(mixtools)
## mixtools package, version 1.2.0, Released 2020-02-05
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.
c2008 <- normalmixEM(udajeEU$Rok2008, k=2, lambda=NULL,mu = NULL, sigma = NULL, maxit=50, epsilon=0.01)  # počet vrcholov je 2 (k=2).Hodnota sa ignoruje, pokiaľ mu aj sigma nie sú nulové (NULL); počiatočná hodnota pomerov miešania (lambda) je nula,  maximálny počet iterácií je 50 (maxit=50), konvergenčné kritérium epsilon=0,01)
## number of iterations= 26
summary(c2008) # vypíše hodnoty lambda, mu, sigma
## summary of normalmixEM object:
##          comp 1    comp 2
## lambda 0.592948  0.407052
## mu     9.870194 29.805457
## sigma  4.694635  2.853055
## loglik at estimate:  -92.64582
plot(c2008, density=TRUE) # vykreslí  histogram a odhadnuté funkcie hustoty

c2020 <- normalmixEM(udajeEU$Rok2020, k=2, lambda=NULL,mu = NULL, sigma = NULL, maxit=50, epsilon=0.01) 
## number of iterations= 2
summary(c2020)
## summary of normalmixEM object:
##           comp 1    comp 2
## lambda  0.591735  0.408265
## mu     13.413666 37.268257
## sigma   4.195649  4.329067
## loglik at estimate:  -95.48336
plot(c2020, density=TRUE)



Lambda je váha jednotlivých komponentov, t.j. f(x) = λ1f1(x) + λ2f2(x).
Pre rok 2008 f1 je funkcia hustoty N(9,865;4,691) a f2 je funkcia hustoty N(29,800504;2,857802), λ1 =0,592793, λ2 =0,407207 .
Pre rok 2020 f1 je funkcia hustoty N(13,40618;4,18988) a f2 je funkcia hustoty N(37,25241;4,34822), λ1 = 0,59137, λ2 = 0,40863.

Literatúra

R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Adrian Dragulescu and Cole Arendt (2020). xlsx: Read, Write, Format Excel 2007 and Excel 97/2000/XP/2003 Files. R package version 0.6.5. https://CRAN.R-project.org/package=xlsx

Martin Maechler (2021). diptest: Hartigan’s Dip Test Statistic for Unimodality - Corrected. R package version 0.76-0. https://CRAN.R-project.org/package=diptest

Ameijeiras-Alonso J, Crujeiras RM, Rodr'iguez-Casal A (2021). “multimode: An R Package for Mode Assessment.” Journal of Statistical Software, 97(9), 1-32. doi: 10.18637/jss.v097.i09 (URL: https://doi.org/10.18637/jss.v097.i09).

Ameijeiras-Alonso, J., Crujeiras, R.M., Rodríguez-Casal, A. (2021). multimode: An R Package for Mode Assessment, Journal of Statistical Software, 97, 1–32.

Cheng, M. Y., Hall, P. (1998). Calibrating the excess mass and dip tests of modality, Journal of the Royal Statistical Society. Series B, 60, 579–589.

Fisher, N.I., Marron, J. S. (2001). Mode testing via the excess mass estimate, Biometrika, 88, 419–517.

Hall, P., York, M. (2001). On the calibration of Silverman’s test for multimodality, Statistica Sinica, 11, 515–536.

Hartigan, J. A., Hartigan, P. M. (1985). The Dip Test of Unimodality, Journal of the American Statistical Association, 86, 738–746.

Hartigan, J. A.,Hartigan, P. M. (1985). The dip test of unimodality. The annals of Statistics, 70-84.

Silverman, B. W. (1981). Using kernel density estimates to investigate multimodality, Journal of the Royal Statistical Society. Series B, 43, 97–99.

Xu, L., Bedrick, E. J., Hanson, T., Restrepo, C. (2014). A comparison of statistical tools for identifying modality in body mass distributions. Journal of data science, 12(1), 175-196.