2

给定一个 n*2 数据矩阵 X 我想计算每个观察的二元经验 cdf，即对于 1:n 中的每个 i，返回第一个元素不大于 X[i,1] 和第二个的观察百分比元素不大于 X[i,2]。

由于涉及嵌套搜索，即使在将其移植到 Fortran 之后，它在 n ~ 100k 时也会变得非常慢。有谁知道是否有更好的方法来处理这样的样本大小？

编辑：我相信这个问题（就复杂性而言）类似于找到 Kendall 的 tau，它的顺序为 O(n^2)。在这种情况下，Knight (1966) 有一个算法可以将其减少到 O(n log(n))。只是想知道是否有任何 O(n*log(n)) 算法可以找到已经存在的二元 ecdf。

编辑 2：这是我在 Fortran 中的代码，根据要求。这在 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
4

1 回答 1

0

我认为我的第一个努力并不是真正的 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))

因此，任务是对我怀疑是 O(N*3) 的 N 个项目进行两部分逻辑测试的单一总和。如果在 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)

在此处输入图像描述

于 2015-12-19T22:46:12.097 回答