У меня есть набор случайных величин, каждая из которых имеет среднее значение и стандартное отклонение. Я хочу найти все переменные, которые являются или не являются доминируемыми, т. е. имеют наивысшее среднее значение своих стандартных отклонений. Другими словами, если X1 имеет среднее значение 10 и стандартное отклонение 10, а X2 имеет среднее значение 9 и стандартное отклонение 11, то X2 доминирует над X1. По сути, я хочу отказаться от него. Это распространенная проблема при оптимизации портфеля.
Легко, когда скорость не имеет значения:
aredominated1 <- function( smean, svar ) {
PN <- length(smean)
kill <- rep(FALSE, PN)
for (i in 1:PN) {
for (j in 1:PN) {
if (i==j) next
kill[i] <- if ((smean[i] < smean[j]) & (svar[i] > svar[j])) 1 else kill[i]
if (kill[i]) break ## no longer any need to check whether i is dominated by a subsequent j
}
}
kill
}
aredominated2 <- function( smean, svar ) {
PN <- length(smean)
kill <- rep(FALSE, PN)
for (i in 1:PN) {
kill[i] <- any( (smean[i] < smean) & (svar[i] > svar) )
}
kill
}
####
T <- 100000
smean <- rnorm(T)
svar <- abs(rnorm(T))
print( system.time( { v1 <- aredominated1( smean, svar ) } ) )
print( system.time( { v2 <- aredominated2( smean, svar ) } ) )
stopifnot( all( v1 == v2 ) )
Грубо говоря, aredominated1
занимает 1,5 секунды, а aredominated2
(более R-подобный) занимает 60 секунд. Есть ли гораздо более быстрый способ написать такую функцию в R, или я близок к собственному пределу R?
Ваш первый метод быстрый, потому что вы не проверяете все возможности.
Вы можете сократить время вычислений еще больше (более чем вдвое), проверяя только часть «svar» при условии, что первая часть истинна, используя &&
вместо &
и опуская if
.
aredominated3 <- function( smean, svar ) {
PN <- length(smean)
kill <- rep(FALSE, PN)
for (i in 1:PN) {
for (j in 1:PN) {
if (i==j) next
## modify this line
kill[i] <- smean[i] < smean[j] && svar[i] > svar[j]
if (kill[i]) break ## no longer any need to check
}
}
kill
}
print( system.time( { v1 <- aredominated1( smean, svar ) } ) )
#user system elapsed
#2.66 0.00 3.02
print( system.time( { v3 <- aredominated3( smean, svar ) } ) )
#user system elapsed
#0.85 0.00 0.92
stopifnot( all( v1 == v3 ) )
Ваш текущий алгоритм — O(n^2). Вы можете перейти к O(n log(n)) сначала отсортировав данные.
library(data.table)
aredominated1_dt <- function(smean, svar) {
dt <- data.table(smean, svar)[, row := .I]
setorder(dt, -smean, svar)
dt[, kill := as.numeric(svar > cummin(svar))]
dt$kill[order(dt$row)]
}
data.table здесь не является строго необходимым. Логику можно легко перевести на базовый R или dplyr (но если производительность критична, я бы провел сравнение с коллапсом).
Тесты:
T <- 100000
smean <- rnorm(T)
svar <- abs(rnorm(T))
microbenchmark::microbenchmark(
aredominated1_dt(smean, svar),
aredominated1(smean, svar),
check = "identical",
times = 3
)
# Unit: milliseconds
# expr min lq mean median
# aredominated1_dt(smean, svar) 24.8443 24.8443 24.8443 24.8443
# aredominated1(smean, svar) 1801.3821 1801.3821 1801.3821 1801.3821
И чтобы получить представление о масштабировании, мы увеличиваем входные данные в четыре раза (T <- 400000).
# Unit: milliseconds
# expr min lq mean median
# aredominated1_dt(smean, svar) 60.3873 60.3873 60.3873 60.3873
# aredominated1(smean, svar) 9213.9352 9213.9352 9213.9352 9213.9352
Возможно, вы можете оптимизировать цикл, как показано ниже.
tic <- function(smean, svar) {
k <- order(smean)
v <- svar[k]
l <- length(v)
kill <- vector(length = l)
for (i in 1:(l - 1)) {
for (j in (i + 1):l) {
if (v[i] > v[j]) {
kill[i] <- 1
break
}
}
}
kill[order(k)]
}
ivo1 <- function(smean, svar) {
PN <- length(smean)
kill <- rep(FALSE, PN)
for (i in 1:PN) {
for (j in 1:PN) {
if (i == j) next
kill[i] <- if ((smean[i] < smean[j]) & (svar[i] > svar[j])) 1 else kill[i]
if (kill[i]) break ## no longer any need to check whether i is dominated by a subsequent j
}
}
kill
}
ivo2 <- function(smean, svar) {
PN <- length(smean)
kill <- rep(FALSE, PN)
for (i in 1:PN) {
kill[i] <- any((smean[i] < smean) & (svar[i] > svar))
}
+kill
}
tic <- function(smean, svar) {
k <- order(smean)
v <- svar[k]
l <- length(v)
kill <- vector(length = l)
for (i in 1:(l - 1)) {
for (j in (i + 1):l) {
if (v[i] > v[j]) {
kill[i] <- 1
break
}
}
}
kill[order(k)]
}
sbaldur <- function(smean, svar) {
dt <- data.table(smean, svar)[, row := .I]
setorder(dt, -smean, svar)
dt[, kill := as.numeric(svar > cummin(svar))]
dt$kill[order(dt$row)]
}
microbenchmark(
ivo1 = ivo1(smean, svar),
ivo2 = ivo2(smean, svar),
tic = tic(smean, svar),
sbaldur = sbaldur(smean, svar),
times = 10,
check = "equal",
unit = "relative"
)
шоу
Unit: relative
expr min lq mean median uq max
ivo1 22.371324 19.695451 17.558189 18.104484 18.385273 11.2895773
ivo2 439.136523 390.277007 340.305299 355.019021 333.328700 220.5935180
tic 1.562904 1.397621 1.255117 1.340289 1.262332 0.7891665
sbaldur 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000
neval
10
10
10
10
фу
sbaldur
функция просто невероятная...