biểu đồ ggplot2 có đường cong mật độ tổng cộng tới 1

Question

Vẽ biểu đồ với đường cong mật độ tổng cộng là 1 cho dữ liệu không chuẩn hóa là ridiculously khó. Có rất nhiều câu hỏi về vấn đề này, nhưng không có giải pháp nào làm việc cho dữ liệu của tôi. Có cần phải là một giải pháp đơn giản mà chỉ hoạt động. Tôi không thể tìm thấy câu trả lời với một giải pháp đơn giản hoạt động.

Một số ví dụ:

giải pháp duy nhất làm việc với dữ liệu bình thường chuẩn ggplot2: Overlay histogram with density curve

với dữ liệu rời rạc và không có đường cong mật độ ggplot2 density histogram with width=.5, vline and centered bar positions

có câu trả lời Overlay density and histogram plot with ggplot2 using custom bins

mật độ không tổng hợp để 1 trên dữ liệu của tôi Creating a density histogram in ggplot2?

không tổng hợp tới 1 trên dữ liệu của tôi giải thích ggplot2 density histogram with custom bin edges

dài ở đây với các ví dụ, nhưng mật độ không phải là 1 với dữ liệu của tôi "Density" curve overlay on histogram where vertical axis is frequency (aka count) or relative frequency?

-

Một số mã ví dụ:

#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() 

Vì vậy, rõ ràng mật độ không tổng hợp để 1.

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

Nó đã làm thay đổi một cái gì đó, nhưng không chính xác.

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

Không thay đổi ở đó.

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

Tuy nhiên sai đường cong mật độ, nhưng bây giờ histogram là sai quá.

Để chắc chắn, tôi đã dành 4 giờ để thử tất cả các loại kết hợp ..từ .. và ..sum .. và .. mật độ .., nhưng vì tôi không thể tìm thấy bất kỳ tài liệu nào về cách chúng có nghĩa vụ phải làm việc, đó là phiên tòa bán mù và lỗi.

Vì vậy, tôi đã từ bỏ và tránh sử dụng ggplot2 để tóm tắt dữ liệu.

Vì vậy, đầu tiên chúng ta cần phải nhận được tỷ lệ data.frame đúng, và đó không phải là đơn giản như vậy:

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) 

này cung cấp cho các loại dữ liệu tóm tắt chúng ta muốn:

> 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 

Bây giờ chúng ta chỉ cần vẽ nó. Nên dễ dàng ...

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

Umm, không hoàn toàn những gì tôi muốn. Để chắc chắn, tôi đã thử mà không có stat = "identity" trong geom_density, tại thời điểm đó nó phàn nàn về việc không có một y.

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

Thậm chí kỳ lạ hơn.

Được rồi, có thể chúng ta hãy từ bỏ việc nhận đường cong mật độ từ dữ liệu tóm tắt. Có lẽ chúng ta cần phải kết hợp các phương pháp một chút ...

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

Ok, ít nhất là hình là ngay bây giờ. Bây giờ, chúng ta cần phải bằng cách nào đó quy mô nó xuống.

#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') 

Hình như chúng ta có một người chiến thắng. Ngoại trừ số được mã hóa cứng.

#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 

Vâng, tôi gần như mong đợi nó hoạt động. Bây giờ tôi đã thử thêm một số .. 's ở đây và ở đó, cũng .. số .. và ..sum .., người đầu tiên mà đã cho một kết quả sai, thứ hai mà đã ném một lỗi. Tôi cũng đã thử sử dụng một số nhân (với 1/20), không may mắn.

#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') 

Vì vậy, cuối cùng tôi nhận được con số đúng (tôi nghĩ, tôi hy vọng).

Hãy cho tôi biết cách thực hiện điều này dễ dàng hơn.

PS. Thủ thuật get() dường như không hoạt động trong một hàm. Tôi đã có thể đặt một chức năng làm việc ở đây để sử dụng trong tương lai, nhưng điều đó cũng không dễ dàng như vậy.

Answer 1

Trước tiên, hãy đọc Wickham về mật độ trong R, lưu ý các điểm yếu và các tính năng của từng gói/chức năng.

Tổng mật độ là 1, nhưng điều đó không có nghĩa là đường/điểm đường cong sẽ không vượt quá 1.

Các chương trình sau cả này và thiếu chính xác của (ít nhất) giá trị mặc định của density khi so sánh với, nói, KernSmooth::bkde (sử dụng lô cơ sở cho ngắn gọn đánh máy):

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 

Làm tương tự cho bản phân phối 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 và integrate.xy đều sử dụng quy tắc hình thang nhưng tôi chạy chúng cho cả hai thấy và để hiển thị các kết quả từ hai chức năng khác nhau.

Vấn đề là mật độ thực tế tính tổng là 1, mặc dù các giá trị trục y dẫn bạn tin rằng chúng không có. Tôi không chắc chắn những gì bạn đang cố gắng giải quyết với các thao tác của bạn.