Nothing
## R implementation of subplex. This exists only for reference, and
## is not exported. Some of the auxilliary functions (partition,
## simplex) are used by a real/integer mixed optimisation algorithm
## (mle-mixed.R)
do.mle.search.subplexR <- function(func, x.init, control, lower,
upper) {
control <- modifyList(list(reltol=.Machine$double.eps^0.25,
parscale=rep(.1, length(x.init))),
control)
check.bounds(lower, upper, x.init)
if ( any(is.finite(lower) | is.finite(upper)) )
func2 <- invert(boxconstrain(func, lower, upper))
else
func2 <- invert(func)
ans <- subplexR(x.init, func2, control)
ans$value <- -ans$value
names(ans)[names(ans) == "value"] <- "lnLik"
if ( ans$convergence != 0 )
warning("Convergence problems in find.mle (subplex): ",
tolower(ans$message))
ans
}
subplexR <- function(par, fn, control=list(), hessian=FALSE, ...) {
sign2 <- function(x) (-2*(x <= 0) + 1)
constrain2 <- function(f, p, i) {
function(x) {
p[i] <- x
f(p)
}
}
if ( hessian )
.NotYetUsed("hessian")
nx <- length(par) # dimension
if ( nx < 2 )
stop("Need at least two dimensions for subplexR")
control <- modifyList(list(maxit=10000, parscale=rep(1, nx),
reltol=.Machine$double.eps), control)
max.eval <- control$maxit
scale <- control$parscale
tolerance <- control$reltol
func <- function(x) fn(x, ...)
psi <- 0.25
omega <- .1
if ( length(scale) == 1 )
scale <- rep(scale, nx)
else if ( length(scale) != nx )
stop("Invalid scale argument")
subspc.min <- min(2, nx)
subspc.max <- min(5, nx)
cur.x <- par
cur.y <- func(par)
ne <- 1
delta.x <- step.size <- scale
flag <- 0
rep <- 0
repeat {
rep <- rep + 1
prev.x <- cur.x
delta.x <- abs(delta.x)
order.x <- order(delta.x, decreasing=TRUE)
dims <- partition2(nx, order.x, delta.x, subspc.min, subspc.max)
## Simplex loop.
for ( i in seq_along(dims) ) {
j <- dims[[i]]
res <- tryCatch(simplex(constrain2(func, cur.x, j), cur.x[j],
step.size[j], max.eval-ne, y.init=cur.y,
dx=control$dx[j]),
simplexPrecisionError=function(e) e$ret)
cur.x[j] <- res$par
cur.y <- res$val
ne <- ne + res$ne
flag <- res$flag
if ( flag != 0 )
break # will return
}
if ( flag != 0 )
break # will return
delta.x <- cur.x - prev.x
done <- any(pmax(abs(delta.x), abs(step.size)*psi)/
pmax(abs(cur.x),1) <= tolerance)
if ( done ) { ## Tolerance is satisfied
break # will return
} else {
## in this case, we rescale a little, then go back simplex
if ( length(dims) > 1 )
step.factor <- min(max(sum(abs(delta.x))/sum(abs(step.size)),
omega), 1/omega)
else
step.factor <- psi
step.size <- abs(step.size * step.factor) * sign2(delta.x)
}
}
msg <- list("number of function evaluations exceeds `maxit'",
NULL,
"limit of machine precision reached",
"fstop reached")[[flag+2]]
list(par=cur.x, value=cur.y, count=ne, convergence=flag, message=msg)
}
## Simplified simplex, following the code directly.
simplex <- function(func, x.init, scale, max.eval=1000, alpha=1.0,
beta=0.5, gamma=2.0, delta=0.5, psi=0.25,
y.init=NULL, dx=NULL) {
simplexPrecisionError <- function(x, y) {
e <- simpleError("Limit of machine precision reached")
e$ret <- list(par=x, val=y, flag=1, ne=ne)
class(e) <- c("simplexPrecisionError", class(e))
stop(e)
}
newpt <- function(coef, x.base, x.old) {
x.new <- x.base + coef * (x.base - x.old)
## Enforces grid.
x.new[on.grid] <-
round_to_grid(x.new[on.grid], x.init[on.grid], dx[on.grid])
if ( identical(x.new, x.base) || identical(x.new, x.old) )
simplexPrecisionError(xx[ind.lo,], yy[ind.lo])
x.new
}
nx <- length(x.init)
ne <- 0
if ( is.null(dx) ) {
on.grid <- rep(FALSE, nx)
} else {
on.grid <- !is.na(dx) & dx > 0
scale[on.grid] <- pmax(round_to_grid(scale[on.grid], 0,
dx[on.grid]), dx[on.grid])
}
if ( length(scale) != nx )
stop("Invalid length scale")
if ( is.null(y.init) ) {
y.init <- func(x.init)
ne <- 1
}
xx <- rbind(x.init, t(x.init + diag(nx) * scale), deparse.level=0)
xl <- matrix.to.list(xx)
if ( any(unlist(lapply(xl[-1], identical, x.init))) )
simplexPrecisionError(x.init, y.init)
yy <- c(y.init, unlist(lapply(xl[-1], func)))
ne <- ne + nx
ord <- order(yy)
ind.lo <- ord[1]
ind.sec <- ord[nx]
ind.hi <- ord[nx+1]
tolerance <- psi^2 * sum((xx[ind.hi,]-xx[ind.lo,])^2)
repeat {
centroid <- colMeans(xx[-ind.hi,])
## "Reflection"
x.ref <- newpt(alpha, centroid, xx[ind.hi,])
y.ref <- func(x.ref)
ne <- ne + 1
if ( y.ref < yy[ind.lo] ) { # new best - keep pushing
## "Expansion"
x.exp <- newpt(-gamma, centroid, x.ref)
y.exp <- func(x.exp)
ne <- ne + 1
if ( y.exp < y.ref ) {
xx[ind.hi,] <- x.exp
yy[ind.hi] <- y.exp
} else {
xx[ind.hi,] <- x.ref
yy[ind.hi] <- y.ref
}
} else if ( y.ref < yy[ind.sec] ) {
## Neither the best nor the worst (in the remaining simplex):
## accept current point
xx[ind.hi,] <- x.ref
yy[ind.hi] <- y.ref
} else { ## Would be the current worst point.
## "Contraction"
x.base <- if ( yy[ind.hi] < y.ref ) xx[ind.hi,] else x.ref
x.con <- newpt(-beta, centroid, x.base)
y.con <- func(x.con)
ne <- ne + 1
if ( y.con < min(y.ref, yy[ind.hi]) ) {
## Successful contraction!
xx[ind.hi,] <- x.con
yy[ind.hi] <- y.con
} else {
## "Massive Contraction"
for ( i in seq_len(nx + 1)[-ind.lo] ) {
xx[i,] <- newpt(-delta, xx[ind.lo,], xx[i,])
yy[i] <- func(xx[i,])
}
ne <- ne + nx
}
}
ord <- order(yy)
ind.lo <- ord[1]
ind.sec <- ord[nx]
ind.hi <- ord[nx+1]
if ( ne > max.eval ) {
flag <- -1
break
} else if ( sum((xx[ind.hi,]-xx[ind.lo,])^2) <= tolerance ) {
flag <- 0
break
}
}
list(par=xx[ind.lo,], val=yy[ind.lo], flag=flag, ne=ne)
}
## This is a direct copy of the Matlab partition.m file's function
## partitionx().
## Description of what is going on:
## "Partition the space defined by deltax[order] into a group of
## sets, such that subs_max <= k <= subs_max for all sets of size k.
## Where possible, pick combinations of points which maximize a-b,
## where a is the norm of the new set /k and b is the norm of the
## remainder over (nx-k)."
partition <- function(nx, order, dx, subs.min, subs.max) {
num.subspaces <- 0
nused <- 0
nleft <- nx
asleft <- sum(dx)
subsp.dims <- rep(NA, nx)
while ( nused < nx ) {
num.subspaces <- num.subspaces + 1
as1 <- sum(dx[order[(nused+1):(nused+subs.min-1)]])
gapmax <- -1
for ( ns1 in subs.min:(min(subs.max,nleft)) ) {
as1 <- as1 + dx[order[nused+ns1]]
ns2 <- nleft - ns1
if ( ns2 > 0 ) {
if ( ns2 >= (trunc((ns2-1)/subs.max+1)*subs.min) ) {
as2 <- asleft - as1
gap <- as1/ns1 - as2/ns2
if ( gap > gapmax ) {
gapmax <- gap
subsp.dims[num.subspaces] <- ns1
as1max <- as1
}
}
} else {
if ( as1/ns1 > gapmax ) {
subsp.dims[num.subspaces] <- ns1
return(list(num.subspaces, subsp.dims[1:num.subspaces]))
}
}
}
nused <- nused + subsp.dims[num.subspaces]
nleft <- nx - nused
asleft <- asleft - as1max
}
list(num.subspaces, subsp.dims[1:num.subspaces])
}
## This just wraps the above partition() function in a more R-ish way
## (at the cost of some speed due to the call to split()), which makes
## the book-keeping easier.
partition2 <- function(nx, order.x, delta.x, subspc.min, subspc.max) {
tmp <- partition(nx, order.x, delta.x, subspc.min, subspc.max)
ret <- split(order.x[seq_len(nx)], rep(seq_len(tmp[[1]]), tmp[[2]]))
names(ret) <- NULL
ret
}
## xmid may be any point that is on the grid.
round_to_grid <- function(x, xmid, dx)
round((x - xmid)/dx)*dx + xmid
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.