如何根据 R 中的密度进行分组汇总统计

How to make grouped summary statistics based off of densities in R

目标:我想为每个组生成分组百分位数 (hrzn)

我有以下数据

# A tibble: 3,500 x 3
    hrzn parameter density
   <dbl>     <dbl>   <dbl>
 1     1    0.0183 0.00914
 2     1    0.0185 0.00905
 3     1    0.0187 0.00897
 4     1    0.0189 0.00888
 5     1    0.0191 0.00880
 6     1    0.0193 0.00872
 7     1    0.0194 0.00864
 8     1    0.0196 0.00855
 9     1    0.0198 0.00847
10     1    0.0200 0.00839

hrzn是组,parameter是参数的网格space,density是[=14]中值的密度=]列。

我想通过 hrzn 生成 10 到 90 乘以 10 的统计百分位数摘要。我正在努力保持这种计算效率。我知道我可以用密度作为权重对参数进行采样,但我很好奇是否有一种更快的方法可以在不进行采样的情况下从密度生成百分位数。

数据可以通过以下方式获取

df <- readr::read_csv("https://raw.githubusercontent.com/alexhallam/density_data/master/data.csv")

当我从你的 csv 加载数据时,5 个组中的每一个都具有相同的参数和密度值:

df
#># A tibble: 3,500 x 3
#>    hrzn parameter density
#>   <int>     <dbl>   <dbl>
#> 1     1    0.0183 0.00914
#> 2     1    0.0185 0.00905
#> 3     1    0.0187 0.00897
#> 4     1    0.0189 0.00888
#> 5     1    0.0191 0.00880
#> 6     1    0.0193 0.00872
#> 7     1    0.0194 0.00864
#> 8     1    0.0196 0.00855
#> 9     1    0.0198 0.00847
#>10     1    0.0200 0.00839
#># ... with 3,490 more rows

sapply(1:5, function(x) all(df$parameter[df$hrzn == x] == df$parameter[df$hrzn == 1]))
# [1] TRUE TRUE TRUE TRUE TRUE

sapply(1:5, function(x) all(df$density[df$hrzn == x] == df$density[df$hrzn == 1]))
# [1] TRUE TRUE TRUE TRUE TRUE

我不确定这是不是一个错误,但很明显,如果你担心计算,你想在所有组上做的任何事情都可以通过只在一个上做快 5 倍组.

总之,要得到每个hrzn的第10个和第90个百分位数,您只需要查看累积分布函数上哪个参数与0.1和0.9相邻即可。让我们概括为所有组解决这个问题,以防数据出现问题或者您想用不同的数据重复它:

library(dplyr)

df %>% 
  mutate(hrzn = factor(hrzn)) %>%
  group_by(hrzn) %>% 
  summarise(centile_10 = parameter[which(cumsum(density) > .1)[1]],
            centile_90 = parameter[which(cumsum(density) > .9)[1]] )

#># A tibble: 5 x 3
#>  hrzn  centile_10 centile_90
#>  <fct>      <dbl>      <dbl>
#>1 1         0.0204      0.200
#>2 2         0.0204      0.200
#>3 3         0.0204      0.200
#>4 4         0.0204      0.200
#>5 5         0.0204      0.200

当然,由于上述原因,它们都是一样的。

如果您担心计算时间(即使上面只需要几毫秒),并且您不介意不透明的代码,您可以利用 cut [=整个 density 列的 15=] 以 0.1 的步长在 0 和 5 之间,以获得所有第 10 个百分位数,如下所示:

summary <- df[which((diff(as.numeric(cut(cumsum(df$density), seq(0,5,.1))) - 1) != 0)) + 1,]
summary <- summary[-(1:5)*10,]
summary$centile <- rep(1:9*10, 5)
summary
#> # A tibble: 45 x 4
#>     hrzn parameter density centile
#>    <int>     <dbl>   <dbl>   <dbl>
#>  1     1    0.0204 0.00824      10
#>  2     1    0.0233 0.00729      20
#>  3     1    0.0271 0.00634      30
#>  4     1    0.0321 0.00542      40
#>  5     1    0.0392 0.00453      50
#>  6     1    0.0498 0.00366      60
#>  7     1    0.0679 0.00281      70
#>  8     1    0.103  0.00199      80
#>  9     1    0.200  0.00114      90
#> 10     2    0.0204 0.00824      10
#> # ... with 35 more rows

也许我误会你了,你实际上是在5维参数space上工作,想知道5d密度的第10个和第90个百分位数的参数值。在这种情况下,您可以利用所有组都相同的事实来计算 5 维密度的第 10 和第 90 个百分位数,只需取这两个百分位数的第 5 个根即可:

df %>% 
  mutate(hrzn = factor(hrzn)) %>%
  group_by(hrzn) %>% 
  summarise(centile_10 = parameter[which(cumsum(density) > .1^.2)[1]],
            centile_90 = parameter[which(cumsum(density) > .9^.2)[1]] )

#> # A tibble: 5 x 3
#>   hrzn  centile_10 centile_90
#>   <fct>      <dbl>      <dbl>
#> 1 1         0.0545      0.664
#> 2 2         0.0545      0.664
#> 3 3         0.0545      0.664
#> 4 4         0.0545      0.664
#> 5 5         0.0545      0.664