If you use square the deviation rather than using absolute value, you
might get better numerical performance, e.g., from an algorithm that
assumes the function is differentiable at the optimum!
Spencer Graves
Liaw, Andy wrote:
[cc to s-news as I supposed is intended.]
Besides, why square beta(x, y) - 0.5 and then take square root? Why not
just use abs()?
Andy
-----Original Message-----
From: Spencer Graves [mailto:spencer.graves@PDF.COM]
Sent: Wednesday, February 19, 2003 11:39 AM
To: Liaw, Andy
Subject: Re: [S] beta function,numerical methods,amoeba
I agree with Andy Liaw that it is better to write the beta
function in
terms of lgamma. I suggest you try carrying this one step further:
Cast your equation for the beta function into one for the log beta
function, e.g., using
lbeta <- function(a, b) lgamma(a) + lgamma(b) -lgamma(a+b)
Beyond this, numerical methods often work best if you provide
reasonable
starting values. Reasonably starting values can often be
obtained using
Stirling's approximation for the log of the gamma function:
lgamma(a) = (a-.5)*ln(a)-a+.5*ln(2*pi)+O(1/a),
where the error O(z) is bounded by a constant time z as z -> Inf.
If you ignore the error and play with this, you might be able
to find a
fairly simple formula that will produce a reasonable first guess.
Spencer Graves
Liaw, Andy wrote:
See if this helps (optim() is in the most recent version of MASS):
beta <- function(a, b) {
+ exp( lgamma(a) + lgamma(b) - lgamma(a+b) )
+ }
beta.f <- function(coef){
+ x<-coef[1]
+ y<-coef[2]
+ sqrt((beta(x,y)-0.5)^2)
+ }
optim(c(3,2), beta.f)
$par:
[1] 3.1686527 0.7762833
$value:
[1] 6.884248e-009
$counts:
function gradient
89 NA
$convergence:
[1] 0
$message:
NULL
The Nelder-Mead simplex (default for optim) was faster than
I can blink.
Andy
-----Original Message-----
From: Steve Su [mailto:s.su@qut.edu.au]
Sent: Wednesday, February 19, 2003 2:40 AM
To: s-news@lists.biostat.wustl.edu
Subject: [S] beta function,numerical methods,amoeba
Dear All,
I am having great difficulties with using numerical methods to solve
equations that involve beta function. For illustration
purposes I will show
only the beta function part of the equation.
Currently S+ NuOpt does not support integrate nor gamma
functions so it does
not seem to be possible to solve this through S+ NuOpt.
On the other hand I found there is Nelder Mead simplex method
(included in
this e mail and originally developed by Bill Clark and
modified by others in
the process) which appears to work with beta function. But it
only works if
you define the beta function through an integral (does not
work well if you
define beta function as gamma functions). It works as follows:
# Define beta function through an integral:
beta.function<-function(a,b)
{ integrate(function(x, a, b)
{
(x^(a - 1)) * (1 - x)^(b - 1)
}
, lower = 0, upper = 1, a =a, b =b, subdivisions = 1000)$integral
}
# Define function to be minimized:
beta.f<-function(coef){
x<-coef[1]
y<-coef[2]
sqrt((beta.function(x,y)-0.5)^2)
}
# Using Numerical Method:
neldmead(beta.f,c(3,2))
This method works as long as beta.function is not singular
but it is very
slow. I have looked around the web and it appears the
numerical recipes in
C++ appears to be the alternative solution. However, I am not
sure exactly
how to interface this with Splus as I have no knowledge of
C++ programming.
Can anyone advise me how I can get this into Splus to work so I can
incorporate any function I might need to evaluate? It would
also help if a
Nelder Simplex method with constraints version is available
so I can avoid
the function being stopped because the beta function is singular?
Thanks for your help.
########################
neldmead<-
function(fobj, startvec, print.level = 0, tol = 1e-006, ...)
{
# Standard tuning parameters.
alpha <- 1.
# alpha >= 1
beta <- 0.5
# 0 < beta < 1
gamma <- 2.
# gamma > alpha
# Build initial simplex by scaling single elements of
starting vector by
# a factor exp(+/-lambda).
if(print.level == 2) cat("Initializing simplex...\n")
lambda <- 2
npar <- length(startvec)
nverts <- npar + 1
verts <- matrix(startvec, nverts, npar, byrow = T)
fvals <- rep(0, nverts)
for(iv in 1:npar)
repeat {
lamb <- lambda
{
verts[iv, iv] <- exp(lamb) * startvec[iv]
if(verts[iv, iv] == 0)
verts[iv, iv] <- 1
fvals[iv] <- fobj(verts[iv, ], ...)
if(is.finite(fvals[iv]))
break
verts[iv, iv] <- exp( - lamb) * verts[iv, iv]
fvals[iv] <- fobj(verts[iv], ...)
if(is.finite(fvals[iv]))
break
lamb <- lamb/2
}
}
verts[nverts, ] <- startvec
fvals[nverts] <- fobj(startvec, ...)
scatter.init <- apply(verts, 2, stdev)
centroid <- startvec
# Start iteration.
iter <- 0
repeat {
# Check for convergence.
if(print.level == 2) {
cat("Press <return> to continue, q<return> to quit\n")
out <- readline()
if(out == "q")
break
}
iter <- iter + 1
scatter <- apply(verts, 2, stdev)
conv.test <- max(scatter/scatter.init)
fmax <- max(fvals)
iv.fmax <- min((1:nverts)[fvals == fmax])
fnext <- max(fvals[(1:nverts) != iv.fmax])
iv.fnext <- min((1:nverts)[fvals == fnext])
fmin <- min(fvals)
iv.fmin <- min((1:nverts)[fvals == fmin])
if(conv.test < tol)
break
if(print.level > 0) {
fcen <- fobj(centroid, ...)
cat("Iteration", iter, "... convergence test =",
signif(conv.test, 3),
"fval =", signif(fcen, 8), "\n")
}
# Make a move. First try projecting across centroid from
high point.
# The centroid referred to here is the centroid of the face
# across from the point where the function is largest.
centroid <- if(npar == 1) verts[iv.fmin] else
apply(verts[(1:nverts) !=
iv.fmax, ], 2, mean)
proj.vert <- centroid + alpha * (centroid - verts[iv.fmax, ])
proj.fval <- fobj(proj.vert, ...)
if(print.level == 2) {
cat("Vertices:\n")
for(iv in 1:nverts)
cat(iv, signif(verts[iv, ]), "\n")
cat("Function values:\n", signif(fvals), "\n")
cat("Projected vertex:\n", signif(proj.vert), "\n")
cat("Projected function value:", signif(proj.fval), "\n")
}
# Make sure function values are different.
if(all(fvals == fmax)) {
cat(paste("nldrmd error: Function values same at all
vertices.", "Change
startvec.\n"))
return(NA)
}
# Decide what to do based on projected function value.
if(proj.fval >= fmin & proj.fval <= fnext) {
if(print.level == 2)
cat("Taking projection.\n")
verts[iv.fmax, ] <- proj.vert
fvals[iv.fmax] <- proj.fval
next
}
if(proj.fval < fmin) {
if(print.level == 2)
cat("Performing expansion.\n")
#### HERE IS THE PROBLEM>>> its fixed, ... use proj.vert instead
exp.vert <- centroid - gamma * (centroid - proj.vert)
exp.fval <- fobj(exp.vert, ...)
if(print.level == 2) {
cat("Expansion vertex:\n", signif(exp.vert), "\n")
cat("Expansion function value:", signif(exp.fval), "\n")
}
if(exp.fval < fmin) {
verts[iv.fmax, ] <- exp.vert
fvals[iv.fmax] <- exp.fval
}
else {
verts[iv.fmax, ] <- proj.vert
fvals[iv.fmax] <- proj.fval
}
next
}
if(proj.fval > fnext) {
if(proj.fval < fvals[iv.fmax]) {
verts[iv.fmax, ] <- proj.vert
fvals[iv.fmax] <- proj.fval
}
cont.vert <- (1 - beta) * centroid + (beta) * verts[iv.fmax, ]
cont.fval <- fobj(cont.vert, ...)
if(cont.fval < fmax) {
if(print.level == 2)
cat("Performing 1D contraction.\n")
if(print.level == 2) {
cat("Contraction vertex:\n", signif(cont.vert), "\n")
cat("Contraction function value:", signif(cont.fval), "\n")
}
verts[iv.fmax, ] <- cont.vert
fvals[iv.fmax] <- cont.fval
next
}
# Didn't work. Shrink toward best point.
if(print.level == 2) cat("Performing wholesale contraction.\n")
for(iv in 1:nverts) {
verts[iv, ] <- beta * verts[iv.fmin, ] + (1 - beta) *
verts[iv, ]
fvals[iv] <- fobj(verts[iv, ], ...)
}
}
}
# Convergence test succeeded.
return(x = verts[iv.fmin, ], fmin = min(fvals), iter = iter - 1)
}
**************************************************************
**************
**********
Steve Su (s.su@qut.edu.au)
PhD student.
School of Accountancy
School of Mathematical Sciences
Queensland University of Technology
Postal Address: Steve Su, School of Accountancy, QUT, PO Box 2434,
Brisbane,
Queensland, Australia, 4000.
Phone: +61 7 3864 2017
Fax: +61 7 3864 1812
Mobile: 0421 840 586
.
_--_|\
/ QUT
\_.--._/
v
**************************************************************
**************
**********
--------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu. To
unsubscribe send e-mail to
s-news-request@lists.biostat.wustl.edu with
the BODY of the message: unsubscribe s-news
--------------------------------------------------------------
----------------
Notice: This e-mail message, together with any
attachments, contains information of Merck & Co., Inc.
(Whitehouse Station, New Jersey, USA) that may be
confidential, proprietary copyrighted and/or legally
privileged, and is intended solely for the use of the
individual or entity named in this message. If you are not
the intended recipient, and have received this message in
error, please immediately return this by e-mail and then delete it.
==============================================================
================
--------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu. To
unsubscribe send e-mail to
s-news-request@lists.biostat.wustl.edu with
the BODY of the message: unsubscribe s-news
------------------------------------------------------------------------------
Notice: This e-mail message, together with any attachments, contains information of
Merck & Co., Inc. (Whitehouse Station, New Jersey, USA) that may be
confidential, proprietary copyrighted and/or legally privileged, and is intended
solely for the use of the individual or entity named on this message. If you are
not the intended recipient, and have received this message in error, please
immediately return this by e-mail and then delete it.
==============================================================================
--------------------------------------------------------------------
This message was distributed by s-news@lists.biostat.wustl.edu. To
unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
the BODY of the message: unsubscribe s-news
|