r - 如何使用统计信息来加速 data.frame 上的逐行计算?
问题描述
我有一个包含 10,000 行和 40 列的数据框。我正在尝试对这些行中的每一行应用一个函数。对于每一行,我期望返回一个标量,它是我在这个函数中计算的统计数据的值。以下是我到目前为止所做的;
library(sandwich)
# Creating example data #
nrows=10000
ncols=40
n1=20
n2=20
df=data.frame(t(replicate(nrows, rnorm(ncols, 100, 3))))
cov=data.frame(group=as.factor(rep(c(1,2),c(n1,n2))))
# Function to evaluate on each row of df #
get_est= function(x){
mod = lm(x~cov$group)
vcov = vcovHC(mod)
coef = as.numeric(mod$coefficients[2])
se = sqrt(as.numeric(diag(vcov)[2]))
stats = coef/se
return(stats)
}
# Applying above function to full data #
t1=Sys.time()
estimates=apply(df, 1, function(x) get_est(x))
t2=Sys.time()-t1
# Time taken by apply function
Time difference of 32.10623 secs
有没有办法显着减少在完整数据上实现 get_est() 所需的时间?我需要在单个 df 上加快计算速度的主要原因是因为我还有 1000 个具有相同维度的数据帧,并且我必须将此函数同时应用于每个数据帧的每一行。为了说明,下面是我正在处理的更广泛的情况;
# Creating example data
set.seed(1234)
nrows = 10000
ncols = 40
n1 = 20
n2 = 20
df.list = list()
for(i in 1:1000){
df.list[[i]] = data.frame(t(replicate(nrows, rnorm(ncols, 100, 3))))
}
# Applying get_est() to each row and to each of data frame in df.list #
all.est = foreach(j = 1:length(df.list), .combine = cbind, .packages = 'sandwich') %dopar% {
cov=data.frame(group=as.factor(rep(c(1,2),c(n1,n2))))
est = apply(df.list[[j]], 1, function(x) get_est(x))
return(est)
}
即使在并行化之后,也需要数小时才能完成。我的最终目标是显着缩短获取“all.est”的时间,该“all.est”将包含 10000 行和 1000 列,其中每列都有相应数据集的统计估计值。任何帮助深表感谢!!提前致谢!
解决方案
您的函数get_est
使用了一些“昂贵”的函数,例如lm
,vcovHC
等。如果你想到 OLS 方程,
$$ \hat{\beta} = (X^TX)^{-1}X^Ty, $$
然后你可以看到第一部分$(X^TX)^{-1}X^T$在你的模拟中没有改变,所以设计矩阵是恒定的。为了利用这一点,我在开始模拟之前计算了 $(X^TX)^{-1}X^T$ 。这种方法还需要使用公式手动计算 HC3 标准误差(参见例如此处)
$$ \widehat{\text{Cov}}_{\text{HC3}}(\hat{\beta}) = (X^TX)^{-1}X^T \text{diag} \left[ \ frac{e_i^2}{(1-h_{ii})^2} \right] X(X^TX)^{-1}。$$
在您的模拟迭代中,除了残差之外的所有内容都是恒定的,因此可以预先计算。一旦我实施了这些技巧,我的速度大约提高了 50 倍。
(注意:lm
使用 QR 分解,在这里也可以类似地实现。也许您可以通过并行化代码来提高速度。)
nrows = 10000
ncols = 40
n1 = 20
n2 = 20
df = data.frame(t(replicate(nrows, rnorm(ncols, 100, 3))))
cov = data.frame(group=as.factor(rep(c(1,2),c(n1,n2))))
# old function
get_est_old = function(x){
mod = lm(x~cov$group)
vcov = sandwich::vcovHC(mod)
coef = as.numeric(mod$coefficients[2])
se = sqrt(as.numeric(diag(vcov)[2]))
stats = coef/se
return(stats)
}
# new function
# first construct design matrix
X = matrix(c(rep(1, ncols), rep(0, ncols / 2), rep(1, ncols / 2)), ncol = 2)
# these quantities will be used below
inv = solve(crossprod(X)) %*% t(X)
h = diag(X %*% inv)
get_est_new= function(x){
coef = (inv %*% x)
resid = x - (X %*% coef)
bread = (resid^2 / (1 - h)^2)[,1]
hc3 = inv %*% diag(bread) %*% t(inv)
se = sqrt(hc3[2,2])
stats = coef[2,1]/sqrt(hc3[2,2])
}
# Applying above function to full data #
system.time({
estimates_old = apply(df, 1, function(x) get_est_old(x))
})
#> user system elapsed
#> 7.876 0.042 7.929
system.time({
estimates_new = apply(df, 1, function(x) get_est_new(x))
})
#> user system elapsed
#> 0.141 0.016 0.158
# check
all.equal(estimates_old, estimates_new)
#> [1] TRUE
由reprex 包(v2.0.1)于 2021-09-04 创建
这些帖子可能很有趣:
推荐阅读
- c - C程序不读取键盘输入
- r - 导致此特定错误消息的语法错误是什么?
- apache-kafka - 如果在 Apache Flink 中操作时发生异常,则不提交来自 Apache Kafka 的消息
- javascript - 如何使用WEB API Dot net core实现文件上传?
- javascript - 如何使用jQuery从字符串中获取特定值
- python - 在 groupby 聚合函数中有条件地连接字符串
- spring-boot - Spring Boot:请求参数中的自定义验证
- oracle - 将 SQL 的 first_value 和 partition by 转换为 SAS
- c++ - 我们可以在 C++17 中检测到“微不足道的可重定位性”吗?
- docker - 错误:作业失败(系统故障):无法在 unix:///var/run/docker.sock 连接到 Docker 守护程序。docker 守护进程是否正在运行?在 Windows 10 上