R中积极支持的函数的卷积

2023年4月11日 195次阅读

我想要在[0,Inf)上定义的两个函数的卷积,比方说

f=function(x)
    (1+0.5*cos(2*pi*x))*(x>=0)

和

g=function(x)
    exp(-2*x)*(x>0)

使用R的集成功能可以做到这一点,

cfg=function(x)
    integrate(function(y) f(y)*g(x-y),0,x)$value

通过搜索网络,似乎有更有效(更准确)的方法(比如使用fft()或convolve()).任何有这种经历的人都可以解释一下吗？
谢谢！

最佳答案 convolve或fft解决方案是获得离散结果,而不是你在cfg中定义的函数.他们可以在一些常规的离散输入上为您提供cfg的数值解决方案.

fft用于定期功能(仅限),因此无法提供帮助.但是,convolve有一种称为“open”的操作模式,它模拟cfg正在执行的操作.

请注意,使用type =“open”,您必须反转第二个序列(请参阅？convolve,“Details”).您还必须仅使用结果的前半部分.下面是c(2,3,5)与c(7,11,13)卷积结果的图像示例,如卷积(c(2,3,5),rev(c(7, 11,13)),type =’open’)：

             2  3  5       2  3  5   2  3  5    2  3  5      2  3  5
      13 11  7         13 11  7     13 11  7      13 11  7        13  11  7
 Sum:       14              43         94           94            65

请注意,评估前三个元素与您的集成结果类似.最后三个将用于反向卷积.

这是与您的功能的比较.你的功能,矢量化,绘制

y <- seq(0,10,by=.01)
plot(y, Vectorize(cfg)(y), type='l')

并使用以下代码绘制卷积的应用程序.请注意,y中每单位间隔有100个点,因此除以100是合适的.

plot(y, convolve(f(y), rev(g(y)), type='open')[1:1001]/100, type='l')

这些并不完全一致,但卷积要快得多：

max(abs(Vectorize(cfg)(y) - convolve(f(y), rev(g(y)), type='open')[1:1001]/100))
## [1] 0.007474999

benchmark(Vectorize(cfg)(y), convolve(f(y), rev(g(y)), type='open')[1:1001]/100, columns=c('test', 'elapsed', 'relative'))
##                                                   test elapsed relative
## 2 convolve(f(y), rev(g(y)), type = "open")[1:1001]/100   0.056        1
## 1                                    Vectorize(cfg)(y)   5.824      104
点赞