时间戳

首页 > 解决方案 > 加速 R 代码以计算马尔可夫链中的（高阶）转换

r - 加速 R 代码以计算马尔可夫链中的（高阶）转换

问题描述

我写了一个小 R 片段来检查一个包含来自马尔可夫链的实现的向量，并返回给定顺序的观察到的转换。具体而言，假设我们对状态空间 $\mathcal{S}$ 的 2 阶转换感兴趣。最终目标是以方便的形式存储计数 $n_{ijk}$, $i, j, k \in \mathcal{S}$ 以供以后使用。

get_contingency_array <- function(seq, order){
  N <- length(seq)
  states <- sort(unique(seq))
  nstates <- length(states)
  inds <- seq_along(states)
  K <- order + 1
  out <- array(0, dim = rep(nstates, K))
  for (z in K:N){
    pos <- match(seq[(z - K + 1):z], states)
    s1 <- paste0("out[", paste(pos, collapse = ","), "]")
    s2 <- paste(s1, "<-", s1, "+ 1")
    eval(parse(text = s2))
  }
  return(out)
}

set.seed(666)
X1 <- rbinom(1E5, size = 1, prob = .8)
get_contingency_array(X1, 3)

给

, , 1, 1

     [,1] [,2]
[1,]  148  624
[2,]  613 2567

, , 2, 1

     [,1]  [,2]
[1,]  601  2583
[2,] 2550 10310

, , 1, 2

     [,1]  [,2]
[1,]  613  2527
[2,] 2578 10327

, , 2, 2

      [,1]  [,2]
[1,]  2590 10310
[2,] 10304 40752

我的问题是：如何才能大大加快速度？对于较大的尺寸，它会稍微拖动一点。

标签： rperformance

解决方案

除了最极端的情况外，eval(parse(.))应避免在所有情况下使用。

第一次切割使用一种鲜为人知的方法来索引带有matrix. 为了演示，我将在第一个for循环中中断函数调用，并在修改后的数组上显示索引：

debug(get_contingency_array)
get_contingency_array(X1, 3)
### step through until 'for' loop, first pass:
z
# [1] 4
### temporarily re-assign `out` just to show the indexing
out <- array(1:16, dim = rep(nstates, K))
out
# , , 1, 1
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4
# , , 2, 1
#      [,1] [,2]
# [1,]    5    7
# [2,]    6    8
# , , 1, 2
#      [,1] [,2]
# [1,]    9   11
# [2,]   10   12
# , , 2, 2
#      [,1] [,2]
# [1,]   13   15
# [2,]   14   16
pos <- match(seq[(z - K + 1):z], states)
pos
# [1] 2 2 1 2
matrix(pos, nrow = 1)
#      [,1] [,2] [,3] [,4]
# [1,]    2    2    1    2
out[matrix(pos, nrow = 1)]
# [1] 12

有了这个，我们可以实现相当大的速度提升。

get_contingency_array2 <- function(seq, order){
  N <- length(seq)
  states <- sort(unique(seq))
  nstates <- length(states)
  inds <- seq_along(states)
  K <- order + 1
  out <- array(0, dim = rep(nstates, K))
  for (z in K:N){
    pos <- matrix(match(seq[(z - K + 1):z], states), nrow = 1)
    out[pos] <- out[pos] + 1
  }
  return(out)
}

尽管维度名称不同，但另一种替代方法更快一些并产生相同的值：

get_contingency_array3 <- function(seq, order){
  N <- length(seq)
  states <- sort(unique(seq))
  nstates <- length(states)
  inds <- seq_along(states)
  K <- order + 1
  ind <- sapply(K:N, function(z) (z-K+1):z)
  ind[] <- seq[ind]
  do.call(table, asplit(ind, 1))
}

计时结果：

set.seed(666)
X1 <- rbinom(1E5, size = 1, prob = .8)
system.time(res1 <- get_contingency_array(X1, 3))
#    user  system elapsed 
#    5.83    0.01    5.86 
system.time(res2 <- get_contingency_array2(X1, 3))
#    user  system elapsed 
#     0.5     0.0     0.5 
system.time(res3 <- get_contingency_array3(X1, 3))
#    user  system elapsed 
#    0.16    0.00    0.17 
identical(res1, res2)
# [1] TRUE
all(res1 == res3)
# [1] TRUE

唯一的区别res3是它是类"table"并且具有不同的维度名称。如果您需要它们，解决这个问题并不难。

更彻底的时序基准：

bench::mark(
  get_contingency_array(X1, 3),
  get_contingency_array2(X1, 3),
  get_contingency_array3(X1, 3),
  check = FALSE
)
# Warning: Some expressions had a GC in every iteration; so filtering is disabled.
# # A tibble: 3 x 13
#   expression                         min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result memory                       time         gc              
#   <bch:expr>                    <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list> <list>                       <list>       <list>          
# 1 get_contingency_array(X1, 3)     7.25s    7.25s     0.138  215.76MB     2.07     1    15      7.25s <NULL> <Rprofmem[,3] [409,247 x 3]> <bch:tm [1]> <tibble [1 x 3]>
# 2 get_contingency_array2(X1, 3) 424.43ms 439.58ms     2.27     1.38MB     3.41     2     3   879.16ms <NULL> <Rprofmem[,3] [18,887 x 3]>  <bch:tm [2]> <tibble [2 x 3]>
# 3 get_contingency_array3(X1, 3) 196.88ms 207.82ms     4.81    37.66MB     4.81     3     3   624.32ms <NULL> <Rprofmem[,3] [944 x 3]>     <bch:tm [3]> <tibble [3 x 3]>

所以最后一个函数有显着的速度提升（在itr/sec现场清楚地显示）。

推荐阅读