AW: [R] Full enumeration, how can this be vectorized
Uwe Ligges
ligges at statistik.uni-dortmund.de
Tue Nov 26 11:13:06 CET 2002
Daniel Hoppe wrote:
>
> >
> > Without looking closer, I think this is a very nice example
> > for a piece
> > of code that should be implemented in either C or Fortran (on the one
> > hand it is easy to implement, on the other hand its "S interpreted"
> > execution takes quite a long time). OK, quite a few possible
> > improvements of the R code are obvious, but I don't think this is the
> > best point to start ...
> >
>
> Probably you are right about C and Fortran. But coming from the combination
> of Java and Windows it's hard to overcome one's inhibitions and get started
> with C :-). Beside that I'm specifically trying to improve my R-skills a
> little because I wrote quite a lot of code for my thesis with R and
> therefore I'm curious for better "code patterns". In my R-code I vectorized
> the function fun (thanks Patrick), seems to run somehow faster already.
>
> > Nevertheless, the code looks like you are searching on a grid. Are you
> > really sure whether there is no method in optim() which
> > performs better?
> >
>
> I tried them but was not successful yet. The search range is constrained,
> and Nelder-Mead doesn't use constraints. L-BFGS-B handles the constraints,
> but it isn't able to handle functions with poles. The other ones did not
> qualify for the same reasons.
>
> I wonder if a grid-search is a very unusual case or if there might be a need
> for a fast and generic implementation of such an algorithm.
>
> Daniel
>
Just for fun two steps on the way to faster code:
fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
{
res <- Inf
seqpar <- seq(along = par)
best <- rep(-1, 3)
x <- seq(.01, 2, by = stepSize[1])
y <- seq(.01, 2, by = stepSize[2])
z <- seq(.01, 15, by = stepSize[3])
for (xi in x)
for (yj in y)
for (zk in z){
# Pass the current best result to the function being minimized
to allow for early termination
newRes <- fun(xi, yj, zk, par, seqpar, res)
## UL: I think (newRes < res) is more often FALSE than
!is.na(newRes),
## UL: so turn around these two statements in that case:
if (!is.na(newRes) && newRes < res){
res <- newRes
best <- c(xi, yj, zk)
}
}
return (best)
}
fun <- function(x, y, z, par, seqpar, currentBestValue = Inf)
{
ym <- y - 1
result <- sum((((x*z) / ym * (1 - (z / (z + seqpar))^ym)) - par)^2)
# Allow for short-cut evaluation if running in complete enumeration
mode
if(is.na(result) || currentBestValue < result)
return(Inf)
return(result)
}
[Not tested!!!]
==========================================================
Next step: The following will be much faster, I guess [code not
tested!], but fun() is combined with fillEnumeration, which is not
reusable for arbitrary "funs" know....
fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
{
res <- Inf
seqpar <- seq(along = par)
best <- rep(-1, 3)
x <- seq(.01, 2, by = stepSize[1])
y <- seq(.01, 2, by = stepSize[2])
z <- seq(.01, 15, by = stepSize[3])
for (zk in z){
zterm <- (zk / (zk + seqpar))
for (yj in (y-1)){
yterm <- yj / (1 - zterm^yj)
for (xi in x){
# Pass the current best result to the function being minimized
to allow for early termination
newRes <- sum((((xi*zk) / yterm) - par)^2)
## UL: I think (newRes < res) is more often FALSE than
!is.na(newRes),
## UL: so turn around these two statements in that case:
if (!is.na(newRes) && newRes < res){
res <- newRes
best <- c(xi, yj, zk)
}
}
}
}
return (best)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list