Манипуляции с большой размерной матрицей в R

Question

Я хотел бы рассчитать следующую функцию для каждой строки матрицы M размерности 3e + 07x4.

func <- function(x){

   (dmultinom(c(x[c(1,2)],50-sum(x[c(1,2)])), size = NULL, rep(1/3,3), log = FALSE))/(x[3]^2+x[4]^3)

}

Я использую следующий код

as.numeric(unlist(apply(M, 1, function(v) func(v))))

К сожалению, это занимает много времени. Я хотел бы сделать это в ближайшее время.

Answer 1

Примерно так, но я не проверил правильность. Ваш:

func <- function(x) {
  (dmultinom(c(x[c(1,2)],50-sum(x[c(1,2)])), size = NULL, rep(1/3,3), log = FALSE))/(x[3]^2+x[4]^3)
}

может быть записан как:

func <- function(x) {
  a <- x[c(1,2)]
  b <- 50 - (a[1] + a[2])
  d <- c(a, b)
  e <- dmultinom(d, size = NULL, rep(1/3,3), log = FALSE)
  f <- x[3]^2 + x[4]^3
  e / f
}

Часть d, которую вы можете векторизовать с помощью матричных вычислений, как:

A <- M[, 1:2]
B <- 50 - (A[,1] + A[,2])
D <- cbind(A, B)

Без погружения в dmultinom(), часть e может быть вычислена через apply() как:

prob <- rep(1/3, times = 3L)
E <- apply(D, MARGIN = 1L, function(d) {
  dmultinom(d, size = NULL, prob = prob, log = FALSE)
})

Часть f, которую вы можете векторизовать с помощью матричных вычислений, как:

F <- M[,3]^2 + M[,4]^3

, которая дает следующее:

Y <- apply(M, 1, function(v) func(v))

может быть записано как:

Y <- E / F

Отказ от ответственности: не проверено, но вы должны получить представление о том, как векторизовать и избежать дублирования вещей.

PS. Если вы посмотрите на dmultinom(), я думаю, вы также можете векторизовать это аналогичным образом. Не исключено, что вы сможете избавиться от оставшегося звонка apply().

Answer 2

К счастью, lgamma - это примитивная функция, и поэтому существует возможность векторизации dmultinom самостоятельно. Здесь опция в сочетании с data.table для более быстрой скорости

set.seed(0L)
nr <- 3e7 #3e7
size <- 50L
DT <- data.table(X1=sample(1:20, nr, TRUE), X2=sample(1:20, nr, TRUE), X3=3, X4=4)

system.time({
    DT[, paste0("lgX", 1L:3L) := c(lapply(1+.SD, lgamma), .(lgamma(1+size-X1-X2))), .SDcols=X1:X2][,
       dmn := exp(lgamma(size + 1L) + log(1/3) * size - (lgX1 + lgX2 + lgX3)) / (X3^2 + X4^2)]
    DT$dmn
})

#   user  system elapsed 
#   7.44    0.17    7.64