[R] Full enumeration, how can this be vectorized

[Daniel Hoppe]

Hi all,

I am currently using the code snippet below to minimize a function which I
wasn't able to minimize properly with optim (cliffy surface, constraints,
poles). Obviously this iteration is probably the poorest way of implementing
such a function and it is expectedly slow. I've already written some
vectorized functions, but I'm afraid that I'm missing the "big picture" in
this case, I don't know how to start. Probably I could build a 3 x
(length(x) * length(y) * length(z)) matrix holding each possible parameter
combinations. With my current step-size and value range that would be
slightly more than 2,000,000 triples, so that should not be a problem. Then
comes the part I'm totally clueless about, that is applying the target
function to each of the triples in the matrix in an efficient way. From the
resulting vector I could then easily take the minimum value index and get
the parameters from the parameter matrix.

Could someone kindly give me a hint how this could be implemented?

Thanks a lot,

Daniel


Code Snippet:
fullEnumeration <- function(par, stepSize = c(.03, .03, .03))
{
   res <- Inf

   best.x <- -1
   best.y <- -1
   best.z <- -1

   stepx <- stepSize[1]
   stepy <- stepSize[2]
   stepz <- stepSize[3]
   x <- seq(.01, 2, by=stepx)
   y <- seq(.01, 2, by=stepy)
   z <- seq(.01, 15, by=stepz)

   for (i in 1:length(x))
    for (j in 1:length(y))
        for (k in 1:length(z))
        {
		# Pass the current best result to the function being minimized to allow
for early termination
            newRes <- fun(c(x[i], y[j], z[k]), par, res)
            if (!is.na(newRes) && newRes < res)
            {
                res = newRes
                best.x <- x[i]
                best.y <- y[j]
                best.z <- z[k]
            }

        }
   return (c(
        best.x,
        best.y,
        best.z))

}

fun <- function(x, par, currentBestValue=Inf)
{
  x = x[1]
  y = x[2]
  z = x[3]
  TMax <- length(par)

  result <- 0
  for (myT in 1:TMax)
  {
    result <- result +
    (
    ((x * z) / (y-1) * (1 - (z / (z + myT))^(y-1) ))
    - par[myT]
    )^2

    # Allow for short-cut evaluation if running in complete enumeration mode

    if (is.na(result) || currentBestValue < result)
    {
        return (Inf)
    }

  }
  return (result)

}



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