кодирование двойной суммы в R

Question

Я хотел бы вычислить их в количествах

a12=sum_(i from 1 to m)sum_(j1<j2)(I(X[i]>Y[j1] and X[i]>Y[j2]))

a13=sum_(j from 1 to n)sum_(i1<i2)(I(X[i1]>Y[j] and X[i2]>Y[j]))

, где I - функция индикатора.

Итак, я пришел с этим кодом R

a12=0; a13=0

for (l in 1:(length(Z1)-1)){

 for (m in  1:(length(Z2)-1)){

 a12<-a12+(Z1[l]<Z2[m])*(Z1[l+1]<Z2[m])*1

 a13<-a13+(Z1[l]<Z2[m])*(Z1[l]<Z2[m+1])*1

        } # closing m

          } # closing l

    a12=a12+sum((Z1[-length(Z1)]<Z2[length(Z2)])*(Z1[-1]<Z2[length(Z2)])*1)

    a13=a13+sum((Z1[length(Z1)]<Z2[-length(Z2)])*(Z1[length(Z1)]<Z2[-1])*1)


a12;
a13

К сожалению, это не только очень медленно, но я не получаю того, что должен получить.

Не могли бы вы мне помочь, пожалуйста!

Спасибо,

Роланд

Answer 1

Я предполагаю (для a12), что вы хотите сделать следующее. У вас есть два вектора x (длины m) и y, и для каждого элемента x[i] из x вы рассчитываете число различных индексные пары j1, j2 из y такие, что x[i] превосходит y[j1] и y[j2], и затем вы суммируете это количество по всем i. Вот быстрый способ сделать a12 (другой будет оставлен в качестве упражнения). Сначала обратите внимание, что вы можете перевернуть порядок суммирования:

a12 = Sum_(j1 < j2) Sum_(i=1:m) I( X[i] > Y[j1] & X[i] > Y[j2] ),

т.е. для каждой отдельной пары индексов j1,j2 мы вычисляем количество элементов x, которые превышают как y[j1], так и y[j2], а затем суммируем эту величину по всем этим отдельным парам индексов. Теперь вычисление внутренней суммы для пар j1,j2 похоже на умножение матриц. Действительно, предположим, что мы определяем векторы x и y:

set.seed(1)
x <- sample(1:5,5,T)
y <- sample(1:5,10,T)

тогда мы можем использовать функцию outer для получения матрицы y_x, чья запись [i,j] равна TRUE, если и только если y[i] < x[j]:

y_x <- outer(y,x,FUN = '<')

Теперь мы получаем внутренние суммы, выполняя

z <- y_x %*% t(y_x)

, где z[i,j] - количество элементов x, которое превышает как y[i], так и y[j]. Поскольку мы хотим суммировать z[i,j] только для различных i < j, мы получаем конечный результат, взяв сумму нижнего треугольника z, используя

a12 <- sum( z[lower.tri( z )])

> a12
[1] 72