Double For Loop для вычисления средних значений и сохранения их в матрице

У меня возникли проблемы с запуском этого двойного цикла for, чтобы правильно сохранить вычисленные значения в матрице (упомянутой ниже). Причина, по которой я решил использовать двойной цикл For Loop, а не 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
  }
}

Здесь я попытался пропустить весь вектор A через весь вектор B, исключив избыточность, так что 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

Какие-либо предложения?


person Provisional.Modulation    schedule 22.02.2015    source источник
comment
Почему вас интересует только половина матрицы? Например, в случае colA <- 1:3 и colB <- 13:11 выходная матрица становится асимметричной (например, A[1] + B[3] != A[3] + B[1]), и вы потеряете информацию, глядя только на половину матрицы.   -  person Marat Talipov    schedule 22.02.2015
comment
@MaratTalipov Меня интересует половина, потому что я хочу взять эти значения и поместить их в столбец, чтобы я мог сравнить их с другими в ggplot. Если есть избыточность, то она будет отражать результат графика.   -  person Provisional.Modulation    schedule 22.02.2015


Ответы (4)


[Исходное решение (более быстрые решения см. в Обновлении 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
comment
Спасибо за помощь, Марат! - person Provisional.Modulation; 22.02.2015
comment
Фактически f.m2 возвращает всю матрицу и, следовательно, не дает желаемого результата. Вам, вероятно, следует изменить его на f.m2 <- function(Q1) { Res <- outer(Q1$A,Q1$B,'+')*0.5; Res[lower.tri(Res, diag = TRUE)] <- NA; Res }, чтобы выполнить требования, но он все равно будет быстрее, чем f.m, хотя - person David Arenburg; 22.02.2015
comment
Вы также можете добавить 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 }, который также побьет f.m, но, вероятно, будет медленнее, чем f.m2 - person David Arenburg; 22.02.2015
comment
Спасибо @DavidArenburg, я добавил ваши предложения, а также решение на основе Rcpp - person Marat Talipov; 23.02.2015

Второй цикл for должен быть:

 for(j in (i+1):length(Q1$B))
person klash    schedule 22.02.2015

вы хотите использовать ключевое слово next, чтобы пропустить ненужные операции, например:

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)){
    if(i >= j)
        next
    mean[i,j]<-sum(Q1$A[i]+Q1$B[j])/2
}

или вы можете сделать итерацию внутреннего цикла for зависимой от значения внешнего индекса, как в:

mean<-matrix(nrow=5, ncol = 5)
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
    }
}

или вы можете использовать outer(), как в:

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)
            }))
person Jthorpe    schedule 22.02.2015

Не совсем двойной цикл For Loop, но вы можете просто использовать функцию outer для вычисления средних значений.

outer(Q1$Col1, Q1$Col2, "+")/2
person Community    schedule 06.03.2015