Vektorisierung der temperierten fraktionalen Differenzrechnung

Question

Ich versuche, diese Annäherung der temperierten fraktionalen Differenzierung zu beschleunigen. Dies steuert den langen/quasi-langen Speicher einer Zeitreihe. Da die erste for-Schleife iterativ ist, kann ich sie nicht vektorisieren. Außerdem ist die Ausgabe der versuchten Vektorisierung etwas von dem unveränderten Rohcode entfernt. Danke für Ihre Hilfe.

Rohcode

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)     } 

Versuchte Vektorisierung

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

Für PI, können Sie cumprod:

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

TPI können ohne Indizes ausgedrückt werden:

TPI <- exp(-eta*k)*PI 

Und ydiff ist x plus die Faltung von x und TPI:

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

So setzen sie alle zusammen:

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]) 
} 

Testfallbeispiel

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

Für PI [k], Reduce ist hilfreich

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