[R] Computation of contour values - Speeding up computation

[Uwe Ligges]

Andreas Wittmann wrote:
> Dear R useRs, 
> 
> i have the following code to compute values needed for a contour plot
> 
> ############################################################
> 
> "myContour" <- function(a, b, plist, veca, vecb, dim)
> {
>   tmpb <- seq(0.5 * b, 1.5 * b, length=dim)
>   tmpa  <- seq(0.5 * a, 1.5 * a, length=dim)
> 
>   z <- matrix(0, nrow=dim, ncol=dim)
>   for(i in 1:dim)
>   {
>     for(j in 1:dim)
>     {
>       z[i, j] <- posteriorPdf(a=tmpa[j], b=tmpb[i], 
>                               plist=plist, veca=veca, vecb=vecb) 
>     }
>   }
> }
> 
> "posteriorPdf" <- function(a, b, plist, veca, vecb)
> {
>   res <- sum(plist[, 1] *       
>              exp(vecb[, 1] * log(vecb[, 2]) + (vecb[, 1] - 1.0) * log(b) 
>                  - vecb[, 2] * b - lgamma(vecb[, 1])) *
>              exp(veca[, 1] * log(veca[, 2]) + (veca[, 1] - 1.0) * log(a) 
>                  - veca[, 2] * a - lgamma(veca[, 1])))
>           
>   return(res)          
> }
> 
> plist <- matrix(0, 100, 3)
> plist[, 1] <- runif(100) 
> veca <- vecb <- matrix(0, 100, 2)
> 
> veca[, 1] <- seq(20, 50, len=100)
> veca[, 2] <- seq(10, 20, len=100)
> 
> vecb[, 1] <- seq(50, 200, len=100)
> vecb[, 2] <- seq(1000, 400000, len=100)
> 
> myContour(a=20, b=0.01, plist=plist, veca=veca, vecb=vecb, dim=50)
> 
> ############################################################
> 
> this is part of my other computations which i do with R. Here i recognized, that my functions myContour and posteriorPdf took a long time of my computations. The key to speed this up is to avoid the two for-loops in myContour, i think. I tried a lot to do this with apply or something like that, but i didn't get it. 
> 
> If you have any advice how i can to this computations fast, i would be very thankful, one idea is to use external c-code?
> 


It takes 0.8 seconds on my machine. Not worth working on it, is it?

Your problem is that you are applying many calculations for all 
iterations of the inner loop, even if the result won't change, example:
lgamma(veca[, 1]) will be calculated dim^2 times!

Hence you can improve your loop considerably:



myContour <- function(a, b, plist, veca, vecb, dim)
{
   tmpb <- seq(0.5 * b, 1.5 * b, length=dim)
   tmpa  <- seq(0.5 * a, 1.5 * a, length=dim)

   z <- matrix(0, nrow=dim, ncol=dim)
   plist1 <- plist[, 1]
   vecb1l2 <- vecb[, 1] * log(vecb[, 2])
   vecb11 <- vecb[, 1] - 1
   vecb1lg <- lgamma(vecb[, 1])
   vecb2 <- vecb[, 2]
   veca1l2 <- veca[, 1] * log(veca[, 2])
   veca11 <- veca[, 1] - 1
   veca2 <- veca[, 2]
   veca1lg <- lgamma(veca[, 1])

   for(i in 1:dim)
   {
     for(j in 1:dim)
     {
       z[i, j] <- sum(plist1 * exp(vecb1l2 + vecb11 * log(tmpb[i]) - 
vecb2 * tmpb[i] - vecb1lg) * exp(veca1l2 + veca11 * log(tmpa[j]) - veca2 
* tmpa[j] - veca1lg))
     }
   }
   z
}

Uwe Ligges






> best regards
> Andreas 
> --
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.