R - Vectorize嵌套循环

Question

对不起，另一个“矢量化for循环”的问题，但我一直没能弄清楚如何做到这一点。我试图编写的功能很简单：

对于enroll.in中的每一行，首先使用hasMedClaims逻辑模型输出作为响应的概率。

生成随机数并使用它来确定是否应建模响应。

如果是，则为响应建模。如果不是，只需输入0.对于每个enroll.in nsim时间行重复一次。

simMedClaims.loop<-function(hasMedClaims.in, MedClaims.in, enroll.in, nsim = 100){ 

    set.seed(100) 
    #dataframe to hold results 
    results<-matrix(0, ncol = nsim, nrow = nrow(enroll.in)) 
    results<-data.frame(results) 

    hasclaims<-predict(hasMedClaims.in, newdata = enroll.in, type = "response") 
    means<-predict(MedClaims.in, newdata = enroll.in, type="response") 
    for(ii in 1:nrow(enroll.in)) 
    { 
    for(jj in 1:nsim){ 
     unif.rand<-runif(1) 
     results[ii,jj]<-ifelse(unif.rand < hasclaims[ii], exp(rnorm(1,mean = means[ii], sd = sqrt(MedClaims.in$sig2))), 0) 
    } 

    } 

    return(results) 
} 

set.seed(100) 
dummy<-data.frame(hasresponse = rbinom(100000, 1, .5), response = rnorm(100000, mean = 5, sd = 1), x1 = runif(100000, 0, 60), x2 = as.factor(rbinom(100000, 1, .5)+1)) 
dummy$response<-dummy$hasresponse*dummy$response 
hasresponse_gam<-mgcv::gam(hasresponse ~ s(x1,bs="ps", by=x2)+x2, data=dummy, family = binomial(link="logit"), method="REML") 
response<-mgcv::gam(response ~ s(x1,bs="ps", by=x2)+x2, data=dummy[dummy$hasresponse==1,]) 
dummyEnroll<-data.frame(x1 = runif(10, 20, 50), x2 = as.factor(rbinom(10, 1, .5)+1)) 
system.time(result<-simMedClaims.loop(hasresponse_gam, response, dummyEnroll, 1000)) 

user system elapsed 
38.66 0.00 39.35 

我已经尝试了很多不同的想法，但我得到了每个人不同的问题。

hasMedClaims.in和MedClaims.in都是使用mgcv gam函数拟合的GAM。

澄清我为什么问这个问题：正如输出显示，每个主题需要几秒钟才能运行1000次模拟。我将在具有数万个主题的数据集上使用它，并且我想运行至少50,000次模拟。我目前的代码有效，但速度太慢。我的目标是优化我的功能以更快地运行。

尝试在@芭菲的FUNC2

simMedClaims2<-function(hasMedClaims.in, MedClaims.in, enroll.in, nsim = 100){ 
    set.seed(100) 
    hasclaims<-predict(hasMedClaims.in, newdata = enroll.in, type = "response") 
    means<-predict(MedClaims.in, newdata = enroll.in, type="response") 
    results<-data.frame(t(vapply(seq(nrow(enroll.in)), function(ii, jj){ 
    ifelse(runif(jj) < hasclaims[ii],1,0)*exp(rnorm(nsim,mean = means[ii], sd = sqrt(MedClaims.in$sig2))) 
    },numeric(nsim),seq(nsim)))) 
    return(results) 
} 

结果看起来合理的，虽然我还没有完全核实他们没有。我还编辑了我的原始循环函数来计算循环外部的方法。快得多

> system.time(result<-simMedClaims.loop(hasresponse_gam, response, dummyEnroll, 100)) 
    user system elapsed 
    0.06 0.00 0.13 
> system.time(result2<-simMedClaims2(hasresponse_gam, response, dummyEnroll, 100)) 
    user system elapsed 
    0.02 0.00 0.02 

但是，运行all.equal(result, result2)表明输出不等效。我无法弄清楚为什么。

Answer 1

考虑传递两个矢量参数在sapply或vapply避免嵌套for循环和需要初始化结果数据帧。当然，它仍然是值得商榷的，如果apply family is truly vectorized：

simMedClaims.loop <- function(hasMedClaims.in, MedClaims.in, enroll.in, nsim = 100){ 

    hasclaims <- predict(hasMedClaims.in, newdata = enroll.in, type = "response") 

    results <- data.frame(t(vapply(seq(nrow(enroll.in)), function(ii,jj) { 
             unif.rand <- runif(jj) 
             ifelse(unif.rand < hasclaims[ii], ..., 0) 
            numeric(nsim), seq(nsim))))  
} 

另外，考虑一个expand.grid()方法中，在结束争吵为多列的需要的格式。虽然没有数据争议，但这将是矢量化的（不使用R循环，但可能是C循环）。

simMedClaims.loop <- function(hasMedClaims.in, MedClaims.in, enroll.in, nsim = 100){ 

    hasclaims <- predict(hasMedClaims.in, newdata = enroll.in, type = "response") 

    # LONG FORMAT 
    df <- expand.grid(1:nrow(enroll.in), 1:nsim) 
    df$unif.rand <- runif(nrow(df)) 
    df$val <- ifelse(df$unif.rand < hasclaims[ii], ..., 0) 

    # WIDE FORMAT 
    results <- data.frame(t(sapply(seq(1, nrow(df), by=nsim), function(i) 
           df$random_num[i:(i+(nsim-1))]))) 

} 

---

以上方法进行了测试用随机数据并返回相同的结果嵌套for循环（不包括OP的predict或ifelse由于没有reproducible example）：

数据

enroll.in <- sapply(1:5, function(i) rnorm(15)) 
nsim <- 100 

方法

func1 <- function() {  
    set.seed(98) 
    results1<-matrix(0, ncol = nsim, nrow = nrow(enroll.in)) 
    results1<-data.frame(results1) 

    for(ii in 1:nrow(enroll.in)) 
    { 
    for(jj in 1:nsim){ 

    results1[ii,jj] <- runif(1) 
    } 
    } 
    return(results1) 
} 

func2 <- function() { 
    set.seed(98) 
    results2 <- data.frame(t(vapply(seq(nrow(enroll.in)), function(ii,jj) 
             runif(jj), 
            numeric(nsim), seq(nsim)))) 
} 

func3 <- function() { 
    set.seed(98) 
    df <- expand.grid(1:nrow(enroll.in), 1:nsim) 
    df$random_num <- runif(nrow(df)) 

    results3 <- data.frame(t(sapply(seq(1, nrow(df), by=nsim), function(i) 
            df$random_num[i:(i+(nsim-1))]))) 
} 

成果

all.equal(func1(), func2()) 
# [1] TRUE 
all.equal(func2(), func3()) 
# [1] TRUE 

和基准，表明至少在小数据，处理不是方法之间的任何好得多。注：大纳秒处理是由于函数'set.seed()为了比较随机生成的数据。因此，古老的谚语认为：有什么错for循环：

library(microbenchmark) 

microbenchmark(func1) 
# Unit: nanoseconds 
# expr min lq mean median uq max neval 
# func1 30 32 37.07  32 33 461 100 

microbenchmark(func2) 
# Unit: nanoseconds 
# expr min lq mean median uq max neval 
# func2 29 31 39.41  32 33 729 100 

microbenchmark(func3) 
# Unit: nanoseconds 
# expr min lq mean median uq max neval 
# func3 30 31 35.6  32 33 370 100