Быстрое вычисление CDF/подвижного соединения на нескольких столбцах

Question

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

library(data.table) # v 1.9.7 

set.seed(2016) 
dt <- data.table(x=rnorm(1000), y=rnorm(1000), z=rnorm(1000)) 
dt 
      x  y  z 
    1: -0.91474 2.07025 -1.7499 
    2: 1.00125 -1.80941 -1.3856 
    3: -0.05642 1.58499 0.8110 
    4: 0.29665 -1.16660 0.3757 
    5: -2.79147 -1.75526 1.2851 
    ---       
996: 0.63423 0.13597 -2.3710 
997: 0.21415 1.03161 -1.5440 
998: 1.15357 -1.63713 0.4191 
999: 0.79205 -0.56119 0.6670 
1000: 0.19502 -0.05297 -0.3288 

Я хочу, чтобы подсчитать количество образцов, таких, что (х < = Х, у < = Y, Z < = Z) для некоторой сетки (X, Y, Z) верхние границы, как

bounds <- CJ(X=seq(-2, 2, by=.1), Y=seq(-2, 2, by=.1), Z=seq(-2, 2, by=.1)) 
bounds 
     X Y Z 
    1: -2 -2 -2.0 
    2: -2 -2 -1.9 
    3: -2 -2 -1.8 
    4: -2 -2 -1.7 
    5: -2 -2 -1.6 
    ---   
68917: 2 2 1.6 
68918: 2 2 1.7 
68919: 2 2 1.8 
68920: 2 2 1.9 
68921: 2 2 2.0 

Теперь я понял, что я элегантно могу сделать это (используя не-следу присоединяется)

dt[, Count := 1] 
result <- dt[bounds, on=c("x<=X", "y<=Y", "z<=Z"), allow.cartesian=TRUE][, list(N.cum = sum(!is.na(Count))), keyby=list(X=x, Y=y, Z=z)] 
result[, CDF := N.cum/nrow(dt)] 
result 
     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 

Но этот метод действительно неэффективно и становится очень медленным s Я начинаю увеличивать счетчик bin. Я думаю, что многовариантная версия функциональности скользящего соединения data.table сделала бы трюк, но это пока невозможно. Любые предложения по ускорению этого?

Answer 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 

Answer 2

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

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] 

Answer 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)) 
}) 

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

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

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

Answer 4

Update

Ниже я представил общий 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.

Выяснение того, как заполнить 3-мерный целочисленный массив 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's new solution is lightning fast 
    user system elapsed 
0.603 0.063 0.668 

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

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.