R中循环的替代

Mor*_*ten 3 for-loop r

可能重复:
加速R中的循环操作

我有几个关于循环的问题.我知道R使用矢量化计算可以更快地工作,我想更改下面的代码来利用这一点.在论坛上查看其他一些答案,sapply函数似乎能够替换inside for循环,但我生成一个零向量,所以有一个错误.陶仍然是1000,我认为这是造成问题的原因.

我主要关心的是速度,因为我需要围绕整个算法创建一个循环,并绘制不同的V和n尺寸以进行进一步分析.

谢谢你的帮助

替代循环

tao = 1000
L = (tao - 1)   
n = 10      
V = 5               
I = 10000                       
V_s = matrix(rnorm(I), I, 1)
V_b = matrix(rnorm(I), I, 1)

signal <- matrix(0, L, 1)  

for( j in (n:L)){

    sapply(((j-n+1):j),function (tao) signal[j] = signal[j] + abs(V_s[tao] - V_b[tao]))

    signal[j] = (signal[j] / (n * V) )

} 
Run Code Online (Sandbox Code Playgroud)

原始循环

tao = 1000
L = (tao - 1)   
n = 10      
V = 5               
I = 10000                       
V_s = matrix(rnorm(I), I, 1)
V_b = matrix(rnorm(I), I, 1)

signal <- matrix(0, L, 1)  

for( j in (n:L)){

    for( tao in ((j-n+1):j))    {

        signal[j] = (signal[j] + abs(V_s[tao] - V_b[tao]))

    }
        signal[j] = (signal[j] / (n * V) )

}
Run Code Online (Sandbox Code Playgroud)

Thi*_*ilo 12

使用过滤器,即使没有任何循环也可以进行计算(并且sapply只不过是一个隐藏循环).

absdif <- abs(V_s - V_b)
signal <- filter(absdif[1:L], rep(1/(n*V), n), sides=1)
signal[is.na(signal)] <- 0
Run Code Online (Sandbox Code Playgroud)

但是,当您不习惯过滤时,了解第二行中发生的事情并非易事.让我们仔细看看:

首先,我们计算循环经常使用的V_s和的绝对差异V_b.然后过滤器.您的计算只不过是n在每个时间值总结过去的值j.因此,我们有类似的东西

signal[j] <- sum(absdif[j-n+1:j])
Run Code Online (Sandbox Code Playgroud)

这正是卷积滤波器所做的 - 总结一些值 - 通过乘以一些权重的一般形式.1/(n*V)对于所有值,我们选择的权重与您在外部循环中执行的规范化相对应.最后一个参数,sides=1只是告诉过滤器只从过去获取值(sides=2意味着sum(absdif[(j-n/2):(j+n/2)])).

最后一行只是填充了NA开头的值(过滤器没有足够的数据来计算总和 - 这等于跳过第一个n值).

最后,一些时间:

您的完整循环解决方案:

   User      System       total 
  0.037       0.000       0.037 
Run Code Online (Sandbox Code Playgroud)

朱巴的解决方案:

   User      System       total 
  0.007       0.000       0.008 
Run Code Online (Sandbox Code Playgroud)

使用过滤器的解决方案

   User      System       total 
  0.000       0.000       0.001 
Run Code Online (Sandbox Code Playgroud)

请注意,过滤器的概念非常好,可以非常快速地完成.

编辑:如上所述?filter,R不使用标准filter命令的快速傅里叶变换.通常,FFT是实现卷积的最有效方式.但是,即使这样也可以通过替换filter命令来完成

signal <- convolve(absdif[1:L], rep(1/(n*V), n), type='filter')
Run Code Online (Sandbox Code Playgroud)

请注意,现在第一个n条目被剥离而不是设置为NA.但结果是一样的.这次的时间没有用 - 总时间低于system.time......的三位数输出.但是,请注意R帮助中的以下注释filter:

convolve(,type ="filter")使用FFT进行计算,因此对于单变量序列上的长过滤器可能更快,但它不返回时间序列(因此时间对齐不清楚),也不处理缺失值.例如,对于长度为1000的系列的长度为100的滤波器,滤波器更快