Ответ 1

Выяснил это.

# Step1 - map each sample to the nearest X, Y, and Z above it. (In other words, bin the data.)

X <- data.table(X=seq(-2, 2, by=.1)); X[, x := X]
Y <- data.table(Y=seq(-2, 2, by=.1)); Y[, y := Y]
Z <- data.table(Z=seq(-2, 2, by=.1)); Z[, z := Z]

dt <- X[dt, on="x", roll=-Inf, nomatch=0]
dt <- Y[dt, on="y", roll=-Inf, nomatch=0]
dt <- Z[dt, on="z", roll=-Inf, nomatch=0]

# Step2 - aggregate by unique (X, Y, Z) triplets and count the samples directly below each of these bounds.
bg <- dt[, .N, keyby=list(X, Y, Z)]

# Step4 - Get the count of samples directly below EVERY (X, Y, Z) bound
bounds <- CJ(X=X$X, Y=Y$Y, Z=Z$Z)
kl <- bg[bounds, on=c("X", "Y", "Z")]
kl[is.na(N), N := 0]

# Step5 (the tricky part) - Consider a single (Y, Z) pair. X will be in ascending order. So we can do a cumsum on X for each (Y, Z) to count x <= X | Y,Z. Now if you hold X and Z fixed, you can do a cumsum on Y (which is also in ascending order) to count x <= X, y <= Y | Z. And then just continue this process.
kl[, CountUntil.XgivenYZ := cumsum(N), by=list(Y, Z)]
kl[, CountUntil.XYgivenZ := cumsum(CountUntil.XgivenYZ), by=list(X, Z)]
kl[, CountUntil.XYZ := cumsum(CountUntil.XYgivenZ), by=list(X, Y)]

# Cleanup
setnames(kl, "CountUntil.XYZ", "N.cum")
kl[, CDF := N.cum/nrow(dt)]

Обобщение

Для тех, кто этого хочет, я обобщил это для работы с любым количеством переменных и сбрасывал функцию в мой пакет R, mltools.

Например, чтобы решить эту проблему, вы можете сделать

library(mltools)

bounds <- list(x=seq(-2, 2, by=.1), y=seq(-2, 2, by=.1), z=seq(-2, 2, by=.1))
empirical_cdf(x=dt, ubounds=bounds)
        x  y    z N.cum   CDF
    1: -2 -2 -2.0     0 0.000
    2: -2 -2 -1.9     0 0.000
    3: -2 -2 -1.8     0 0.000
    4: -2 -2 -1.7     0 0.000
    5: -2 -2 -1.6     0 0.000
   ---                       
68917:  2  2  1.6   899 0.899
68918:  2  2  1.7   909 0.909
68919:  2  2  1.8   917 0.917
68920:  2  2  1.9   924 0.924
68921:  2  2  2.0   929 0.929

Ответ 2

Обновление

Ниже я предоставил общее решение base R (оно будет работать на неравномерных сетках). Он был быстрее, чем самое быстрое опубликованное решение, предоставляемое OP (подробнее об этом позже). Поскольку ОП указывает, генерация столбца N.cum является настоящим узким местом, поэтому я сосредоточил свои усилия только на этой задаче (т.е. Создание CDF является тривиальной задачей после получения N.cum).

JoeBase <- function(dtt, s) {
    m <- matrix(c(dtt$x, dtt$y, dtt$z), ncol = 3)
    N.Cum <- array(vector(mode = "integer"), dim = rev(sapply(s, length)))
    for (i in seq_along(s[[1]])) {
        t1 <- m[,1] <= s[[1]][i]
        for (j in seq_along(s[[2]])) {
            t2 <- t1 & (m[,2] <= s[[2]][j])
            for (k in seq_along(s[[3]])) {
                N.Cum[k,j,i] <- sum(t2 & (m[,3] <= s[[3]][k]))
            }
        }
    }
    as.vector(N.Cum)
}

В приведенном выше алгоритме используются векторизованные операции, в частности создание и использование логических векторов t1 и t2. Этот вектор используется для получения числа строк, удовлетворяющих критериям для всех 3 столбцов в исходной таблице данных. Мы просто полагаемся на внутреннее принуждение через R от логического вектора к интегральному вектору действием sum.

Выяснение того, как заполнить трехмерный массив целых чисел N.cum, было немного сложной задачей, так как позже он будет преобразован в вектор через as.vector. Для изучения поведения as.vector потребовалось немного проб и ошибок. К моему удивлению, "последнее" и "первое" измерение должно быть перестроено для того, чтобы принуждение к вектору происходило верно (первые несколько раз, я использовал N.Cum [i, j, k] вместо N.Cum [K, J, I]).

Сначала давайте проверим равенство:

library(data.table)
## Here is the function I used to test against. I included the generation
## of "bounds" and "bg" as "result" depends on both of these (N.B. "JoeBase" does not)
BenDT <- function(dt, s) {
    X <- data.table(X=s[[1]]); X[, x := X]
    Y <- data.table(Y=s[[2]]); Y[, y := Y]
    Z <- data.table(Z=s[[3]]); Z[, z := Z]

    dt <- X[dt, on="x", roll=-Inf, nomatch=0]
    dt <- Y[dt, on="y", roll=-Inf, nomatch=0]
    dt <- Z[dt, on="z", roll=-Inf, nomatch=0]
    bg <- dt[, .N, keyby=list(X, Y, Z)]

    bounds <- CJ(X=X$X, Y=Y$Y, Z=Z$Z)

    kl <- bg[bounds, on=c("X", "Y", "Z")]
    kl[is.na(N), N := 0]

    # Counting
    kl[, CountUntil.XgivenYZ := cumsum(N), by=list(Y, Z)]
    kl[, CountUntil.XYgivenZ := cumsum(CountUntil.XgivenYZ), by=list(X, Z)]
    kl[, CountUntil.XYZ := cumsum(CountUntil.XYgivenZ), by=list(X, Y)]

    # Cleanup
    setnames(kl, "CountUntil.XYZ", "N.cum")
    kl[, CDF := N.cum/nrow(dt)]
    kl
}

t1 <- BenDT(dt, seq(-2,2,0.1))
t2 <- JoeBase(dt, seq(-2,2,0.1))

all.equal(t1$N.cum, t2)
[1] TRUE

Теперь мы проверяем скорость. Сначала мы скомпилируем обе функции с помощью cmpfun из пакета compiler. Первый контрольный показатель отражает эффективность на более мелких примерах.

library(compiler)
c.JoeBase <- cmpfun(JoeBase)
c.BenDT <- cmpfun(BenDT)
c.OldBenDT <- cmpfun(OldBenDT)  ## The previous best solution that Ben contributed

st <- list(seq(-2, 2, 0.1), seq(-2, 2, 0.1), seq(-2, 2, 0.1))
microbenchmark(c.BenDT(dt, st), c.OldBenDT(dt, st), c.JoeBase(dt, st), times = 10)
Unit: milliseconds
              expr        min         lq       mean    median         uq        max neval cld
   c.BenDT(dt, st)   34.24872   34.78908   38.87775   37.4924   43.37179   46.12859    10 a  
c.OldBenDT(dt, st) 1485.68178 1532.35878 1607.96669 1593.9813 1619.58908 1845.75876    10  b 
 c.JoeBase(dt, st) 1880.71648 1962.38160 2049.43985 2007.4880 2169.93078 2281.02118    10   c

Ниже приведен старый тест.
Однако, когда количество ящиков увеличивается, c.JoeBase действительно начинает доминировать (более чем в 5 раз быстрее).

st <- list(seq(-5, 5, 0.1), seq(-5, 5, 0.1), seq(-5, 5, 0.1))
microbenchmark(c.JoeBase(dt, st), c.OldBenDT(dt, st), times = 5)
Unit: seconds
              expr       min        lq      mean    median        uq       max neval cld
 c.JoeBase(dt, st)  23.50927  23.53809  29.61145  24.52748  30.81485  45.66759     5  a 
c.OldBenDT(dt, st) 110.60209 123.95285 133.74601 124.97929 125.96186 183.23394     5   b

После проведения дальнейших тестов у меня есть некоторые опасения относительно результатов (@Ben указал на подобное мнение в комментариях). Я уверен, что c.JoeBase работает быстрее только из-за ограничений моего старого компьютера. Как отметил в своем ответе @stephematician, исходное решение интенсивно запоминается, и если вы просто выполните system.time на c.OldBenDT, вы увидите, что большую часть времени тратится в категории system и категория user сравнима с категорией user c.JoeBase. Мой 6-летний Mac имеет только 4 ГБ оперативной памяти, и я предполагаю, что с этими операциями происходит большое количество свопов памяти. Обратите внимание:

## test with very tiny buckets (i.e. 0.025 instead of 0.1 above)
st <- list(seq(-1.5, 1.5, 0.025), seq(-1.5, 1.5, 0.025), seq(-1.5, 1.5, 0.025))  
system.time(c.JoeBase(dt, st))
  user  system elapsed 
36.407   4.748  41.170

system.time(c.OldBenDT(dt, st))
   user  system elapsed 
49.653  77.954 475.304

system.time(c.BenDT(dt, st))  ## Ben new solution is lightning fast
  user  system elapsed 
 0.603   0.063   0.668

Несмотря на это, последнее решение @Ben намного превосходит. Ознакомьтесь с этими новыми критериями:

st <- list(seq(-5, 5, 0.1), seq(-5, 5, 0.1), seq(-5, 5, 0.1))
microbenchmark(c.JoeBase(dt, st), BenDT(dt, st), times = 5)
Unit: milliseconds
             expr        min         lq       mean     median        uq        max neval cld
c.JoeBase(dt, st) 26517.0944 26855.7819 28341.5356 28403.7871 29926.213 30004.8018     5   b
    BenDT(dt, st)   342.4433   359.8048   400.3914   379.5319   423.336   496.8411     5  a

Еще одна победа data.table.

Ответ 3

Просто примечание об альтернативе, и все же очевидное решение:

set.seed(2016)
dt <- data.table(x=rnorm(20000), y=rnorm(20000), z=rnorm(20000))

system.time({
    dt <- t(as.matrix(dt))

    bounds <- as.matrix(expand.grid(z=seq(-2,2,0.1),
                                    y=seq(-2,2,0.1),
                                    x=seq(-2,2,0.1)))

    bounds <- bounds[,ncol(bounds):1]

    n_d <- ncol(bounds)

    x <- apply(bounds,
               1,
               function(x) sum(colSums(dt < x) == n_d))
})

Это решение на моей машине занимает примерно в два раза больше времени, чтобы рассчитать решения JoeBase и OldBenDT. Основное различие? Использование памяти. Он больше связан с процессором.

Я не знаю точного способа сравнения использования памяти в R, но функция memory.size(max=T) сообщила, используя 5 ГБ памяти для этих предыдущих подходов (а не для подхода без привязки), в то время как только с использованием 40 Мб памяти для подхода apply (примечание: я использовал 20000 точек в распределении выборки).

Я думаю, что это имеет важные последствия для масштаба вычислений, которые вы можете выполнить.

Ответ 4

Необходимо быстрее вычислять пропорции и делать соединения за один шаг, поэтому промежуточные результаты не должны быть реализованы:

set.seed(2016)
dt <- data.table(x=rnorm(1000), y=rnorm(1000), z=rnorm(1000))
setkey(dt)

bounds <- CJ(x=seq(-2, 2, by=.1), y=seq(-2, 2, by=.1), z=seq(-2, 2, by=.1))

a <- dt[bounds,.N / nrow(dt),on=c("x<x","y<y","z<z"),
        by=.EACHI,
        allow.cartesian=T]