[Исходное решение (более быстрые решения см. в Обновлении 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) Как заметил Дэвид Аренбург, функция 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) Вот еще один подход, предложенный Дэвидом Аренбургом, который использует функцию 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
person
Marat Talipov
schedule
22.02.2015
colA <- 1:3
иcolB <- 13:11
выходная матрица становится асимметричной (например,A[1] + B[3] != A[3] + B[1]
), и вы потеряете информацию, глядя только на половину матрицы. - person Marat Talipov   schedule 22.02.2015