Как векторизовать эту математическую формулу?

Question

Я оценки this formula (от Székely and Rizzo, 2013, страниц 1263), в которой a и b симметричны, nn матрицы матрицы с размерностью.

Как я могу векторизовать секцию с вложенным циклом for? Я чувствую, что есть, вероятно, умные манипуляции с матрицами, которые я мог бы использовать здесь, возможно, вместе с функцией sweep.

Вот что я сейчас:

u_squared <- function(a, b, n) { 
    a. <- .colSums(a, n, n) 
    b. <- .colSums(b, n, n) 
    a.. <- sum(a.) 
    b.. <- sum(b.) 
    u1 <- 0 
    u2 <- 0 
    u3 <- 0 
    for (i in 1:n) { 
    for (j in 1:n) { 
     u1 <- u1 + a[i, j] * b[i, j] 
     u2 <- u2 + a[i, j] * (b.. - 2*b.[j] - 2*b.[i] + 2*b[i, j]) 
     u3 <- u3 + a[i, j] * (b.[i] - b[i, j]) 
    } 
    } 
    u1 <- u1/(n * (n-1)) 
    u2 <- u2/(n * (n-1) * (n-2) * (n-3)) 
    u3 <- u3/(n * (n-1) * (n-2)) 
    return (u1 + u2 - 2 * u3) 
} 

Например, я знаю, что u1 может быть вычислена просто как u1 <- a * b, но я хотел, чтобы размотать всю формулу, чтобы продемонстрировать основную математику. Я ищу аналогичные векторизации для u2 и u3.

Answer 1

Это может быть легко векторизованы

u_squared1 <- function(a, b, n){ 
    b. <- matrix(.colSums(b, n, n), n, n) 
    B1 <- b.. - 2 * b. - 2 * t(b.) + 2*b 
    B2 <- b. - b 
    u1 <- sum(a*b)/(n * (n-1)) 
    u2 <- sum(a*B1)/(n * (n-1) * (n-2) * (n-3)) 
    u3 <- sum(a*B2)/(n * (n-1) * (n-2)) 
    return (u1 + u2 - 2 * u3) 
}