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