使用VCM提高我的R函数性能

Question

我正在使用变系数模型运行模拟，但有一些调整。没有R包可以做我正在寻找的。我的代码运行得不够快。我期待着让vcm函数运行得更快

###########################################################################
###########################################################################
###                                                                     ###
###                        EPANECHNIKOV FUNCTION                        ###
###                                                                     ###
###########################################################################
###########################################################################
epan <- function(t,h){
idx    = 0.75 * (1 - (t/h)**2) / h
kernal = 0.50 * (abs(idx) + idx)
kernal
}
###########################################################################
###########################################################################
###                                                                     ###
###                             UNPENALIZED                             ###
###                      VARYING COEFFICIENT MODEL                      ###
###                                                                     ###
###########################################################################
###########################################################################

vcm <- function(x,y,z,z0) {
n = dim(x)[1]
p = dim(x)[2]

n0 = length(z0)

Z = outer(z,z0,"-")
Width = sd(z) * n**(-0.2) * 2
H = sapply(X = 1:n0, FUN = function(X) epan(t = Z[,X], h = Width))
diag(H) = 0
W_h = H / rep(colSums(H), each = n0)

G = lapply(X = 1:n0, FUN = function(X) cbind(x, Z[,X]*x))

AB = matrix(NA, n0, 2*p)
II = 1e-4 * diag(2*p)  # to avoid singularity
for(i in 1:n0) {
AB[i,] = solve(crossprod(G[[i]] * W_h[,i], G[[i]]) + II) %*% crossprod(G[[i]] * W_h[,i], y)
}

AB
}

到目前为止我所做的是

• 配置代码并查看慢的部分
• 用sapply和lapply代替for环，但差异不显著

如何使用代码?下面是一个使用代码函数的小模拟。

n = 100000
p = 5
n0 = 1000
z = runif(n)
z0 = seq(0.05, 0.95, length.out = n0)
x = MASS::mvrnorm(n, rep(0,p), diag(p))
gz = cbind(2*sin(2*pi*z), 3*z*(1-2*z), exp(-2*z + z**2), 2*z, 0)
y = apply(x * gz, 1, sum) + rnorm(n)
vvc_m = vcm(x,y,z,z0)

我愿意使用Rcpp或任何其他库，如果它们能显著提高我的代码性能，即使我没有使用Rcpp的经验。

感谢您的帮助!

Answer 1

不需要apply和lapply。而且，G[[i]] * W_h[,i]只需要计算一次。这些更改将减少几秒钟的时间，但大部分时间都花在for循环上。你可能是对的，任何收益都必须与Rcpp/RcppArmadillo有关。

vcm2 <- function(x,y,z,z0) {
n = dim(x)[1]
p = dim(x)[2]

n0 = length(z0)

Z = outer(z,z0,"-")
Width = sd(z) * n**(-0.2) * 2
H = epan(Z, Width)
diag(H) = 0
W_h = H / rep(colSums(H), each = n0)

AB = matrix(NA, n0, 2*p)
II = 1e-4 * diag(2*p)  # to avoid singularity
G = matrix(x, n, 2*p)
idx = (p + 1):(2*p)
for(i in 1:n0) {
G[,idx] = Z[,i]*x
GW_h = G*W_h[,i]
AB[i,] = solve(crossprod(GW_h, G) + II) %*% crossprod(GW_h, y)
}

AB
}
system.time(vvc_m <- vcm(x,y,z,z0))
#>    user  system elapsed 
#>   21.71    5.42   27.14
system.time(vvc_m2 <- vcm2(x,y,z,z0))
#>    user  system elapsed 
#>   19.45    3.52   22.99
identical(vvc_m, vvc_m2)
#> [1] TRUE