MCP*_*tor 10 performance r vectorization
我正在尝试对以下函数进行矢量化以删除 sapply 循环。我正在计算累积偏度。
cskewness <- function(.x) {
skewness <- function(.x) {
sqrt(length(.x)) * sum((.x - mean(.x))^3) / (sum((.x - mean(.x))^2)^(3 / 2))
}
sapply(seq_along(.x), function(k, z) skewness(z[1:k]), z = .x)
}
我的代数没搞对。有这个是错误的：
skewness2 <- function(.x) {
n <- length(.x)
csum <- cumsum(.x)
cmu <- csum / 1:length(.x)
num <- cumsum(.x - cmu)^3
den <- cumsum((.x - cmu)^2)^(3/2)
sqrt(n) * num / den
}
正确的代码会产生：
x <- c(1,2,4,5,8)
> cskewness(x)
[1] NaN 0.0000000 0.3818018 0.0000000 0.4082483
> skewness2(x)
[1] NaN 1.000000 1.930591 3.882748 4.928973
感谢jblood94努力捕捉重大异常值可能造成的危险，让我有机会重新思考和优化代码。
我们可以用它scale来预处理输入向量，这样可以最小化异常值的有害影响，并且偏度也不会受到影响
cskewness_tic <- function(y) {
# rescale `y` to avoid detrimental impacts by outliers
y <- scale(y)
# cumulative length of y
k <- seq_along(y)
# cumulative n-th raw moments of y
m3 <- cumsum(y^3)
m2 <- cumsum(y^2)
m1 <- cumsum(y)
u <- m1 / k
# expand cubic terms and refactor skewness in terms of num/den
num <- (m3 - 3 * u * m2 + 3 * u^2 * m1 - k * u^3) / k
den <- sqrt((m2 + k * u^2 - 2 * u * m1) / k)^3
c(NaN, (num / den)[-1])
}
包含用于基准测试的函数列表（请参阅@jblood94 的答案的详细cumskewness信息）cumskewnessCpp
set.seed(0)
x <- sample(1e6) + 1e10
microbenchmark(
cskewness_tic = cskewness_tic(x),
cumskewness = cumskewness(x),
cumskewnessCpp = cumskewnessCpp(x),
check = "equal",
unit = "relative",
times = 50L
)
节目
Unit: relative
expr min lq mean median uq max neval
cskewness_tic 3.144719 3.268279 3.461112 3.314326 3.846614 2.395820 50
cumskewness 4.647973 4.627621 4.455523 4.585634 4.448119 2.687571 50
cumskewnessCpp 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 50
首先，非矢量化方法并不像您想象的那么糟糕，因此*apply如果您认为它符合您的最终目标并简化您的工作，那么它也是一个不错的选择。无需带着恐惧面对循环。
如果您想从矢量化中受益（例如速度），您可以尝试根据原始矩扩展偏度的数学表达式（请参阅this和this），然后cumsum对每个解耦项进行运算
cskewness_tic <- function(y) {
# cumulative length of y
k <- seq_along(y)
# cumulative n-th raw moments of y
m3 <- cumsum(y^3)
m2 <- cumsum(y^2)
m1 <- cumsum(y)
u <- m1 / k
# expand cubic terms and refactor skewness in terms of num/den
num <- (m3 - 3 * u * m2 + 3 * u^2 * m1 - k * u^3) / k
den <- sqrt((m2 + k * u^2 - 2 * u * m1) / k)^3
num / den
}
这样
> cskewness_tic(x)
[1] NaN 0.0000000 0.3818018 0.0000000 0.4082483
方法列表
cskewness_origin <- function(.x) {
skewness <- function(.x) {
sqrt(length(.x)) * sum((.x - mean(.x))^3) / (sum((.x - mean(.x))^2)^(3 / 2))
}
sapply(seq_along(.x), function(k, z) skewness(z[1:k]), z = .x)
}
cskewness_jblood94 <- function(.x) {
d <- outer(.x, cumsum(.x) / (1:length(.x)), "-")
d[lower.tri(d)] <- 0
sqrt(1:length(.x)) * colSums(d^3) / colSums(d^2)^(3 / 2)
}
cskewness_tic <- function(y) {
# cumulative length of y
k <- seq_along(y)
# cumulative n-th order moments of y
m3 <- cumsum(y^3)
m2 <- cumsum(y^2)
m1 <- cumsum(y)
u <- m1 / k
# expand cubic terms and refactor skewness in terms of num/den
num <- (m3 - 3 * u * m2 + 3 * u^2 * m1 - k * u^3) / k
den <- sqrt((m2 + k * u^2 - 2 * u * m1) / k)^3
num / den
}
通过运行以下基准测试
set.seed(0)
x <- sample(1e3)
microbenchmark(
origin = cskewness_origin(x),
jblood94 = cskewness_jblood94(x),
tic = cskewness_tic(x),
unit = "relative",
check = "equivalent",
times = 50L
)
我们看
Unit: relative
expr min lq mean median uq max neval
origin 252.7504 247.1442 253.5965 254.8663 257.2905 239.0963 50
jblood94 302.8267 309.0854 326.7546 318.7141 310.2365 388.3786 50
tic 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 50
使用单遍for循环将是高效的：
cumskewness <- function(x) {
n <- length(x)
if (n == 0L) {
return(x)
} else if (n == 1L) {
return(0)
}
m2 <- m3 <- term1 <- 0
out <- numeric(n)
out[1] <- NaN
m1 <- x[1]
for (i in 2:n) {
n0 <- i - 1
delta <- x[i] - m1
delta_n <- delta/i
m1 <- m1 + delta_n
term1 <- delta*delta_n*n0
m3 <- m3 + term1*delta_n*(n0 - 1) - 3*delta_n*m2
m2 <- m2 + term1
out[i] <- sqrt(i)*m3/m2^1.5
}
out
}
写起来也很简单Rcpp：
library(Rcpp)
cppFunction('
NumericVector cumskewnessCpp(const NumericVector& x) {
const int n = x.size();
if (n == 0) {
return(x);
} else if (n == 1) {
return(0);
}
int n1;
double m1 = x(0);
double m2, m3, delta, delta_n, term1;
NumericVector out(n);
out(0) = R_NaN;
for (int i = 1; i < n; i++) {
n1 = i + 1;
delta = x(i) - m1;
delta_n = delta/n1;
m1 += delta_n;
term1 = delta*delta_n*i;
m3 += term1*delta_n*(i - 1) - 3*delta_n*m2;
m2 += term1;
out(i) = sqrt(n1)*m3/pow(m2, 1.5);
}
return out;
}
')
cumskewness(x)
#> [1] NaN 0.0000000 0.3818018 0.0000000 0.4082483
cumskewnessCpp(x)
#> [1] NaN 0.000000e+00 3.818018e-01 -2.808667e-17 4.082483e-01
包括来自 @ThomisIsCoding 的矢量化解决方案：
cskewness_tic <- function(y) {
# cumulative length of y
k <- seq_along(y)
# cumulative n-th order moments of y
m3 <- cumsum(y^3)
m2 <- cumsum(y^2)
m1 <- cumsum(y)
u <- m1 / k
# expand cubic terms and refactor skewness in terms of num/den
num <- (m3 - 3 * u * m2 + 3 * u^2 * m1 - k * u^3) / k
den <- sqrt((m2 + k * u^2 - 2 * u * m1) / k)^3
num / den
}
set.seed(0)
x <- sample(1e3)
microbenchmark::microbenchmark(
cskewness_tic = cskewness_tic(x),
cumskewness = cumskewness(x),
cumskewnessCpp = cumskewnessCpp(x),
check = "equal",
unit = "relative"
)
#> Unit: relative
#> expr min lq mean median uq max neval
#> cskewness_tic 2.035272 2.07954 2.006118 2.388216 2.282022 0.7228815 100
#> cumskewness 3.930424 3.96835 4.003956 3.939762 3.711879 1.1339946 100
#> cumskewnessCpp 1.000000 1.00000 1.000000 1.000000 1.000000 1.0000000 100
cskewness_tic灾难性的抵消可能会导致高阶矩的精度误差。当标准差相对平均值较小时就会发生这种情况。展示：
set.seed(0)
x <- sample(1e3) + 1e8
y1 <- cskewness(x)
y2 <- cumskewness(x)
y3 <- cumskewnessCpp(x)
y4 <- cskewness_tic(x); y4[1] <- NaN
all.equal(y1, y2)
#> [1] TRUE
all.equal(y1, y3)
#> [1] TRUE
all.equal(y1, y4)
#> [1] "Mean relative difference: 270.1872"
skewness2 <- function(.x) {
d <- outer(.x, cumsum(.x)/(1:length(.x)), "-")
d[lower.tri(d)] <- 0
sqrt(1:length(.x))*colSums(d^3)/colSums(d^2)^(3/2)
}
x <- c(1,2,4,5,8)
cskewness(x)
#> [1] NaN 0.0000000 0.3818018 0.0000000 0.4082483
skewness2(x)
#> [1] NaN 0.0000000 0.3818018 0.0000000 0.4082483