編碼中的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的y使得x[i]超過既y[j1]和y[j2]j2數，，然後被求和這個數量全部爲i。 這裏有一個快速的方法來做a12（另一個將作爲練習）。你可以翻轉求和的順序首先要注意：

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

即針對每個不同的指數對j1,j2，我們計算出超過兩y[j1]和y[j2]，然後我們在所有這些總結這個量x元素的數量不同的索引對。現在計算j1,j2對的內部總和就像是一個矩陣乘法。事實上，假如我們定義矢量x和y：

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

那麼我們就可以用outer函數產生一個矩陣y_x其[i,j]項爲真當且僅當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，我們通過使用Z1`的`

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

> a12 
[1] 72