如何将频率分布转换为R中的概率分布

Posted 2023-02-16

【中文标题】如何将频率分布转换为R中的概率分布

【英文标题】：How to convert frequency distribution to probability distribution in R

【发布时间】：2015-03-31 19:35:09

【问题描述】：

我有一个包含 n 行观察值的矩阵。观测值是特征的频率分布。我想将频率分布转换为每行之和为 1 的概率分布。因此，矩阵中的每个元素都应除以该元素的行之和。

我编写了以下 R 函数来完成这项工作，但它对于大型矩阵非常慢：

prob_dist <- function(x) 

    row_prob_dist <- function(row) 
       return (t(lapply(row, function(x,y=sum(row)) x/y)))
       

    for (i in 1:nrow(x)) 
       if (i==1) p_dist <- row_prob_dist(x[i,])
       else p_dist <- rbind(p_dist, row_prob_dist(x[i,]))
       
    return(p_dist)


B = matrix(c(2, 4, 3, 1, 5, 7), nrow=3, ncol=2)
B
     [,1] [,2]
[1,]    2    1
[2,]    4    5
[3,]    3    7

prob_dist(B)
     [,1]      [,2]    
[1,] 0.6666667 0.3333333
[2,] 0.4444444 0.5555556
[3,] 0.3       0.7     

您能否建议 R 函数来完成这项工作和/或告诉我如何优化我的函数以更快地执行？

【问题讨论】：

t(apply(B, 1, prop.table))?

一般要点：由于您将第一行设为特殊情况，因此请在循环外计算它并执行 for( in 2:nrow(x)) 并删除循环内的 if/else。接下来，由于您事先知道输出矩阵的维度，因此创建一个空的 p_dist<-matrix(NA,nrow=nrow(x),ncol=ncol(x)) 。所有rbind 都在浪费时间。

@DavidArenburg 您可能要提到prop.table 只是sweep 的快捷方式

【参考方案1】：

不用apply，一行一行的矢量化解：

t(t(B) / rep(rowSums(B), each=ncol(B)))
          [,1]      [,2]
[1,] 0.6666667 0.3333333
[2,] 0.4444444 0.5555556
[3,] 0.3000000 0.7000000

或者：

diag(1/rowSums(B)) %*% B

【讨论】：

非常好！我打算挖一个非循环，非应用，但你的更好。

太棒了！第一个比 @DavidArenburg 提出的 apply 版本快 3 倍。第二个对于大矩阵非常慢。

歌利亚有时很有希望获胜 ;)

【参考方案2】：

我不确定您的函数是否有任何价值，因为您可以使用hist 或density 函数来完成相同的结果。此外，apply 的使用将如前所述。但它可以作为一个合理的编程示例。

您的代码中有几个效率低下的地方。

您使用 for 循环而不是矢量化代码。这是非常昂贵的。您应该使用上述 cmets 中提到的 apply。

您正在使用rbind，而不是为您的输出预先分配空间。这也非常昂贵。

out <- matrix(NA, nrow= n, ncol= ncol(B))
for (i in 1:nrow(B)) 
  out[i,] <- row_prob_dist(B[i,])

【讨论】：

Alex，在这种情况下你会如何使用 hist 或 density？

【参考方案3】：

这是一个尝试，但在数据框而不是矩阵上：

df <- data.frame(replicate(100,sample(1:10, 10e4, rep=TRUE)))

我尝试了dplyr 方法：

library(dplyr)
df %>% mutate(rs = rowSums(.)) %>% mutate_each(funs(. / rs), -rs) %>% select(-rs)

结果如下：

library(microbenchmark) 
mbm = microbenchmark(
dplyr = df %>% mutate(rs = rowSums(.)) %>% mutate_each(funs(. / rs), -rs) %>% select(-rs),
t = t(t(df) / rep(rowSums(df), each=ncol(df))),
apply = t(apply(df, 1, prop.table)),
times = 100
)

#> mbm
#Unit: milliseconds
#  expr       min        lq      mean    median        uq       max neval
# dplyr  123.1894  124.1664  137.7076  127.3376  131.1523  445.8857   100
#     t  384.6002  390.2353  415.6141  394.8121  408.6669  787.2694   100
# apply 1425.0576 1520.7925 1646.0082 1599.1109 1734.3689 2196.5003   100

编辑：@David benchmark 更符合 OP，所以如果你要使用矩阵，我建议你考虑他的方法。

【讨论】：

史蒂文，以前从未遇到过带有 %>% 的符号，谷歌搜索也没有发现任何引用。你能指出一些参考文献吗？

@AndresKull - %>% 是管道运算符（来自包 magrittr）。你可以在这里阅读：cran.r-project.org/web/packages/magrittr/vignettes/…

想发布你用来生成那张出色图表的代码吗？

@CarlWitthoft ggplot2 中有一个微基准对象的自动绘图方法。要重现上图，您可以简单地执行ggplot2::autoplot(mbm)

【参考方案4】：

其实我想了一下，最好的 vecotization 就是简单

B/rowSums(B)
#           [,1]      [,2]
# [1,] 0.6666667 0.3333333
# [2,] 0.4444444 0.5555556
# [3,] 0.3000000 0.7000000

实际上，@Stevens 基准测试具有误导性，因为 OP 有一个矩阵，而 Steven 基准测试是在一个数据框上。

这是一个带有矩阵的基准。所以对于矩阵，两个向量化的解决方案都会比dplyr 更好，后者不适用于矩阵

set.seed(123)
m <- matrix(sample(1e6), ncol = 100)

library(dplyr)
library(microbenchmark) 

Res <- microbenchmark(
  dplyr = as.data.frame(m) %>% mutate(rs = rowSums(.)) %>% mutate_each(funs(. / rs), -rs) %>% select(-rs),
  t = t(t(m) / rep(rowSums(m), each=ncol(m))),
  apply = t(apply(m, 1, prop.table)),
  DA = m/rowSums(m),
  times = 100
)