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