计算时间序列矩阵的近似熵
引入近似熵来量化时间序列中的规律性和波动的不可预测性.
功能
approx_entropy(ts, edim = 2, r = 0.2*sd(ts), elag = 1)
从包中pracma,计算出时间序列的近似熵ts.
我有一个时间序列矩阵(每行一个系列)mat,我会估计每个矩阵的近似熵,将结果存储在一个向量中.例如:
library(pracma)
N<-nrow(mat)
r<-matrix(0, nrow = N, ncol = 1)
for (i in 1:N){
r[i]<-approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
}
但是,如果N很大,这段代码可能会太慢.建议加快速度吗？谢谢!
我还要说并行化，因为应用函数显然没有带来任何优化。
我尝试了该approx_entropy()功能：
- 申请
- 拉普利
- 申请帕利
- foreach（来自@Mankind_008）
- data.table 和 ParApply 的组合
这ParApply似乎比其他两个并行函数稍微高效一些。
由于我没有得到与@Mankind_008相同的时间，所以我用microbenchmark. 这些是 10 次运行的结果：
Run Code Online (Sandbox Code Playgroud)Unit: seconds expr min lq mean median uq max neval cld forloop 4.067308 4.073604 4.117732 4.097188 4.141059 4.244261 10 b apply 4.054737 4.092990 4.147449 4.139112 4.188664 4.246629 10 b lapply 4.060242 4.068953 4.229806 4.105213 4.198261 4.873245 10 b par 2.384788 2.397440 2.646881 2.456174 2.558573 4.134668 10 a parApply 2.289028 2.300088 2.371244 2.347408 2.369721 2.675570 10 a DT_parApply 2.294298 2.322774 2.387722 2.354507 2.466575 2.515141 10 a
完整代码：
library(pracma)
library(foreach)
library(parallel)
library(doParallel)
# dummy random time series data
ts <- rnorm(56)
mat <- matrix(rep(ts,100), nrow = 100, ncol = 100)
r <- matrix(0, nrow = nrow(mat), ncol = 1)
## For Loop
for (i in 1:nrow(mat)){
r[i]<-approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
}
## Apply
r1 = apply(mat, 1, FUN = function(x) approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1))
## Lapply
r2 = lapply(1:nrow(mat), FUN = function(x) approx_entropy(mat[x,], edim = 2, r = 0.2*sd(mat[x,]), elag = 1))
## ParApply
cl <- makeCluster(getOption("cl.cores", 3))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
## Foreach
registerDoParallel(cl = 3, cores = 2)
r4 <- foreach(i = 1:nrow(mat), .combine = rbind) %dopar%
pracma::approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
stopImplicitCluster()
## Data.table
library(data.table)
mDT = as.data.table(mat)
cl <- makeCluster(getOption("cl.cores", 3))
r5 = parApply(cl = cl, mDT, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
## All equal Tests
all.equal(as.numeric(r), r1)
all.equal(r1, as.numeric(do.call(rbind, r2)))
all.equal(r1, r3)
all.equal(r1, as.numeric(r4))
all.equal(r1, r5)
## Benchmark
library(microbenchmark)
mc <- microbenchmark(times=10,
forloop = {
for (i in 1:nrow(mat)){
r[i]<-approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
}
},
apply = {
r1 = apply(mat, 1, FUN = function(x) approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1))
},
lapply = {
r1 = lapply(1:nrow(mat), FUN = function(x) approx_entropy(mat[x,], edim = 2, r = 0.2*sd(mat[x,]), elag = 1))
},
par = {
registerDoParallel(cl = 3, cores = 2)
r_par <- foreach(i = 1:nrow(mat), .combine = rbind) %dopar%
pracma::approx_entropy(mat[i,], edim = 2, r = 0.2*sd(mat[i,]), elag = 1)
stopImplicitCluster()
},
parApply = {
cl <- makeCluster(getOption("cl.cores", 3))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
DT_parApply = {
mDT = as.data.table(mat)
cl <- makeCluster(getOption("cl.cores", 3))
r5 = parApply(cl = cl, mDT, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
}
)
## Results
mc
Unit: seconds
expr min lq mean median uq max neval cld
forloop 4.067308 4.073604 4.117732 4.097188 4.141059 4.244261 10 b
apply 4.054737 4.092990 4.147449 4.139112 4.188664 4.246629 10 b
lapply 4.060242 4.068953 4.229806 4.105213 4.198261 4.873245 10 b
par 2.384788 2.397440 2.646881 2.456174 2.558573 4.134668 10 a
parApply 2.289028 2.300088 2.371244 2.347408 2.369721 2.675570 10 a
DT_parApply 2.294298 2.322774 2.387722 2.354507 2.466575 2.515141 10 a
## Time-Boxplot
plot(mc)
核心数量也会影响速度，更多并不总是更快，因为在某些时候，发送给所有工作人员的开销会消耗掉一些获得的性能。我ParApply用 2 到 7 个核心对该函数进行了基准测试，在我的机器上，用 3 / 4 个核心运行该函数似乎是最佳选择，尽管偏差不是那么大。
Run Code Online (Sandbox Code Playgroud)mc Unit: seconds expr min lq mean median uq max neval cld parApply_2 2.670257 2.688115 2.699522 2.694527 2.714293 2.740149 10 c parApply_3 2.312629 2.366021 2.411022 2.399599 2.464568 2.535220 10 a parApply_4 2.358165 2.405190 2.444848 2.433657 2.485083 2.568679 10 a parApply_5 2.504144 2.523215 2.546810 2.536405 2.558630 2.646244 10 b parApply_6 2.687758 2.725502 2.761400 2.747263 2.766318 2.969402 10 c parApply_7 2.906236 2.912945 2.948692 2.919704 2.988599 3.053362 10 d
完整代码：
## Benchmark N-Cores
library(microbenchmark)
mc <- microbenchmark(times=10,
parApply_2 = {
cl <- makeCluster(getOption("cl.cores", 2))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
parApply_3 = {
cl <- makeCluster(getOption("cl.cores", 3))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
parApply_4 = {
cl <- makeCluster(getOption("cl.cores", 4))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
parApply_5 = {
cl <- makeCluster(getOption("cl.cores", 5))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
parApply_6 = {
cl <- makeCluster(getOption("cl.cores", 6))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
parApply_7 = {
cl <- makeCluster(getOption("cl.cores", 7))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
}
)
## Results
mc
Unit: seconds
expr min lq mean median uq max neval cld
parApply_2 2.670257 2.688115 2.699522 2.694527 2.714293 2.740149 10 c
parApply_3 2.312629 2.366021 2.411022 2.399599 2.464568 2.535220 10 a
parApply_4 2.358165 2.405190 2.444848 2.433657 2.485083 2.568679 10 a
parApply_5 2.504144 2.523215 2.546810 2.536405 2.558630 2.646244 10 b
parApply_6 2.687758 2.725502 2.761400 2.747263 2.766318 2.969402 10 c
parApply_7 2.906236 2.912945 2.948692 2.919704 2.988599 3.053362 10 d
## Plot Results
plot(mc)
随着矩阵变大，使用ParApplywithdata.table似乎比使用矩阵更快。以下示例使用具有 500*500 个元素的矩阵，从而得出这些计时（仅适用于 2 次运行）：
Run Code Online (Sandbox Code Playgroud)Unit: seconds expr min lq mean median uq max neval cld ParApply 191.5861 191.5861 192.6157 192.6157 193.6453 193.6453 2 a DT_ParAp 135.0570 135.0570 163.4055 163.4055 191.7541 191.7541 2 a
尽管最大值几乎相同，但最小值要低得多，该箱线图也很好地说明了这一点：
完整代码：
# dummy random time series data
ts <- rnorm(500)
# mat <- matrix(rep(ts,100), nrow = 100, ncol = 100)
mat = matrix(rep(ts,500), nrow = 500, ncol = 500, byrow = T)
r <- matrix(0, nrow = nrow(mat), ncol = 1)
## Benchmark
library(microbenchmark)
mc <- microbenchmark(times=2,
ParApply = {
cl <- makeCluster(getOption("cl.cores", 3))
r3 = parApply(cl = cl, mat, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
},
DT_ParAp = {
mDT = as.data.table(mat)
cl <- makeCluster(getOption("cl.cores", 3))
r5 = parApply(cl = cl, mDT, 1, FUN = function(x) {
library(pracma);
approx_entropy(x, edim = 2, r = 0.2*sd(x), elag = 1)
})
stopCluster(cl)
}
)
## Results
mc
Unit: seconds
expr min lq mean median uq max neval cld
ParApply 191.5861 191.5861 192.6157 192.6157 193.6453 193.6453 2 a
DT_ParAp 135.0570 135.0570 163.4055 163.4055 191.7541 191.7541 2 a
## Plot
plot(mc)