我无法运行这个双循环来将计算值正确存储到矩阵中(如下所述).我选择使用双For循环而不是apply()或mean()的原因是我想获得两列的唯一组合并消除冗余(如下所述).请参阅下面的示例:
A<-c(1,2,3,4,5)
B<-c(2,3,4,5,6)
Q1<-data.frame(cbind(A,B))
mean<-matrix(nrow=5, ncol = 5)
for(i in 1: length(Q1$A)){
for(j in 2: length(Q1$B)){
mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
}
}
在这里,我尝试通过整个B向量运行整个A向量,同时消除冗余,使得A [1]具有来自B [2]的四个值,并且A [2]具有来自B [3]的三个值.但是,这是我的结果.
[,1] [,2] [,3] [,4] [,5]
[1,] NA 2.0 2.5 3.0 3.5
[2,] NA 2.5 3.0 3.5 4.0
[3,] NA 3.0 3.5 4.0 4.5
[4,] NA 3.5 4.0 4.5 5.0
[5,] NA 4.0 4.5 5.0 5.5
虽然第一列是我的预期,但我有我不想要的价值观.我想要的是以下矩阵输出:
[,1] [,2] [,3] [,4] [,5]
[1,] NA 2.0 2.5 3.0 3.5
[2,] NA NA 3.0 3.5 4.0
[3,] NA NA NA 4.0 4.5
[4,] NA NA NA NA 5.0
[5,] NA NA NA NA NA
有什么建议？
[原始解决方案(有关更快的解决方案,请参阅更新2)]
f.m <- function(Q1) {
z <- matrix(nrow=nrow(Q1),ncol=nrow(Q1))
b <- row(z) < col(z)
z[b] <- (Q1$A[col(z)[b]] + Q1$B[row(z)[b]])/2
z
}
[样本输出]
f.m(Q1)
# [,1] [,2] [,3] [,4] [,5]
# [1,] NA 2 2.5 3.0 3.5
# [2,] NA NA 3.0 3.5 4.0
# [3,] NA NA NA 4.0 4.5
# [4,] NA NA NA NA 5.0
# [5,] NA NA NA NA NA
[基准设置]
f0 <- function(Q1) {
mean<-matrix(nrow=nrow(Q1), ncol = nrow(Q1))
for(i in 1: length(Q1$A)){
for(j in 2: length(Q1$B)){
mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
}
}
mean
}
f1 <- function(Q1) {
mean<-matrix(nrow=nrow(Q1), ncol = nrow(Q1))
for(i in 2: length(Q1$A)){
for(j in i: length(Q1$B)){
mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
}
}
mean
}
# Note that f0() and f1() don't return the desired result for the sample output
f2 <- function(Q1) {
mean<-outer(1: length(Q1$A),
1: length(Q1$B),
Vectorize(function(i,j){
if(i >= j)
return(NA)
else
return(sum(Q1$A[i]+Q1$B[j])/2)
}))
mean
}
library(rbenchmark)
[基准测试结果]
A <- B <- 1:100
Q1<-data.frame(A,B)
benchmark(f0(Q1), f1(Q1), f2(Q1), f.m(Q1), replications = 10)
test replications elapsed relative user.self sys.self user.child sys.child
4 f.m(Q1) 10 0.011 1.000 0.012 0.000 0 0
1 f0(Q1) 10 3.018 274.364 3.007 0.008 0 0
2 f1(Q1) 10 1.477 134.273 1.474 0.003 0 0
3 f2(Q1) 10 1.777 161.545 1.774 0.002 0 0
[更新1]
通过直接计算整个矩阵可以节省另一个运行时间顺序,这避免了与代价高昂的(与求和相比)子集的混乱:
f.m2 <- function(Q1) outer(Q1$A,Q1$B,'+')*0.5
基准测试的另一部分:
A <- B <- 1:1000
Q1<-data.frame(A,B)
#benchmark(f0(Q1), f1(Q1), f2(Q1), f.m(Q1), replications = 10)
benchmark(f.m(Q1), f.m2(Q1), replications = 10)
test replications elapsed relative user.self sys.self user.child sys.child
1 f.m(Q1) 10 1.839 10.274 1.746 0.093 0 0
2 f.m2(Q1) 10 0.179 1.000 0.144 0.035 0 0
[更新2]
1)如David Arenburg所述,函数f.m2()不能精确地产生预期的输出,因为输出的左下三角和主对角线应该用NA填充.f.m2()可以修复该功能,以牺牲性能为代价来产生正确的答案(参见下面的基准测试).
# Suggested by David Arenburg
f.m2.1 <- function(Q1) {
Res <- outer(Q1$A,Q1$B,'+')*0.5;
Res[lower.tri(Res, diag = TRUE)] <- NA;
Res
}
2)这是David Arenburg建议的另一种方法,它利用CJ了data.table包中的功能:
library(data.table)
f.DA <- function(Q1){
Res <- matrix(rowMeans(CJ(Q1$A, Q1$B)), ncol = nrow(Q1))
Res[lower.tri(Res, diag = TRUE)] <- NA
Res
}
3)这是一种Rcpp基于方法的方法:
library(Rcpp)
cppFunction('NumericMatrix fC(NumericVector A, NumericVector B) {
int n = A.size();
NumericMatrix out(n,n);
std::fill( out.begin(), out.end(), NumericVector::get_na() ) ;
for (int i = 0; i < n; i++) {
for (int j = i+1; j < n; j++) {
out(i,j) = 0.5*(A[i] + B[j]);
}
}
return out;
}')
4)另一项基准研究:
A <- B <- 1:3000
Q1<-data.frame(A,B)
benchmark(f.m2(Q1), f.m2.1(Q1), f.DA(Q1), fC(Q1$A, Q1$B), replications = 10)
test replications elapsed relative user.self sys.self user.child sys.child
3 f.DA(Q1) 10 7.442 11.556 6.200 1.209 0 0
2 f.m2.1(Q1) 10 5.111 7.936 4.404 0.661 0 0
1 f.m2(Q1) 10 1.007 1.564 0.733 0.263 0 0
4 fC(Q1$A, Q1$B) 10 0.644 1.000 0.525 0.116 0 0