Эффективное вычисление двумерного эмпирического cdf в R/Fortran

Question

Учитывая n*2 матрицы данных X, я хотел бы рассчитать двумерный эмпирический cdf для каждого наблюдения, то есть для каждого i в 1:n, вернуть процент наблюдений с 1-м элементом, не превышающим X[i,1] и 2-й элемент не больше X[i,2].

Из-за вложенного поиска он становится ужасно медленным для n ~ 100k, даже после его переноса на Фортран. Кто-нибудь знает, есть ли лучший способ обработки размеров выборки, как этот?

Изменить: я считаю, что эта проблема похожа (с точки зрения сложности) на поиск тау Кендалла, который имеет порядок O(n^2). В этом случае Найт (1966) имеет алгоритм, чтобы уменьшить его до O(n log(n)). Просто интересно, есть ли какой-нибудь алгоритм O(n*log(n)) для поиска двумерного ecdf уже там.

Редактировать 2: Это код, который у меня есть на Фортране, в соответствии с просьбой. Это вызывается в R обычным способом, поэтому код R здесь опущен. Код предназначен для произвольных измерений, но для конкретной вещи, которую я делаю, двумерный достаточно.

! Calculates multivariate empirical cdf for each point
! n: number of observations
! d: dimension (>=2)
! umat: data matrix
! outvec: vector of ecdf

subroutine mecdf(n,d,umat,outvec)
    implicit none

    integer :: n, d, i, j, k, tempsum
    double precision, dimension(n) :: outvec
    double precision, dimension(n,d) :: umat
    logical :: flag

    do i = 1,n
        tempsum = 0
        do j = 1,n
            flag = .true.
            do k = 1,d
                if (umat(i,k) < umat(j,k)) then
                    flag = .false.
                    exit
                end if
            end do
            if (flag) then
                tempsum = tempsum + 1
            end if
        end do
        outvec(i) = real(tempsum)/n
    end do
    return
end subroutine

Answer 1

Я думаю, что моя первая попытка была не на самом деле в формате ecdf, хотя она отображала точки на интервале [0,1]. Пример, матрица 25 x 2, сгенерированная с помощью:

#M <- matrix(runif(100), ncol=2)
M <- 
structure(c(0.0468267474789172, 0.296053855214268, 0.205678076483309, 
0.467400068417192, 0.968577065737918, 0.435642971657217, 0.929023026255891, 
0.038406387437135, 0.304360694251955, 0.964778139721602, 0.534192910650745, 
0.741682186257094, 0.0848641532938927, 0.405901980120689, 0.957696850644425, 
0.384813814423978, 0.639882878866047, 0.231505588628352, 0.271994129288942, 
0.786155494628474, 0.349499785574153, 0.279077709652483, 0.206662984099239, 
0.777465222170576, 0.705439242534339, 0.643429880728945, 0.887209519045427, 
0.0794123203959316, 0.849177583120763, 0.704594585578889, 0.736909110797569, 
0.503158083418384, 0.49449566937983, 0.408533290959895, 0.236613316927105, 
0.297427259152755, 0.0677345870062709, 0.623845702270046, 0.139933609170839, 
0.740499466424808, 0.628097783308476, 0.678438259987161, 0.186680511338636, 
0.339367639739066, 0.373212536331266, 0.976724133593962, 0.94558056560345, 
0.610417427960783, 0.887977657606825, 0.663434249348938, 0.447939050383866, 
0.755168803501874, 0.478974275058135, 0.737040047068149, 0.429466919740662, 
0.0021107573993504, 0.697435079608113, 0.444197302218527, 0.108997165458277, 
0.856855363817886, 0.891898229718208, 0.93553287582472, 0.991948011796921, 
0.630414301762357, 0.0604106825776398, 0.908968194155023, 0.0398679254576564, 
0.251426834380254, 0.235532913124189, 0.392070295521989, 0.530511683085933, 
0.319339724024758, 0.534880011575297, 0.92030712752603, 0.138276003766805, 
0.213625695323572, 0.407931711757556, 0.605797187192366, 0.424798395251855, 
0.471233424032107, 0.0105366336647421, 0.625802840106189, 0.524665891425684, 
0.0375960320234299, 0.54812005511485, 0.0105806747451425, 0.438266788609326, 
0.791981092421338, 0.363821814302355, 0.157931488472968, 0.47945317090489, 
0.906797411618754, 0.762243523262441, 0.258681379957125, 0.308056800393388, 
0.91944490163587, 0.412255838746205, 0.347220918396488, 0.68236422073096, 
0.559149842709303), .Dim = c(50L, 2L))

Поэтому задача состоит в том, чтобы провести одно суммирование логического теста, состоящего из двух частей, для N элементов, которые, как я подозреваю, представляют собой O(N*3). Это может быть немного быстрее, если реализовано в Rcpp, но это векторизованные операции.

# Wrong: ecdf2d <- function(m,i,j) { ord <- rank(m[ , 1]^2+m[ , 2]^2)
#            ord[i]/nrow(m)}  # scales to [0,1] interval

ecdf2d.v2 <- function(obj, x, y) sum( obj[,1] < x & obj[,2] < y)/nrow(obj)