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)、a、和给定的生成数据的最快方法是什么？by_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

我通过解决递归想出了一个答案：

## 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 项，而无需计算任何先验项……如果有用的话……

Answer 2

编辑：与您的 while 循环输出进行比较后，可能是这样的：

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

我永远无法Reduce正确获得函数中的 x 和 y。叹。

使用 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 倍的改进。