convolution of positively supported functions in R

Question

I want the convolution of two functions defined on [0,Inf), say

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

and

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

Using the integrate function of R I can do this,

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

By searching the web, it seems that there are more efficient (and more accurate) ways of doing this (say using fft() or convolve()). Can anyone with such experiences explain how please? Thanks!

Answer 1

convolve or fft solutions are to get a discrete result, rather than a function as you have defined in cfg. They can give you the numeric solution to cfg on some regular, discrete input.

fft is for periodic functions (only) so that is not going to help. However, convolve has a mode of operation called "open", which emulates the operation that is being performed by cfg.

Note that with type="open", you must reverse the second sequence (see ?convolve, "Details"). You also have to only use the first half of the result. Here is a pictoral example of the result of convolution of c(2,3,5) with c(7,11,13) as would be performed by convolve(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

Note that evaluation the first three elements is similar to the results of your integration. The last three would be used for the reverse convolution.

Here is a comparison with your functions. Your function, vectorized, plotted with

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

And an application of convolve plotted with the following code. Note that there are 100 points per unit interval in y so division by 100 is appropriate.

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

These do not quite agree, but the convolution is much faster:

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