[R] Computation of contour values - Speeding up computation

[Andreas Wittmann]

Thank you very much, yes maybe its worth working on it.

best regards

Andreas




Uwe Ligges wrote:
>
>
> 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.
>