istogramma ggplot2 con curva di densità che somma a 1

Question

Tracciare un istogramma con una curva di densità che corrisponde a 1 per i dati non standardizzati è ridicolmente difficile. Ci sono già molte domande su questo, ma nessuna delle loro soluzioni funziona per i miei dati. Deve esserci una soluzione semplice che funziona e basta. Non riesco a trovare una risposta con una soluzione semplice che funzioni.

Alcuni esempi:

soluzione funziona solo con dati normali standardizzate ggplot2: Overlay histogram with density curve

con dati discreti e nessuna curva di densità ggplot2 density histogram with width=.5, vline and centered bar positions

risposta Overlay density and histogram plot with ggplot2 using custom bins

densità non ammontino al 1 sui miei dati Creating a density histogram in ggplot2?

non somma da 1 a miei dati ggplot2 density histogram with custom bin edges

lunga spiegazione qui con gli esempi, ma la densità non è uno con i miei dati "Density" curve overlay on histogram where vertical axis is frequency (aka count) or relative frequency?

-

qualche esempio di codice:

#Example code 
set.seed(1) 
t = data.frame(r = runif(100)) 

#first we try the obvious simple solution that should work 
ggplot(t, aes(r)) + 
    geom_histogram() + 
    geom_density() 

Quindi, chiaramente la densità non somma 1.

#maybe geom_histogram needs a ..density.. ? 
ggplot(t, aes(r)) + 
    geom_histogram(aes(y = ..density..)) + 
    geom_density() 

Lo ha fatto cambiare qualcosa, ma non in modo corretto.

#maybe geom_density needs a ..density.. too ? 
ggplot(t, aes(r)) + 
    geom_histogram(aes(y = ..density..)) + 
    geom_density(aes(y = ..density..)) 

Nessun cambiamento lì.

#maybe binwidth = 1? 
ggplot(t, aes(r)) + 
    geom_histogram(aes(y = ..density..), binwidth=1) + 
    geom_density(aes(y = ..density..)) 

curva di densità ancora sbagliata, ma ora l'istogramma è sbagliato troppo.

Per essere sicuro, ho trascorso 4 ore a provare tutti i tipi di combinazioni di ..count .. e ..sum .. e ..densità .., ma poiché non riesco a trovare alcuna documentazione su come questi sono dovrebbe funzionare, è prova ed errore semi-cieco.

Quindi ho rinunciato e ho evitato di utilizzare ggplot2 per riepilogare i dati.

Quindi, prima abbiamo bisogno di ottenere le giuste proporzioni data.frame, e che non era così semplice:

get_prop_table = function(x, breaks_=20){ 
    library(magrittr) 
    library(plyr) 
    x_prop_table = cut(x, 20) %>% table(.) %>% prop.table %>% data.frame 
    colnames(x_prop_table) = c("interval", "density") 
    intervals = x_prop_table$interval %>% as.character 
    fetch_numbers = str_extract_all(intervals, "\\d\\.\\d*") 
    x_prop_table$means = laply(fetch_numbers, function(x) { 
    x %>% as.numeric %>% mean 
    }) 
    return(x_prop_table) 
} 

t_df = get_prop_table(t$r) 

Questo dà il tipo di dati di sintesi che vogliamo:

> head(t_df) 
      interval density means 
1 (0.00859,0.0585] 0.06 0.033545 
2 (0.0585,0.107] 0.09 0.082750 
3 (0.107,0.156] 0.07 0.131500 
4 (0.156,0.205] 0.10 0.180500 
5 (0.205,0.254] 0.08 0.229500 
6 (0.254,0.303] 0.03 0.278500 

Ora dobbiamo solo tracciarlo. Dovrebbe essere facile ...

ggplot(t_df, aes(means, density)) + 
    geom_histogram(stat = "identity") + 
    geom_density(stat = "identity") 

Umm, non proprio quello che volevo. Per essere sicuro, ho provato senza stat = "identity" in geom_density, a quel punto si è lamentato di non averne uno.

#lets try adding ..density.. then 
ggplot(t_df, aes(means, density)) + 
    geom_histogram(stat = "identity") + 
    geom_density(aes(y = ..density..)) 

Ancora più strano.

Ok, forse rinunciamo a ottenere la curva di densità dai dati di riepilogo. Forse abbiamo bisogno di mescolare gli approcci un po '...

#adding together 
ggplot(t_df, aes(means, density)) + 
    geom_bar(stat = "identity") + 
    geom_density(data=t, aes(r, y = ..density..), stat = 'density') 

Ok, almeno la forma è in questo momento. Ora, abbiamo bisogno di ridimensionarlo in qualche modo.

#lets try dividing by the number of bins 
ggplot(t_df, aes(means, density)) + 
    geom_bar(stat = "identity") + 
    geom_density(data=t, aes(r, y = ..density../20), stat = 'density') 

Sembra che abbiamo un vincitore. Tranne che il numero è hardcoded.

#removing the hardcoding? 
divisor = nrow(t_df) 
ggplot(t_df, aes(means, density)) + 
    geom_bar(stat = "identity") + 
    geom_density(data=t, aes(r, y = ..density../divisor), stat = 'density') 

Error in eval(expr, envir, enclos) : object 'divisor' not found 

Bene, quasi mi aspettavo che funzionasse. Ora ho provato ad aggiungere un po '.. qua e là, anche ..count .. and ..sum .., il primo che ha dato un altro risultato sbagliato, il secondo che ha generato un errore. Ho anche provato a usare un moltiplicatore (con 1/20), senza fortuna.

#salvation with get() 
divisor = nrow(t_df) 
ggplot(t_df, aes(means, density)) + 
    geom_bar(stat = "identity") + 
    geom_density(data=t, aes(r, y = ..density../get("divisor", pos = 1)), stat = 'density') 

Così, ho finalmente avuto la figura a destra (credo, spero).

Per favore dimmi che c'è un modo più semplice per farlo.

PS. Il trucco get() apparentemente non funziona all'interno di una funzione. Avrei messo una funzione di lavoro qui per un uso futuro, ma non è stato altrettanto facile.

Answer 1

Innanzitutto, leggere Wickham sulla densità in R, osservando le debolezze e le caratteristiche di ciascun pacchetto/funzione.

Le densità sommano a 1, ma ciò non significa che la linea/punti della curva non supererà 1.

Quanto segue mostra sia presente e l'imprecisione di (almeno) i valori di default density rispetto, per esempio, KernSmooth::bkde (utilizzando piazzole di base per brevità di battitura):

library(KernSmooth) 
library(flux) 
library(sfsmisc) 

# uniform dist 
set.seed(1) 
dat <- runif(100) 

d1 <- density(dat) 
d1_ks <- bkde(dat) 

par(mfrow=c(2,1)) 
plot(d1) 
plot(d1_ks, type="l") 

auc(d1$x, d1$y) 
## [1] 1.000921 

integrate.xy(d1$x, d1$y) 
## [1] 1.000921 

auc(d1_ks$x, d1_ks$y) 
## [1] 1 

integrate.xy(d1_ks$x, d1_ks$y) 
## [1] 1 

fare lo stesso per la distribuzione beta:

# beta dist 
set.seed(1) 
dat <- rbeta(100, 0.5, 0.1) 

d2 <- density(dat) 
d2_ks <- bkde(dat) 

par(mfrow=c(2,1)) 
plot(d2) 
plot(d2_ks, typ="l") 

auc(d2$x, d2$y) 
## [1] 1.000187 

integrate.xy(d2$x, d2$y) 
## [1] 1.000188 

auc(d2_ks$x, d2_ks$y) 
## [1] 1 

integrate.xy(d2_ks$x, d2_ks$y) 
## [1] 1 

auc e integrate.xy entrambi utilizzano la regola del trapezio, ma li ho corse sia per dimostrare che e mostrare i risultati di due diverse funzioni.

Il punto è che le densità di fatto sommano a 1, nonostante i valori dell'asse y che portano a credere che non lo siano. Non sono sicuro di cosa stai cercando di risolvere con le tue manipolazioni.