CODEWITHSUNDEEP
r
Joseph Gardiner
Joseph Gardiner

Reputation: 107

More efficient matrix operation in R

I am trying to compute the convolution of two discrete probability distributions in R. I have two vectors, each containing the probabilities. I need to compute a third vector, which has the combined probabilities of the previous two vectors. The ith position in the third vector contains the sum of probabities (a[j]*b[k]) for all j+k=i. I have the following function for doing this: 

convolute <- function(a, b ){

    out <- rep(0, (length(a)+741))

    for(i in 1:length(a)){

        for (j in 1:length(b)){

            out[i + j] <- out[i+j] + (a[i]*b[j])


        }
    }

    return(out)
}

My problem is that this function needs to be called multiple times (>1000000) and is (relatively) slow. Is there an more efficient way in R to achieve this operation, without using the two for loops? The length of a will either be 741 or 1482, b is always 741.

Thank you

Upvotes: 1

Views: 134

Answers (1)

Jan van der Laan
Jan van der Laan

Reputation: 8105

convolve(a, rev(b), type="open")

Does the same as your function, except that your function starts with a 0 and convolve doesn't:

> a <- runif(1000, 0, 1)
> b <- runif(741, 0, 1)
> c1 <- convolute(a, b)
> c2 <- convolve(a, rev(b), type="open")
> 
> all.equal(c1[-1], c2)
[1] TRUE
> system.time(c1 <- convolute(a, b))
   user  system elapsed 
  4.152   0.000   4.155 
> system.time(c2 <- convolve(a, rev(b), type="open"))
   user  system elapsed 
  0.000   0.000   0.001

Upvotes: 4

Related Questions