Computation of contour values - Speeding up computation

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 --