优化使用Expand.grid(和dplyr?)的模拟

Question

我具有此基本函数'laste_up'（比我的真实函数更简单），它从线性回归中提取一些信息。这两个参数是单个（n）的数量和通过个人（RP）重复的数量。

follow_up <- function(n=10,rp=5){
donnees <- data.frame(Id=rep(1:n, rep(rp,n)),X=rnorm(n*rp),Y=rnorm(n*rp,4,2))
sfit <- summary(lm(Y~X, donnees))
output <- c(sfit$R.chisq,sfit$coeff[1],sfit$coeff[2])
return(output)
}
n=c(5,10)
rp=c(2,4,6)
p=expand.grid(n=n,rp=rp)

我想实现大型模拟，并获得以下内容：

对于N和RP的每种组合（均坐在Expand.grid提供的两个第一列中），我想实现大约1,000次迭代，对每次迭代进行评估" colloct_up"功能，然后放入数据帧的其他列返回的三个组件的平均值（即colloct_up"（即R2和Coeffs的平均值）。

因为我的实际功能更为复杂，并且由于N和RP具有更高的维度，因此我想优化我的代码（例如，如果可能的话，请避免使用RBIND或循环）。谢谢您的帮助。

Answer 1

您可以：

set.seed(1)
follow_up_vectorized <- Vectorize(follow_up)
sims <- replicate(1e3, follow_up_vectorized(p$n, p$rp))
res <- apply(sims, c(1, 2), mean)
#             [,1]         [,2]        [,3]        [,4]        [,5]         [,6]
# [1,]  4.00783364  3.991355959 4.011558264 3.983996744  3.99937381  4.009033518
# [2,] -0.03425608 -0.004379941 0.005743333 0.005036114 -0.01332833 -0.007702833

但我不知道您实际代码的性能瓶颈而不会称其为"优化"。

---

编辑
根据CPAK的要求，输出为新列：

cbind(p, t(res))
#    n rp        1            2
# 1  5  2 4.007834 -0.034256082
# 2 10  2 3.991356 -0.004379941
# 3  5  4 4.011558  0.005743333
# 4 10  4 3.983997  0.005036114
# 5  5  6 3.999374 -0.013328326
# 6 10  6 4.009034 -0.007702833

Answer 2

您的数据

n=c(5,10)
rp=c(2,4,6)
p=expand.grid(n=n,rp=rp)

更改功能返回数据框

follow_up_df <- function(n=10,rp=5){
                  donnees <- data.frame(Id=rep(1:n, rep(rp,n)),X=rnorm(n*rp),Y=rnorm(n*rp,4,2))
                  sfit <- summary(lm(Y~X, donnees))
                  output <- c(sfit$R.chisq, sfit$coeff[1], sfit$coeff[2])
                  df <- data.frame(X1=output[1], X2=output[2])
                  return(df)
                }

平淡的解决方案

CP <- function() {
          require(tidyverse)
          totiter <- 1000
          # Copy p 1000 times
          p1 = p[rep(seq_len(nrow(p)), totiter ), ] %>%                     
                 mutate(ID = seq_len(totiter*nrow(p)))                      # unique ID to join
          # Calculate mean of N iterations
          ans <- map_df(1:nrow(p1), ~follow_up(p1$n[.x], p1$rp[.x])) %>%    # follow_up rowwise
                    mutate(ID = seq_len(totiter*nrow(p))) %>%               # unique ID to join
                    left_join(., p1, by="ID") %>%                           # join with p1
                    group_by(n, rp) %>%     
                    summarise(X1 = mean(X1), X2 = mean(X2)) %>%             # mean per n,rp pair
                    ungroup()
      }

输出

set.seed(1)
      n    rp       X1           X2
1     5     2 4.007834 -0.034256082
2     5     4 4.011558  0.005743333
3     5     6 3.999374 -0.013328326
4    10     2 3.991356 -0.004379941
5    10     4 3.983997  0.005036114
6    10     6 4.009034 -0.007702833

其他解决方案

 Aurele <- function() {
              set.seed(1)
              follow_up_vectorized <- Vectorize(follow_up)
              sims <- replicate(1e3, follow_up_vectorized(p$n, p$rp))
              res <- apply(sims, c(1, 2), mean)
          }    

性能

library(microbenchmark)
microbenchmark(CP(), times=5L)   
    expr      min       lq     mean   median       uq      max neval
    CP() 25.02497 25.58269 25.83376 25.92396 26.26672 26.37044     5
Aurele() 21.31826 21.44110 21.73005 21.79842 21.85301 22.23944     5

结论

Aurele's solution is faster!