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
****************************************************************************
**********
|