R: «векторизация» тройного цикла

Question

Я написал фрагмент кода в R, который вычисляет двойную сумму так называемой статистики рангов.

Мне нужно повторить вычисление Q минимум 1000 раз, но с 3 петлями внутри, для этого требуется довольно много времени, чтобы сделать это только один раз.

Вот мой код:

#u, a - real numbers 
l <- function(u, a) { 
    -sqrt((1-a)/a)*I(u>=0 & u<a) + sqrt(a/(1-a))*I(u>=a & u<=1) 
} 


# r,s - real number, R,S - vectors of real numbers (equal lengths) 
L<-function(r, s, R, S) { 
    n<-length(R) 
    x<-0 
    for (i in 1:n) { 
    x<-x+l(R[i]/(n+1),r) * l(S[i]/(n+1),s) 
    } 
    1/sqrt(n)*x 
} 

# r, s, X, Y - vectors of real numbers; X and Y must be equally long 
Q<-function(r,s,X,Y) { 
    n<-length(X) 
    R<-rank(X) 
    S<-rank(Y) 
    q<-0 
    for (j in 1:length(r)) { 
    for (k in 1:length(s)) { 
     q<-q+L(r[j],s[k],R,S)^2 

    } 
    } 
    q 
} 

я пытался превратить свои функции, используя sapply и применять, но тогда первая функция не удалась, поскольку размеры г и s не могут быть равны (и не должны длины r, s равно длине X (или Y)).

Есть ли способ создать функцию L, которая принимает 4 вектора и создает матрицу, чтобы я избавился от циклов?

Заранее благодарен!

// Edit:

Я написал альтернативную функцию с помощью mapply:

Q1<-function(r,s,X,Y) { 
    n<-length(X) 
    R<-rank(X) 
    S<-rank(Y) 
    rs <- expand.grid(r,s) 
    q<-do.call(mapply, c(function(r,s) L(r,s,R=R,S=S)^2, unname(rs))) 
    sum(q) 
} 

, но это, кажется, еще медленнее. (.)

Answer 1

Если вы хотите создать все значения L для различных значений r и s, метод петли меньше может быть:

rs <- expand.grid(r=r,s=s); rm(r); rm(s) 
    #edit 
    rs$qrs <- with(rs, L(r, s, R, S)^2) 
    q <- sum(rs$qrs) 

Я не уверен, что это будет быстрее , Существует распространенное, но ошибочное представление о том, что петли в R неэффективны. Большая часть эффективности достигается за счет упрощения внутренних функций.

> set.seed(123) 
> r <- runif(4) 
> s <- runif(3) 
> rs <- expand.grid(r=r,s=s) 
> rs 
      r   s 
1 0.2875775 0.9404673 
2 0.7883051 0.9404673 
3 0.4089769 0.9404673 
4 0.8830174 0.9404673 
5 0.2875775 0.0455565 
6 0.7883051 0.0455565 
7 0.4089769 0.0455565 
8 0.8830174 0.0455565 
9 0.2875775 0.5281055 
10 0.7883051 0.5281055 
11 0.4089769 0.5281055 
12 0.8830174 0.5281055 
> rs$qrs <- with(rs, L(r, s, 1:10, 1:10)^2) 
> q <- sum(rs$qrs) 
> q 
[1] 14.39009 
> rs 
      r   s   qrs 
1 0.2875775 0.9404673 0.0004767998 
2 0.7883051 0.9404673 0.0003911883 
3 0.4089769 0.9404673 6.6571168565 
4 0.8830174 0.9404673 0.0017673788 
5 0.2875775 0.0455565 0.0004767998 
6 0.7883051 0.0455565 0.0003911883 
7 0.4089769 0.0455565 6.6571168565 
8 0.8830174 0.0455565 0.0017673788 
9 0.2875775 0.5281055 0.0004767998 
10 0.7883051 0.5281055 0.0003911883 
11 0.4089769 0.5281055 6.6571168565 
12 0.8830174 0.5281055 0.0017673788