可能重复:
加速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的滤波器,滤波器更快