Calcolare le modalità in una distribuzione multimodale in R

Question

Ho misurato le altezze del corpo di tutti i miei figli. Quando ho tracciare tutte le altezze lungo un asse di lunghezze, questo è ciò che il risultato appare come:

Ogni rossa (ragazzi) o violetto (ragazze) tick è un bambino. Se due bambini hanno la stessa altezza del corpo (in millimetri), le zecche si impilano. Attualmente ci sono sette bambini con la stessa altezza del corpo. (L'altezza e la larghezza delle zecche non ha significato. Sono state ridimensionate per essere visibili.)

Come si può vedere, le diverse altezze non sono distribuite uniformemente lungo l'asse, ma si raggruppano intorno a determinati valori.

Una trama Istogramma e la densità dei dati si presenta così (con le due stime di densità tracciati seguenti this answer):

Come si può vedere, si tratta di una distribuzione multimodale.

Come si calcolano le modalità (in R)?

---

Ecco i dati grezzi per voi a giocare con:

mm <- c(418, 527, 540, 553, 554, 558, 613, 630, 634, 636, 645, 648, 708, 714, 715, 725, 806, 807, 822, 823, 836, 837, 855, 903, 908, 910, 911, 913, 915, 923, 935, 945, 955, 957, 958, 1003, 1006, 1015, 1021, 1021, 1022, 1034, 1043, 1048, 1051, 1054, 1058, 1100, 1102, 1103, 1117, 1125, 1134, 1138, 1145, 1146, 1150, 1152, 1210, 1211, 1213, 1223, 1226, 1334) 

Answer 1

ho costruito qualcosa per conto mio usando i dati mm.

Prima di tutto tracciare la densità del mm, al fine di visualizzare i modi:

plot(density(mm)) 

Quindi, possiamo vedere ci sono 2 modi in questa distribuzione. Uno intorno al 600 e uno intorno al 1000. Vediamo come trovarli.

Al fine di trovare gli indici della modalità che ho fatto questa funzione:

find_modes<- function(x) { 
    modes <- NULL 
    for (i in 2:(length(x)-1)){ 
    if ((x[i] > x[i-1]) & (x[i] > x[i+1])) { 
     modes <- c(modes,i) 
    } 
    } 
    if (length(modes) == 0) { 
    modes = 'This is a monotonic distribution' 
    } 
    return(modes) 
} 

Proviamo sulla nostra densità:

mymodes_indices <- find_modes(density(mm)$y) #you need to try it on the y axis 

Ora mymodes_indices contiene gli indici dei nostri modi vale a dire:

> density(mm)$y[mymodes_indices] #just to confirm that those are the correct 
[1] 0.0008946929 0.0017766183 

> density(mm)$x[mymodes_indices] #the actual modes 
[1] 660.2941 1024.9067 

Spero che aiuti!

Answer 2

Ho modificato la risposta di Jeffrey Evans in Peak of the kernel density estimation per consentire di modificare il parametro bw e di conseguenza ottenere più o meno picchi. È necessario per altri casi, che produrranno molti picchi con la risposta accettata. Il parametro signifi consente di gestire le fascette.

library(dplyr) 
library(tidyr) 
get.modes2 <- function(x,bw,signifi) { 
    den <- density(x, kernel=c("gaussian"),bw=bw) 
    den.s <- smooth.spline(den$x, den$y, all.knots=TRUE, spar=0.1) 
    s.1 <- predict(den.s, den.s$x, deriv=1) 
    s.0 <- predict(den.s, den.s$x, deriv=0) 
    den.sign <- sign(s.1$y) 
    a<-c(1,1+which(diff(den.sign)!=0)) 
    b<-rle(den.sign)$values 
    df<-data.frame(a,b) 
    df = df[which(df$b %in% -1),] 
    modes<-s.1$x[df$a] 
    density<-s.0$y[df$a] 
    df2<-data.frame(modes,density) 
    df2$sig<-signif(df2$density,signifi) 
    df2<-df2[with(df2, order(-sig)), ] 
    print(df2) 
    df<-as.data.frame(df2 %>% 
         mutate(m = min_rank(desc(sig))) %>% #, count = sum(n)) %>% 
         group_by(m) %>% 
         summarize(a = paste(format(round(modes,2),nsmall=2), collapse = ',')) %>% 
         spread(m, a, sep = '')) 
    colnames(df)<-paste0("m",1:length(colnames(df))) 
    print(df) 
} 
mm <- c(418, 527, 540, 553, 554, 558, 613, 630, 634, 636, 645, 648, 708, 714, 715, 725, 806, 807, 822, 823, 836, 837, 855, 903, 908, 910, 911, 913, 915, 923, 935, 945, 955, 957, 958, 1003, 1006, 1015, 1021, 1021, 1022, 1034, 1043, 1048, 1051, 1054, 1058, 1100, 1102, 1103, 1117, 1125, 1134, 1138, 1145, 1146, 1150, 1152, 1210, 1211, 1213, 1223, 1226, 1334) 
mmdf<-data.frame(mm=mm) 
library(ggplot2) 
#0.25 defines the number of peaks. 
ggplot(mmdf,aes(mm)) + geom_density(adjust=0.25) 
modes<-get.modes2(mm,bw.nrd0(mm)*0.25,2) 
     m1  m2  m3   m4  m5  m6  m7    m8 
1 1032.08 921.91 1134.09 636.27,826.01 1217.74 548.54 715.84 420.00,1334.04