r - 式 y_(t 1)=ay_t + bx_t を使用して、a、b、y_0、および x_t を指定して生成を高速化

Question

R で最速の再帰方程式:

R で、、および を指定y(t+1)=a*y(t)+b*x(t)して、方程式 を使用してデータを生成する最速の方法は何ですか?aby_0x_t

例：

a=3;b=2;y_0=1
x_t=1:10 
y_t=matrix(rep(0,11),ncol=11) 
y_t[1]=y_0 
tt<-2 
while (tt<=11) { 
   y_t[tt]<-a*y_t[tt-1]+b*x_t[tt-1] 
   tt<-tt+1 
} 
y_t 

Answer 1

編集：whileループの出力と比較した後、おそらくこれです：

Reduce(function(X,Y) { b*Y + a*X}, x_t, init=y_0, acc=TRUE)

関数のxとyをReduce正しく取得することはできません。はぁ。

pkg：microbenchmarkとより長いx_tベクトルの使用：

a=3;b=2;y_0=1; x_t=c(1:10, 10:(-10), (-10):10)
> res <- microbenchmark( Rres=invisible( Reduce(function(X,Y) { b*Y + a*X}, x_t, init=y_0, acc=TRUE)) , times=1000L)
> res
Unit: microseconds
  expr     min      lq   median       uq      max
1 Rres 302.114 310.787 316.5705 328.1195 13770.14

> res <- microbenchmark(  Wres= {a=3;b=2;y_0=1; x_t=c(1:10, 10:(-10), (-10):10)
+  y_t=matrix(rep(0,11),ncol=11) 
+  y_t[1]=y_0 
+  tt<-2 
+  while (tt<=53) { y_t[tt]<-a*y_t[tt-1]+b*x_t[tt-1] 
+  tt<-tt+1 }} , times=1000L)
> res
Unit: microseconds
  expr     min       lq   median       uq      max
1 Wres 461.141 470.7545 477.2865 503.7155 13165.52

必ずしも大きな違い（2倍未満）ではありませんが、Reduceより高速です。これは、速度を求めている場合にパッケージ「Rcpp」を調べる必要がある種類の問題です。インラインC++コードは、通常50〜100倍の改善をもたらします。

Answer 2

再帰を解くことで答えを思いつきました：

## y(1)=y0
## y(2) = a*y(1)+b*x(1) = 5
## y(3) = a^2*y(1)+a*b*x(1)+b*x(2) = 9 + 6 + 4
## y(4) = a^3*y(1)+a^2*b*x(1)+a*b*x(2)+b*x(3)

f_recurse <- function(a=3,b=2,y_0=1,x_t=1:10) {
    n <- length(x_t)
    m <- matrix(nrow=n+1,ncol=n)
    m[] <- ifelse(col(m)>=row(m),0,a^(row(m)-col(m)-1))
    a^(0:n)*y_0 + b*rowSums(sweep(m,2,x_t,"*"))
}

恐ろしいことに、@DWinの回答(テスト目的で複製)よりも遅いことが判明しました:

f_reduce <- function(a=3,b=2,y_0=1,x_t=1:10) {
    Reduce(function(X,Y) { b*Y + a*X}, x_t, init=y_0, acc=TRUE)
}


library(rbenchmark)
benchmark(f_recurse(),f_reduce(),replications=2000)
##          test replications elapsed relative 
## 1 f_recurse()         2000   0.430    1.706 
## 2  f_reduce()         2000   0.252    1.000    

多項式を構築するために何か賢いことをすることで、より速くなるかもしれません。また、前の項を計算せずに、結果の n^th 項を直接取得するのが簡単であるという利点もあります...それが役立つ場合...