テンポラル分数計算のベクトル化

Question

私はこの近似テンポラル分数差を高速化しようとしています。 これは、時系列のlong/quasi-longメモリを制御します。最初のforループが反復的であることを考えると、私はそれをベクトル化する方法を知らない。また、試みられたベクトル化の出力は、変更されていない生コードから少し離れている。ご協力ありがとうございました。

生のコード

tempfracdiff= function (x,d,eta) { 

n=length(x);x=x-mean(x);PI=numeric(n) 
PI[1]=-d;TPI=numeric(n);ydiff=x 

for (k in 2:n) {PI[k]=PI[k-1]*(k-1-d)/k} 
for (j in 1:n) {TPI[j]=exp(-eta*j)*PI[j]} 
for (i in 2:n) {ydiff[i]=x[i]+sum(TPI[1:(i-1)]*x[(i-1):1])} 
return(ydiff)     } 

しようとしましたベクトル化

tempfracdiffFL=function (x,d,eta) { 

n=length(x);x=x-mean(x);PI=numeric(n) 
PI[1]=-d;TPI=numeric(n);ydiff=x 

for (k in 2:n) {PI[k]=PI[k-1]*(k-1-d)/k} 
TPI[1:n]=exp(-eta*1:n)*PI[1:n] 
ydiff[2:n]=x[2:n]+sum(TPI[1:(2:n-1)]*x[(2:n-1):1]) 
return(ydiff)   } 

Answer 1

PIについては、使用することができますcumprod：

k <- 1:n 
PI <- cumprod((k-1-d)/k) 

TPIは、インデックスなしで表現することができる。

だから

ydiff <- x+c(0,convolve(x,rev(TPI),type="o")[1:n-1]) 

すべて一緒にそれを置く：

そしてydiffはxプラスxの畳み込みとTPIある

mytempfracdiff = function (x,d,eta) { 
    n <- length(x) 
    x <- x-mean(x) 
    k <- 1:n 
    PI <- cumprod((k-1-d)/k) 
    TPI <- exp(-eta*k)*PI 
    x+c(0,convolve(x,rev(TPI),type="o")[1:n-1]) 
} 

テストケースの例

set.seed(1) 
x <- rnorm(100) 
d <- 0.1 
eta <- 0.5 

all.equal(mytempfracdiff(x,d,eta), tempfracdiff(x,d,eta)) 
# [1] TRUE 

library(microbenchmark) 
microbenchmark(mytempfracdiff(x,d,eta), tempfracdiff(x,d,eta)) 

 
    Unit: microseconds 
          expr  min  lq  mean median  uq 
    mytempfracdiff(x, d, eta) 186.220 198.0025 211.9254 207.473 219.944 
     tempfracdiff(x, d, eta) 961.617 978.5710 1117.8803 1011.257 1061.816 
      max neval 
     302.548 100 
    3556.270 100 

Answer 2

PI [k]についてはReduceが役に立ちます

n <- 5; d <- .3 
fun <- function(a,b) a * (b-1-d)/b 
Reduce(fun, c(1,1:n), accumulate = T)[-1] # Eliminates PI[0] 

[1] -0.30000000 -0.10500000 -0.05950000 -0.04016250 -0.02972025